% test code
delta = 1e-3;

p = 0.5;
tol = 1e-6;
lambda0 = 1;
mu0 = 1; 
eps0 = 1;
rho = 2; 
theta = 1/rho; 
M = 5;
maxiter = inf;
verbose = 0;
repeat = 10;

% for spgL1
opts = spgSetParms('verbosity',0);

rand('seed',3);
randn('seed',3);
fprintf(' version 9\n')

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

for index = 12:2:30
   % signal length
   N = 512*index;
   % number of spikes to put down
   T = 20*index;
   % number of observations to make
   K = 120*index;
   
   % reinitialize variables for logging
   fvalP = 0; fvallp2 = 0;
   nonzP = 0; nonzSPGL1 = 0; nonzlp2 = 0;
   errP = 0; errSPGL1 = 0; errlp2 = 0;
   timeP = 0; timeSPGL1 = 0; timelp2 = 0;

   lambdaP = 0;
   
   % generate multiple random instances
   for i = 1:repeat
     
     % random Gaussian signal
     x_orig = zeros(N,1);
     q = randperm(N);
     x_orig(q(1:T)) = randn(T,1);
     % measurement matrix
     A = randn(K,N);
     A = orth(A')'; % orthogonalize
     % noisy observations
     e = delta*randn(K,1);
     b = A*x_orig + e;
     sigma = norm(e);
     if sigma >= norm(b)
       error(' failure: contains zero\n')
     end
%      xtrue = x;

     x0 = ones(N,1);
     feas_pt = A\b;
   
     fprintf(' **** Finish generation: N = %g ***\n', N)
   
     [x,fval,residual,iter,nf,lambda,time] = Lp(A,b,sigma,p,tol,lambda0,mu0,eps0,rho,theta,M,maxiter,feas_pt,x0,0);
     fvalP = fvalP + norm(x,p)^p;
     nonzP = nonzP + nnz(x);
     errP = errP + norm(x - x_orig);
     timeP = timeP + time;
     lambdaP = lambdaP + lambda;
     
     [x3,fval3,residual3,iter3,nf3,lambda3,time3] = Lp2(A,b,sigma,p,tol,lambda0,eps0,rho,theta,M,maxiter,x0,verbose);
     
     fvallp2 = fvallp2 + norm(x3,p)^p;
     nonzlp2 = nonzlp2 + nnz(x3);
     errlp2 = errlp2 + norm(x3 - x_orig);
     timelp2 = timelp2 + time3;
   
     tstart = tic;
     [x2,~,~,info] = spg_bpdn(A, b, sigma, opts);
     time2 = toc(tstart);
     residual2 = max(norm(A*x2-b)^2 - sigma^2, 0);
     
     nonzSPGL1 = nonzSPGL1 + nnz(x2);
     errSPGL1 = errSPGL1 + norm(x2 - x_orig);
     timeSPGL1 = timeSPGL1 + time2;
     
     
%      fprintf(fid, '%5d & %10.0f & %2.1e & %10.2f \n', info.iter, nnz(x2), residual2, time2);
     fprintf('\n\n --------------- The computational result for SPGL1 ----------------\n\n')
     fprintf(' Iter = %3i err = %7.3f  nnz=%d   time = %.2f\n', ....
       info.iter, norm(x2 - x_orig), nnz(x2),time2)
     fprintf('\n ----------------------------------------------------------\n\n')
     
   end
   
   % average and save outputs
   
   fprintf(fid,' %13d & %13d & %13d &', K, N, T);
   fprintf(fid, '%10.0f & %2.1e & %10.2f &', nonzSPGL1/repeat, errSPGL1/repeat, timeSPGL1/repeat);
   fprintf(fid, '%3.2e & %10.0f & %2.1e & %10.2f &', fvallp2/repeat, nonzlp2/repeat, errlp2/repeat, timelp2/repeat);
   fprintf(fid, '%3.2e & %10.0f & %2.1e & %10.2f & %3.2e\n', fvalP/repeat, nonzP/repeat, errP/repeat, timeP/repeat, lambdaP/repeat);
end
fclose(fid);