% TEST CLUSTERING:  Main test file for solving the clustering problem:
%                       minimize        < C,X >
%                       subject to:     Trace(X) = k
%                                       X positive-semi-definite
%                                       sum(X,2) = ones
%                                       X elementwise nonnegative
%
% SDP relaxation of clustering problem is solved with the preprocessed
% data of MNIST digits[1]. 
% TODO(mfsahin): Add also synthetic data option and other real data options as
%                shown in [2]
%
% [1]: D. G. Mixon, S. Villar, and R. Ward. Clustering sub- gaussian
%      mixtures by semidefinite programming. In- formation and Inference:
%      A Journal of the IMA, 6 (4):389?415, 2017
% [2]: Mariano Tepper, Anirvan M. Sengupta, and Dmitri B. Chklovskii,
%      The surprising secret identity of the semidefinite relaxation of       
%      k-means: manifold learning, CoRR abs/1706.06028 (2017)
% 
% Author: Mehmet Fatih Sahin
% Date  : 02.11.2018

clearvars
close all

vec = @(x)reshape( x, numel( x ), 1 );
addpath(genpath('./helpers/'));
addpath(genpath('./functions/'));
%%
clustering_load_data();
% NOTE(mfsahin): we might try normalizing C
normC = normest(C);
%% function 
% TODO(mfsahin): remove Aef2(u) - b2 condition since it is redundant.
Aef = @(U) U*sum(U, 1)';
Aef2 = @(U) sum(U,1)*U';

r = 20;

f_grad_u = @(U, beta_, yc) C*U + 1/beta_*( repmat(U*sum(U, 1)' - b, [1, r]).*repmat(sum(U,1), [N, 1]) + repmat((sum(U,1)*U' -b2)*U, [N,1]))...
        + repmat(yc(1:size(b,1)), [1, r]).*repmat(sum(U,1), [N, 1]) + repmat(yc(size(b,1)+1:2*size(b,1))'*U, [N,1]);
f_u = @(u) sum(vec(u .* (C*u)));
proj = @(u) clustering_proj(u, k);
get_feas = @(u) [Aef(u) - b; (Aef2(u)-b2)'];
gamma_rule = @(beta_, yc) .5*beta_/(2*N + beta_*normC + sqrt(2*N)*norm(beta_*yc-1));

%% Define algorithmic parameters and rules
beta0 = 1e-2; % initial value for the penalty parameter
sigma0 = 1e2*beta0; % initial value for dual step size
%%

update_params = @(feas_cur, feas_old, k_t, iter, beta_old, sigma_old, sigma_vec, beta_vec) update_params(feas_cur, feas_old, k_t, iter, beta_old, sigma_old, sigma_vec, beta_vec, sigma0, beta0);

params.averaging = 0;
params.iters = 1e5;% Number of iterations
params.verbose = 1e2;

params.Uinit = randn(N, r);     % In case a random initialization is used
params.ycinit = zeros(2*size(b,1),1);
params.savehist = 1;
params.startfeasible = 0; % this option is for maxcut
%% Run the algorithm
[Udl, output] = linearized_augmented_lagrangian(f_u, f_grad_u, proj, update_params, ...
    gamma_rule, get_feas, params);
%% Plot results
opt_val = 77.190856515019890;
plot_results();