clear;clc
close all;

n = 300; % n dimension
r = 5; % r number of column
mu = .8; % mu sparse parameter

epsilon_set = [1e-1;1e-2;1e-3;1e-4;1e-5]; % mu

seed = round(rand() * 10000000);
fprintf('seed:%d\n', seed);
rng(seed);

randnum = 100; 
tab = zeros(3, length(epsilon_set), randnum);
maxiter = 5000;
outputgap = 10;
numLSfails = 0;

for j = 1 : randnum
    seed = round(rand() * 10000000);
    fprintf('seed:%d\n', seed);
    rng(seed);
    % generate the random data matrix A
    m = 50;
    A = randn(m,n);
    A = A - repmat(mean(A,1),m,1);
    A = normc(A);
    fprintf('n:%d, r:%d, mu:%f, seed:%d\n', n, r, mu, seed);
    % random initialization
    [phi_init, ~] = svd(randn(n,r),0); 
    x0 = phi_init;
    option.n = n; option.r = r; option.mu = mu;
    option.tol = 1e-8*n*r;
    option.maxiter = maxiter;
    option.x0 = x0;
    option.stop = 1e-10;
    option.outputgap = outputgap;
    for i = 1 : length(epsilon_set)
        option.epsilon = epsilon_set(i);

        %% Drive_RPN-H
        [x_rpnH, fs_rpnH, nv_rpnH, iter_rpnH, sparsity_rpnH, time_rpnH, iter_time_rpnH, status] = driver_rpn_h(A, option);
        Fs_rpnH = fs_rpnH(end);
        Ns_rpnH = nv_rpnH(end);
        if Ns_rpnH <= option.stop
            tab(1, i, j) = 1;
        end
    end

    [x_rpnncg, fs_rpnncg, nv_rpnncg, iter_rpnncg, sparsity_rpnncg, time_rpnncg, iter_time_rpnncg] = driver_rpnncg(A, option);
    Ns_rpnncg = nv_rpnncg(end);
    if Ns_rpnncg <= option.stop
        for i = 1 : length(epsilon_set)
            tab(2, i, j) = 1;
        end
    elseif iter_rpnncg < maxiter
        numLSfails = numLSfails + 1;
    end

    [x_rpncg, fs_rpncg, nv_rpncg, iter_rpncg, sparsity_rpncg, time_rpncg, iter_time_rpncg] = driver_rpncg(A, option);
    Ns_rpncg = nv_rpncg(end);
    if Ns_rpncg <= option.stop
        for i = 1 : length(epsilon_set)
            tab(3, i, j) = 1;
        end
    end
end

taball = sum(tab, 3)

numLSfails

% 1: RPN-H, 2:RPN-NCG, 3:RPN-CG
fout = fopen('RPNH_vs_RPNNCG_vs_RPNCG.txt','w');
for i = 1 : size(tab, 1)
    for j = 1 : (size(tab, 2) + 1)
        if(j == 1)
            fprintf(fout, '%d & ', i);
        elseif(j>=2 && j <= size(tab, 2))
            fprintf(fout, '%d\\%% & ', taball(i, j - 1));
        else
            fprintf(fout, '%d\\%% \\\\ ', taball(i, j - 1));
        end
    end
    fprintf(fout,'\n');
end
fclose(fout);


function str = outputfloat(x)
    if(x <= 0)
        sn = '-';
        x = abs(x);
    else
        sn = '';
    end
    p = log(x)/log(10);
    p = - ceil(-p);
    x = round(x * 10^(-p) * 100);
    x = x / 100;
    strx = sprintf('%3.2f', x);
    if(p ~= 0)
        str = [sn strx '_{' num2str(p) '}'];
    else
        str = [sn strx];
    end
end




