% completely mtrix factorization

clear all;
clc

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

fid_1 = fopen('NPG_run_lambda.txt','w');
fid_2 = fopen('AP_run_lambda.txt','w');

repeats = 100; % repeats > 1

maxiter = 5000;
tol_1 = 1e-16;
tol_2 = 1e-15;
lambda_list =[0 0.2 0.4 0.6 0.8 0.9 0.95 0.96 0.97 0.98 0.99];


%% Begin loop for lambda
for i = 1: length(lambda_list)
    lambda = lambda_list(i);
    successtnpg = 0; successtap = 0;
    fmax_1 = 0;  fmin_1 = inf;
    fmax_2 = 0;  fmin_2 = inf;
    r = 12;
    fopen(fid_1);
    fopen(fid_2);
    progress_1 = []; progress_2 = [];
    
    for repeat = 1: repeats
        % preprocessing
        tstart_pre = tic;
        A = generate(lambda);
        [B_0,p] = chol(A,'lower');
        if p>0
            [U,S] = eig(A);
            B_0 = U*diag(sqrt(diag(S)))*U';
        end
        B = expand(B_0,r);
        
        [U_1,V_1] = svd(randn(r));
        Q_0 = U_1*V_1';
        time0 = toc(tstart_pre);
        fprintf('******** Done preprocessing ******\n')
        
        fprintf('******** NPG for CPMF ********\n')
        
        clear opts
        
        opts.maxiter = maxiter;
        opts.freq = 100;
        opts.disp = 0;
        opts.M = 5;
        opts.tol = tol_1;
        tstart = tic;
        [Q_1, successnpg, fnew, error_1, iter_1, ~] = NPG(B, Q_0, opts);
        time_1 = toc(tstart);
        
        successtnpg = successnpg + successtnpg;
        
        if successnpg == 1
            iter_s_1 = iter_1; iter_f_1 = 0;
            time_s_1 = time_1; time_f_1 = 0;
        else
            iter_f_1 = iter_1; iter_s_1 = 0;
            time_f_1 = time_1; time_s_1 = 0;
        end
        
        fmax_1 = max(fmax_1,fnew);
        fmin_1 = min(fmin_1,fnew);
        
        progress_1 = [progress_1; iter_s_1  iter_f_1 time_s_1 time_f_1];
        
        fprintf('********AP for CPMF ********\n')
        
        clear opts
        
        opts.maxiter = maxiter;
        opts.freq = 100;
        opts.disp = 0;
        opts.tol = tol_2;
        tstart = tic;
        [Q_2,successap, fval, error_2, iter_2, ~] = AP(B, Q_0, opts);
        time_2 = toc(tstart);
        
        successtap = successap + successtap;
        
        if successap == 1
            iter_s_2 = iter_2; iter_f_2 = 0;
            time_s_2 = time_2; time_f_2 = 0;
        else
             iter_f_2 = iter_2; iter_s_2 = 0;
             time_f_2 = time_2; time_s_2 = 0;
        end
        
        fmax_2 = max(fmax_2,fval);
        fmin_2 = min(fmin_2,fval);
        
        progress_2 = [progress_2; iter_s_2 iter_f_2 time_s_2 time_f_2];
        
     end
   
    fprintf(fid_1, '%3.2f  & %5d  & %5.0f & %1.0e & %1.0e & %5.0f & %5.0f & %5.4f & %5.4f\n',lambda, r, successtnpg/repeats*100, fmax_1, fmin_1, sum(progress_1(:,1))/ successtnpg,sum(progress_1(:,2))/ (repeats-successtnpg),sum(progress_1(:,3))/ successtnpg,sum(progress_1(:,4))/ (repeats-successtnpg));
    fprintf(fid_2, '%3.2f  & %5d  & %5.0f & %1.0e & %1.0e & %5.0f & %5.0f & %5.4f & %5.4f\n',lambda, r,  successtap/repeats*100, fmax_2, fmin_2, sum(progress_2(:,1))/ successtap,sum(progress_2(:,2))/ (repeats-successtap),sum(progress_2(:,3))/ successtap,sum(progress_2(:,4))/ (repeats-successtap));
end
fclose(fid_1);
fclose(fid_2);