% completely mtrix factorization for random generated instace
% using random generated initial points

clear all;
clc

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

a = clock;
fname1 = ['NPG_run_random-instance-randinits' '-' date '-' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fname2 = ['AP_run_random-instance-randinits' '-' date '-' int2str(a(4)) '-' int2str(a(5)) '.txt'];

fid_1 = fopen(fname1,'w');
fid_2 = fopen(fname2,'w');

repeats = 50; % repeats > 1

maxiter = 5000;
tol_1 = 1e-16;
tol_2 = 1e-15;
n_list = [10 20 30 40 100 200 300 400];

for i = 1: length(n_list)
    n = n_list(i);
    r = 1.5*n_list(i);
    if n <= 50 
        init_no = 100;
    else
        init_no = 10;
    end
    successtnpg = 0; successtap = 0;
    fmax_1 = 0;  fmin_1 = inf;
    fmax_2 = 0;  fmin_2 = inf;
    
    fopen(fid_1);
    fopen(fid_2);
    progress_1 = []; progress_2= [];
    
    for repeat = 1: repeats
        % preprocessing
        tstart_pre = tic;
        A_0 = abs(randn(n,2*n));
        A = A_0*A_0';
        [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);
        Q_0 = randn(r);
        [U_Q,~,V_Q] = svd(Q_0);
        Q_init = U_Q*V_Q';
        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;
        iter_1 = 0;
        Q_0 = Q_init;
        tstart = tic;
        for ips_1 = 1:init_no
            [Q_1, successnpg, fnew, error_1, iter_1_i, ~] = NPG(B, Q_0, opts);
            iter_1 = iter_1 + iter_1_i;
            if successnpg == 1 || ips_1 == init_no
                break;
            else
                Q = randn(r);
                [U_i, ~, V_i] = svd(Q);
                Q_0 = U_i*V_i';
            end
        end
        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;
            ips_s_1 = ips_1; ips_f_1  = 0;
        else
            iter_f_1 = iter_1; iter_s_1 = 0;
            time_f_1 = time_1; time_s_1 = 0;
            ips_f_1 = ips_1; ips_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 ips_s_1];
        
        fprintf('********AP for CPMF ********\n')
        
        clear opts
        
        opts.maxiter = maxiter;
        opts.freq = 100;
        opts.disp = 0;
        opts.tol = tol_2;
        iter_2 = 0;
        Q_0 = Q_init;
        tstart = tic;
        for ips_2  = 1:init_no
            [Q_2,successap, fval, error_2, iter_2_i, ~] = AP(B, Q_0, opts);
            iter_2 = iter_2 + iter_2_i;
            if successap == 1 || ips_2 == init_no
                break;
            else 
                Q = randn(r);
                [U_i, ~, V_i] = svd(Q);
                Q_0 = U_i*V_i';
            end
        end 
        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;
            ips_s_2 = ips_2; ips_f_2 = 0;
        else
             iter_f_2 = iter_2; iter_s_2 = 0;
             time_f_2 = time_2; time_s_2 = 0;
             ips_f_2 = ips_2; ips_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  ips_s_2];
        
     end
   
    fprintf(fid_1, '%5.0f  & %5d  & %5.0f & %1.0e & %1.0e & %5.0f & %5.0f & %5.4f & %5.4f & %5.1f \n',n, 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), sum(progress_1(:,5))/successtnpg);
    fprintf(fid_2, '%5.0f  & %5d  & %5.0f & %1.0e & %1.0e & %5.0f & %5.0f & %5.4f & %5.4f & %5.1f \n', n, 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), sum(progress_2(:,5))/successtap);
end
fprintf(fid_1, '\n\n tol = %1.0e', tol_1);
fprintf(fid_2, '\n\n tol = %1.0e', tol_2);
fclose(fid_1);
fclose(fid_2);

