% Test script for FBE


clear all
clc

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

fprintf('new\n')

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

repeats = 10;

path('./FBE_alg_new',path);

mu = 1e-3;
sigma = 1e-2; % noise level; fixed
lambda = 1;
M = 10;


% parameters for different FBE
beta = 0.5;
beta_1 = 0.7;
beta_2 = 0.9;

totalcount = 0;

fopen(fid);

for index = 1:10
    
  iter_FBE = 0;
  time_FBE = 0;
  fval_FBE = 0;
  
  iter_FBE1 = 0;
  time_FBE1 = 0;
  fval_FBE1 = 0;
  
  iter_FBE2 = 0;
  time_FBE2 = 0;
  fval_FBE2 = 0;
  
 
  
  time_lambda = 0;
  
  % signal length
  n = 2560*index;
  % number of spikes to put down
  k = 160*index; % originally 80*index
  % number of observations to make
  m = 720*index;
  
  fprintf(' Size: m, n, k = %g, %g, %g\n', m, n, k)
  
  for repeat = 1:repeats
    
    I = randperm(n);
    J = I(1:k);
    x = zeros(n,1);
    x(J) = randn(k,1);
    
    A = randn(m,n);
    for j = 1:n
        A(:,j) = A(:,j)/norm(A(:,j));
    end
    b = A*x + sigma*randn(m,1);
    
    if m > 2000
        clear opts
        opts.issym = 1;
        tstart = tic;
        nmA = eigs(A*A',1,'LM',opts);
        time_lambda = time_lambda + toc(tstart);
    else
        tstart = tic;
        nmA = norm(A*A');
        time_lambda = time_lambda + toc(tstart);
    end
    fprintf('\n Lipschitz constant L = %g\n', nmA)
    
    fprintf(' ********** L1-L2 FB_envelope with L-BFGS **********\n')
    clear opts
    opts.x0 = zeros(2*n,1);
    opts.freq = inf;
    opts.maxiter = inf;
    opts.disp = 0;
    opts.L = (nmA + sqrt(nmA^2 + 4*mu^2*lambda^2))/2*0.95/beta;
    opts.M = M;
    opts.tol = 1e-6;
    
    [x,iter,time,count] = FB_env_ls_l1_l2_lbfgs(A,b,mu,lambda,opts);
    iter_FBE = iter_FBE + iter;
    time_FBE = time_FBE + time;
    fval = 0.5*norm(A*x(1:n) - b)^2 + mu*norm(x(1:n),1) - mu*lambda*norm(x(1:n));
    fval_FBE = fval_FBE + fval;
    
    totalcount = totalcount + count;
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter, time, fval);    
    
    opts.L = (nmA + sqrt(nmA^2 + 4*mu^2*lambda^2))/2*0.95/beta_1;
    [x1,iter1,time1,count1] = FB_env_ls_l1_l2_lbfgs(A,b,mu,lambda,opts);
    iter_FBE1 = iter_FBE1 + iter1;
    time_FBE1 = time_FBE1 + time1;
    fval1 = 0.5*norm(A*x1(1:n) - b)^2 + mu*norm(x1(1:n),1) - mu*lambda*norm(x1(1:n));
    fval_FBE1 = fval_FBE1 + fval1;
    
    totalcount = totalcount + count1;
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter1, time1, fval1);   
    
    opts.L = (nmA + sqrt(nmA^2 + 4*mu^2*lambda^2))/2*0.95/beta_2;
    [x2,iter2,time2,count2] = FB_env_ls_l1_l2_lbfgs(A,b,mu,lambda,opts);
    iter_FBE2 = iter_FBE2 + iter2;
    time_FBE2 = time_FBE2 + time2;
    fval2 = 0.5*norm(A*x2(1:n) - b)^2 + mu*norm(x2(1:n),1) - mu*lambda*norm(x2(1:n));
    fval_FBE2 = fval_FBE2 + fval2;
    
    totalcount = totalcount + count2;
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter2, time2, fval2); 
    
  end
  
  fprintf(fid,' %7.0f & %7.0f & %7.0f & %5.0f & %5.0f & %5.0f & %5.1f & %5.1f & %5.1f & %6.5e & %6.5e & %6.5e\n', m, n, k,...
      iter_FBE/repeats, iter_FBE1/repeat, iter_FBE2/repeat, time_FBE/repeats, time_FBE1/repeats, time_FBE2/repeats, fval_FBE/repeats, fval_FBE1/repeats, fval_FBE2/repeats);  

end

fclose(fid);

save totalcount totalcount