clear;clc
close all;

seed = round(rand()*10000);
% seed = 4170
rng(seed);

r = 4; % r number of column
mu = 3; % mu sparse parameter



dat = importdata('dataBCT.txt');
D = (dat.data)';
[m, n] = size(D);
A = D - ones(size(D,1),1)*mean(D, 1);
tmp = sqrt(sum(A .* A));
A = A ./ repmat(tmp, m, 1);

seed = round(rand() * 10000000);
seed = 1;
rng(seed);


% random initialization
% x0 = randn(n,r);
% x0 = x0/norm(x0,'fro');
[U, S, V] = svd(A', 0);
x0 = U(:, 1 : r);

fid = 1;
maxiteration = 5000;
outputgap = 100;

fprintf('n:%d, r:%d, m:%d, mu:%f\n', n, r, m, mu);

%% Drive_ManPG
option.n = n; option.r = r; option.mu = mu;
option.tol = 1e-8*n*r;
option.maxiter = maxiteration;
option.x0 = x0;
option.stop = 1e-10;
option.outputgap = outputgap;

[x_manpg, fs_manpg, nv_manpg, iter_manpg,...
    sparsity_manpg, time_manpg,iter_time_manpg] = driver_ManPG(A, option);
Fs_manpg = fs_manpg(end);
Ns_manpg = nv_manpg(end);

%% Drive_ManPG_Ada
[x_manpg_ada, fs_manpg_ada, nv_manpg_ada,...
    iter_manpg_ada, sparsity_manpg_ada, time_manpg_ada,iter_time_manpg_ada] = driver_ManPG_ada(A, option);
Fs_manpg_ada = fs_manpg_ada(end);
Ns_manpg_ada = nv_manpg_ada(end);


%% ManPQN
option.type =0;
option.inner_iter = 200;
M = 5;
[X_pqn, F_pqn,F_pqn_list,sp_pqn,t_pqn,...
    maxit_att_pqn,succ_flag_pqn,lins_pqn,in_av_pqn,nv_pqn, iter_time_pqn]= manpqn_orth_sparse(A,option,M,1);
succ_no_pqn = 1;
iter_pqn = sum(maxit_att_pqn)/succ_no_pqn;
time_pqn = sum(t_pqn)/succ_no_pqn;
Fs_pqn   = sum(F_pqn)/succ_no_pqn;
sparsity_pqn   = sum(sp_pqn)/succ_no_pqn;
Ns_pqn   = nv_pqn(end);

%% Drive_RPN-CG
[x_rpncg, fs_rpncg, nv_rpncg, iter_rpncg,...
    sparsity_rpncg, time_rpncg, iter_time_rpncg] = driver_rpncg(A, option);
Fs_rpncg = fs_rpncg(end);
Ns_rpncg = nv_rpncg(end);
% return;

fprintf(fid,' Alg ****        Iter *****  Fval ******* sparsity ***** cpu ******normv**\n');
print_format =  'ManPG     : &  %e   & %1.5e  &   %1.2f    &   %3.2f   & %1.3e\n';
fprintf(fid, print_format, iter_manpg , Fs_manpg, sparsity_manpg, time_manpg, Ns_manpg);
print_format =  'ManPG-Ada : &  %e   & %1.5e  &   %1.2f    &   %3.2f   & %1.3e\n';
fprintf(fid, print_format, iter_manpg_ada , Fs_manpg_ada, sparsity_manpg_ada, time_manpg_ada, Ns_manpg_ada);
print_format =  'ManPQN    : &  %e   & %1.5e  &   %1.2f    &   %3.2f   & %1.3e\n';
fprintf(fid,print_format, iter_pqn , Fs_pqn, sparsity_pqn, time_pqn, Ns_pqn);
print_format =  'RPN-CG    : &  %e   & %1.5e  &   %1.2f    &   %3.2f   & %1.3e \n';
fprintf(fid, print_format, iter_rpncg,  Fs_rpncg, sparsity_rpncg, time_rpncg, Ns_rpncg);

clf;
set(0,'defaultaxesfontsize',15, ...
   'defaultaxeslinewidth',0.7, ...
   'defaultlinelinewidth',.8,'defaultpatchlinewidth',0.8);
set(0,'defaultlinemarkersize',10)

subplot(2,1,1)

semilogy(1:length(nv_manpg), nv_manpg,'r-o','LineWidth',1.3, 'MarkerIndices', 1:300:length(nv_manpg))
hold on
semilogy(1:length(nv_manpg_ada), nv_manpg_ada,'k-s','LineWidth',1.3, 'MarkerIndices', 1:300:length(nv_manpg_ada))
hold on
semilogy(1:length(nv_pqn), nv_pqn,'m->', 'LineWidth',1.3, 'MarkerIndices', 1:300:length(nv_pqn))
hold on
sspa = [1:50:length(nv_rpncg), length(nv_rpncg)-20:length(nv_rpncg)];
semilogy(1:length(nv_rpncg), nv_rpncg,'b-+', 'LineWidth',1.3, 'MarkerIndices', sspa)
legend('ManPG','ManPG-Ada','ManPQN','RPN-CG')
xlabel('Iter');
ylabel('$\|v(x_k)\|$','Interpreter','latex')
title(' DNA data ')
 
subplot(2,1,2)

semilogy(iter_time_manpg, nv_manpg,'r-o','LineWidth',1.3, 'MarkerIndices', 1:300:length(iter_time_manpg))
hold on
semilogy(iter_time_manpg_ada, nv_manpg_ada,'k-s','LineWidth',1.3, 'MarkerIndices', 1:300:length(iter_time_manpg_ada))
hold on
semilogy(iter_time_pqn, nv_pqn,'m->', 'LineWidth',1.3, 'MarkerIndices', 1:300:length(iter_time_pqn))
hold on
sspa = [1:50:length(iter_time_rpncg), length(iter_time_rpncg)-20:length(iter_time_rpncg)];
semilogy(iter_time_rpncg, nv_rpncg,'b-+', 'LineWidth',1.3, 'MarkerIndices', sspa)
legend('ManPG','ManPG-Ada','ManPQN','RPN-CG')
xlabel('CPU time(s)');
ylabel('$\|v(x_k)\|$','Interpreter','latex')
title(' DNA data ')


