diff --git a/MaxCut.jl/LICENSE.md b/MaxCut.jl/LICENSE.md
new file mode 100644
index 0000000000000000000000000000000000000000..0fab6ae24acb4f128c5576ffb5a78fadaf81f303
--- /dev/null
+++ b/MaxCut.jl/LICENSE.md
@@ -0,0 +1,28 @@
+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
diff --git a/MaxCut.jl/Project.toml b/MaxCut.jl/Project.toml
new file mode 100644
index 0000000000000000000000000000000000000000..c08a1789f34eb3b14c110fb76c9f1a3091812af1
--- /dev/null
+++ b/MaxCut.jl/Project.toml
@@ -0,0 +1,12 @@
+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"
diff --git a/MaxCut.jl/README.md b/MaxCut.jl/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..563bfc3a281b1eecf9e167879b357f2266e3833d
--- /dev/null
+++ b/MaxCut.jl/README.md
@@ -0,0 +1,46 @@
+**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
diff --git a/MaxCut.jl/src/MaxCut.jl b/MaxCut.jl/src/MaxCut.jl
new file mode 100644
index 0000000000000000000000000000000000000000..4cd806f1a6430ac3d7de1320855fa23a665f6681
--- /dev/null
+++ b/MaxCut.jl/src/MaxCut.jl
@@ -0,0 +1,71 @@
+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
diff --git a/algo/maxcut.cpp b/algo/maxcut.cpp
index 1983bb3bb64e9babbf7c22b68253e9c116d849a3..09e3854ebab089b2533ab96a787495233c437938 100644
--- a/algo/maxcut.cpp
+++ b/algo/maxcut.cpp
@@ -1,10 +1,11 @@
#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;