% Run code for simultaneous sparse and low rank

clear all
clc

randn('seed',1);
rand('seed',1);

fid = fopen('lrsp_final.txt','w');

%% Parameters for data
m_list = [1000 2000 3000];
n = 500; k = 10;
p = 0.1; 

lambda_init = 0.1;
repeats = 10;

fopen(fid);
for sigma = [0.005 0.01 0.02]
    for index = 1:3
        iter_rank = 0; time_rank = 0; err_rank = 0; Fval_rank = 0;
        iter_sparse = 0; time_sparse = 0; err_sparse = 0; Fval_sparse = 0;
        
        m = m_list(index);s = m*n*p; num_s = (1-p)*m;
        fprintf('m = %g, n = %g, k = %g, p = %g, sigma = %g, Start Now!\n',m, n, k, p, sigma);
        
        for repeat = 1: repeats
            %% Generate sparse and low rank matrix with noise
            I = randperm(m);
            M1 = randn(m,k); M1(I(1:num_s),:) = 0;
            M2 = randn(k,n);
            M = M1*M2; Mnoise = M + sigma*randn(m,n);
            
            t_rank = tic;
            fprintf('****** SDCAM for low rank constraint ******\n')
            clear opts
            opts.X0 = zeros(m,n); opts.tol = 1e-5; opts.freq = 100000; opts.maxiter = 10000;
            opts.disp = 1; opts.nonM = 5; lambda = lambda_init;
            
            Fval_X0 = 0.5*norm(opts.X0 - Mnoise,'fro')^2;
            Xfeas = opts.X0; Fval_feas = Fval_X0;
            iter_1 = 0;
            while 1 == 1
                if Fval_feas < Fval_X0
                    fprintf('swap initial point, lambda = %g\n', lambda)
                    opts.X0 = Xfeas;
                end
                [X,err,fval_ls,iter,time] = LRS_rank(Mnoise,k,s,lambda,opts);
                iter_1 = iter_1 + iter;
                fprintf(' One inner loop done, lambda = %g tol = %g\n', lambda, opts.tol)
                lambda = lambda/10;
                opts.X0 = X;
                opts.tol = max(opts.tol/1.5,1e-6);
                if err < 1e-6*norm(X,'fro')
                    break
                end
                Fval_X0 = fval_ls + err^2/(2*lambda);
            end
            iter_rank = iter_rank + iter_1;
            time_rank = time_rank + toc(t_rank);
            err_rank = err_rank + err;
            Fval_rank = Fval_rank + fval_ls;
            
            
            
            t_sparse = tic;
            fprintf('****** SDCAM for sparse constraint ******\n')
            clear opts
            opts.X0 = zeros(m,n); opts.tol = 1e-5; opts.freq = 100000; opts.maxiter = 10000;
            opts.disp = 1; opts.nonM = 5; lambda = lambda_init;
            
            Fval_X0 = 0.5*norm(opts.X0 - Mnoise,'fro')^2;
            Xfeas = opts.X0; Fval_feas = Fval_X0;
            iter_2 = 0;
            while 1 == 1
                if Fval_feas < Fval_X0
                    fprintf('swap initial point, lambda = %g\n', lambda)
                    opts.X0 = Xfeas;
                end
                [X_1,err,fval_ls,iter,time] = LRS_sparse(Mnoise,k,s,lambda,opts);
                iter_2 = iter_2 + iter;
                fprintf(' One inner loop done, lambda = %g tol = %g\n', lambda, opts.tol)
                lambda = lambda/10;
                opts.X0 = X_1;
                opts.tol = max(opts.tol/1.5, 1e-6);
                if err < 1e-6*norm(X_1,'fro')
                    break
                end
                Fval_X0 = fval_ls + err^2/(2*lambda);
            end
            iter_sparse = iter_sparse + iter_2;
            time_sparse = time_sparse + toc(t_sparse);
            err_sparse = err_sparse + err;
            Fval_sparse = Fval_sparse + fval_ls;
        end
        
        fprintf(fid, ' %4.3f & %g & %5.0f & %5.0f & %5.1f & %5.1f & %5.4e & %5.4e\n', sigma, m, iter_rank/repeats, iter_sparse/repeats,...
            time_rank/repeats, time_sparse/repeats, err_rank/repeats, err_sparse/repeats);
    end
end

fclose(fid);


