% Example script of adding concave quadratic cuts for max-cut instance
% Copyright (c) 2015. Jaehyun Park and Stephen Boyd

function maxcut_example()
    [n, W] = read_graph('01.mc');
    n = 15; W = W(1:n, 1:n);

    tic;
    % Goemans-Williamson SDP relaxation
    cvx_solver mosek % Needs CVX professional license
    cvx_begin quiet
        variables X(n, n) x(n)
        expression Z
        Z = [X, x; x', 1];
        maximize (sum(W*x) - trace(W*X))
        subject to
            Z == semidefinite(n+1);
            diag(X) == x;
    cvx_end
    ub1 = cvx_optval;
    X = (X + X')/2; % Force symmetry
    
    % Add concave quadratic cuts
    A = []; B = [];
    for iters = 1:10000
        a = zeros(n, 1);
        i = randi(n); j = randi(n); k = randi(n);
        a(i) = a(i) + sign(randn());
        a(j) = a(j) + sign(randn());
        a(k) = a(k) + sign(randn());
        b = floor(a'*x);
        if a'*X*a - (2*b+1)*a'*x + b*(b+1) < 0 % found cut
           A = [A, a]; B = [B; b];
        end
    end

    % Use CVX to solve SDP
    cvx_begin quiet
        variables X(n, n) x(n)
        expression Z
        Z = [X, x; x', 1];
        maximize (sum(W*x) - trace(W*X))
        subject to
            Z == semidefinite(n+1);
            diag(X) == x;
            % Constraints corresponding to concave quadratic cuts
            diag(A'*X*A)-(2*B+1).*(A'*x)+B.*(B+1) >= 0;
    cvx_end
    ub2 = cvx_optval;
    
    % Randomized algorithm
    mu = x; Sigma = X - x*x';
    % Using eigenvalue decomposition because chol() complains about
    % near-zero negative eigenvalues
    [V, D] = eig(Sigma);
    A = V*sqrt(max(D, 0));
    
    K = 3*n;
    lb = -inf; xhat = zeros(n, 1);
    for k = 1:K
        x = round(mulrandn_cached(mu, A));
        x = max(min(x, 1), 0);
        [x, val] = ONEOPT(x, W);
        if lb < val
            lb = val; xhat = x;
        end
    end
    elapsed = toc;
    
    fprintf('Elapsed time: %.5f seconds\n', elapsed);
    fprintf('Goemans-Williamson bound   : %.5f\n', ub1);
    fprintf('SDP lower bound (with cuts): %.5f\n', ub2);
    fprintf('Lower bound                : %.5f\n', lb);
end

% read the graph data file and returns the edge weight matrix
function [n, A] = read_graph(filename)
    fid = fopen(filename);
    if fid == -1
        n = 0; A = [];
        return;
    end
    n = fscanf(fid, '%d', 1);
    m = fscanf(fid, '%d', 1);
    edges = (fscanf(fid, '%d %d %d\n', [3, Inf]))';
    fclose(fid);
    A = zeros(n, n);
    for i = 1:m
        x = edges(i, 1); y = edges(i, 2); w = edges(i, 3);
        A(x, y) = A(x, y) + w;
        A(y, x) = A(y, x) + w;
    end
end


% Subroutine for multivariate normal distribution
%
% Samples from a multivariate normal distribution with mean mu
% and square root matrix of covariance, A (covariance = A*A').
%
% Reference:
% http://en.wikipedia.org/wiki/Multivariate_normal_distribution

function x = mulrandn_cached(mu, A)
    n = size(mu, 1);
    z = randn(n, 1);
    x = mu + A*z;
end


% 1-opt greedy descent subroutine
%
% Starting from x, greedily chooses and optimizes over a single coordinate 
% to maximize the objective sum(W x) - x^T W x the most. Stops when no
% single coordinate change can reduce the objective. Returns the 1-opt
% point and the function value.

function [x, val] = ONEOPT(x, W)
    val = sum(W*x) - x'*W*x;
    v = diag(W);
    while true
        dels = (W*(1-2*x)).*(1-2*x) - v;
        [del, k] = max(dels);
        if del <= 0
            break;
        end
        x(k) = 1 - x(k);
        val = val + del;
    end
end
