Skip to content
Snippets Groups Projects
Commit a24c2a85 authored by Jiří Kalvoda's avatar Jiří Kalvoda
Browse files

algo/maxcut_naive: Požití Goemans-Williamsonova algoritmu

parent 14f75a55
No related branches found
No related tags found
No related merge requests found
The MaxCut.jl package is licensed under the MIT "Expat" License:
> Copyright (c) 2021: Eric Proffitt.
>
> Permission is hereby granted, free of charge, to any person obtaining
> a copy of this software and associated documentation files (the
> "Software"), to deal in the Software without restriction, including
> without limitation the rights to use, copy, modify, merge, publish,
> distribute, sublicense, and/or sell copies of the Software, and to
> permit persons to whom the Software is furnished to do so, subject to
> the following conditions:
>
> The above copyright notice and this permission notice shall be
> included in all copies or substantial portions of the Software.
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
> EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
> MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
> IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
> CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
> TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
> SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
This software uses supporting software with their own licenses:
https://github.com/JuliaLang/julia/blob/master/LICENSE.md<br />
https://github.com/jump-dev/Convex.jl/blob/master/LICENSE.md<br />
https://github.com/jump-dev/SCS.jl/blob/master/LICENSE.md
name = "MaxCut"
uuid = "6af15267-3226-46fb-9d72-8bee1ad40623"
version = "1.0.1"
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
[compat]
julia = "1"
Convex = "~0.14.16"
**The Goemans-Williamson algorithm for the MAXCUT problem.**
Partition a graph into two disjoint sets such that the sum of the edge weights which cross the partition is as large as possible (known to be NP-hard).
A cut of a graph can be produced by assigning either 1 or -1 to each vertex. The Goemans-Williamson algorithm relaxes this binary constraint to allow for vector assignments drawn from the (n-1)-sphere (choosing an n-1 dimensional space will ensure seperability). This relaxation can then be written as a semidefinite program. Once the optimal vector assignments are found, origin centered hyperplanes are generated and their corresponding cuts evaluated. After `iter` trials, or when the desired tolerance is reached, the hyperplane with the highest corresponding binary cut is used to partition the vertices.
### Installation
```julia
(@v1.7) pkg> add https://github.com/ericproffitt/MaxCut.jl
```
### Dependencies
```julia
LinearAlgebra
Convex
SCS
```
### Arguments
```julia
W: Adjacency matrix.
tol: (keyword arg) Maximum acceptable distance between a cut and the MAXCUT upper bound (default=0).
iter: (keyword arg) Maximum number of hyperplane iterations before a cut is chosen (default=100).
```
### Example
```julia
using MaxCut
W = [0 5 2 1 0;
5 0 3 2 0;
2 3 0 0 0;
1 2 0 0 4;
0 0 0 4 0];
max_cut, max_partition = maxcut(W);
@show max_cut;
## max_cut = 14
@show max_partition;
## max_partition = ([1, 3, 4], [2, 5])
```
### Reference
http://www.sfu.ca/~mdevos/notes/semidef/GW.pdf
\ No newline at end of file
module MaxCut
using LinearAlgebra
using Convex
using SCS
"The Goemans-Williamson algorithm for the MAXCUT problem."
function maxcut(W::Matrix{<:Real}; iter::Int=100, tol::Real=0)
"Partition a graph into two disjoint sets such that the sum of the edge weights
which cross the partition is as large as possible (known to be NP-hard)."
"A cut of a graph can be produced by assigning either 1 or -1 to each vertex. The Goemans-Williamson
algorithm relaxes this binary condition to allow for vector assignments drawn from the (n-1)-sphere
(choosing an n-1 dimensional space will ensure seperability). This relaxation can then be written as
a semidefinite program. Once the optimal vector assignments are found, origin centered hyperplanes are generated and
their corresponding cuts evaluated. After 'iter' trials, or when the desired tolerance is reached,
the hyperplane with the highest corresponding binary cut is used to partition the vertices."
"W: Adjacency matrix."
"tol: Maximum acceptable distance between a cut and the MAXCUT upper bound."
"iter: Maximum number of hyperplane iterations before a cut is chosen."
LinearAlgebra.checksquare(W)
issymmetric(W) || throw(ArgumentError("Adjacency matrix must be symmetric."))
all(W .>= 0) || throw(ArgumentError("Adjacency matrix must be nonnegative."))
all(iszero.(diag(W))) || throw(ArgumentError("Diagonal of adjacency matrix must be zero (no self loops)."))
(tol >= 0) || throw(ArgumentError("The tolerance must be nonnegative."))
(iter > 0) || throw(ArgumentError("The number of iterations must be a positive integer."))
"This is the standard SDP Relaxation of the MAXCUT problem, a reference can be found at,
http://www.sfu.ca/~mdevos/notes/semidef/GW.pdf"
k = size(W, 1)
S = Semidefinite(k)
expr = dot(W, S)
constr = [S[i,i] == 1.0 for i in 1:k]
problem = minimize(expr, constr...)
solve!(problem, SCS.Optimizer; silent_solver=true)
## ensure symmetric positive-definite
A = 0.5 * (S.value + S.value')
A += (max(0, -eigmin(A)) + 1e-14) * Matrix(I, size(A, 1), size(A, 1))
X = Matrix(cholesky(A))
## a non-trivial upper bound on MAXCUT
upperbound = (sum(W) - dot(W, S.value)) / 4
## random origin-centered hyperplanes, generated to produce partitions of the graph
max_cut = -1
max_partition = nothing
for i in 1:iter
gweval = X' * randn(k)
partition = (findall(gweval .>= 0), findall(gweval .< 0))
cut = sum(W[partition...])
if cut > max_cut
max_partition = partition
max_cut = cut
end
(upperbound - max_cut < tol) && break
end
return max_cut, max_partition
end
export maxcut
end
#include "../main.h"
#include <stdlib.h>
#include <assert.h>
#include <cstring>
using std::pair;
const int VERSION = 5;
const int VERSION = 1;
const bool ALLOW_MISTAKES = 1;
struct Nd {
......@@ -13,7 +14,7 @@ struct Nd {
};
void calc_maxcut(Nd * nd, int nd_size)
void calc_maxcut_naive(Nd * nd, int nd_size)
{
fo(i, nd_size) nd[i].color = false;
bool change = true;
......@@ -39,11 +40,65 @@ void calc_maxcut(Nd * nd, int nd_size)
}
}
void calc_maxcut_gw(Nd * nd, int nd_size)
{
auto [p_stdin, p_stdout] = lib::dualpopen("julia");
{
vector<vector<int>> matrix(nd_size, vector<int>(nd_size, 0));
fo(i, nd_size)
for(auto [target, weight] : nd[i].edges)
matrix[i][target-nd] += weight;
fprintf(p_stdin,
"include(\"MaxCut.jl/src/MaxCut.jl\");\
using Random;\
Random.seed!(3);\
MaxCut.maxcut([");
fo(i,nd_size)
{
fo(j, nd_size) fprintf(p_stdin, " %d", matrix[i][j]);
fprintf(p_stdin, ";");
//fo(j, nd_size) fprintf(stderr, " %d", matrix[i][j]);
//fprintf(stderr, ";\n");
}
fprintf(p_stdin, "])[2][1]\n");
}
fflush(p_stdin);
char * line;
int elements;
assert(fscanf(p_stdout, "%d%m[^\n]", &elements, &line)==2);
assert(!strcmp(line, "-element Vector{Int64}:"));
delete [] line;
fo(i, nd_size) nd[i].color = false;
fo(i, elements)
{
int val;
assert(fscanf(p_stdout, "%d", &val)==1);
val--;
//fprintf(stderr, "-> %d\n", val);
nd[val].color = true;
}
fclose(p_stdin);
fclose(p_stdout);
}
void calc_maxcut(Nd * nd, int nd_size, const char * algorithm)
{
if(!strcmp("naive", algorithm)) calc_maxcut_naive(nd, nd_size);
else
if(!strcmp("gw", algorithm)) calc_maxcut_gw(nd, nd_size);
else assert(false);
}
vector<bool> solve(int n, vector<int> input,
int argc, char const* const* argv)
{
int cars_per_vertex = argc>=1?atoi(argv[0]):3;
assert(argc==2);
int cars_per_vertex = atoi(argv[0]);
assert(cars_per_vertex>0);
const char * maxcut_algorithm = argv[1];
vector<int> other = lib::gen_arr_other_car_of_type(n, input);
vector<bool> out(2*n);
int nd_size = (2*n -1)/cars_per_vertex + 1;
......@@ -57,7 +112,7 @@ vector<bool> solve(int n, vector<int> input,
nd[a].edges.push_back({nd + b, 1});
}
calc_maxcut(nd, nd_size);
calc_maxcut(nd, nd_size, maxcut_algorithm);
fo(i, 2*n)
out[i] = nd[i/cars_per_vertex].color;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment