% Test script for FBE


clear all
clc

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

fprintf('new\n')

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

repeats = 10;

path('./FBE_alg_new',path);

mu = 5e-4;
sigma = 1e-2; % noise level; fixed
lambda = 1;

M = 10; % changed on June 2, 2016; changed back on June 4, 2016
M2 = 5;
M3 = 5;

fopen(fid);
count10 = 0;

for index = 1:10
    
  iter_FBE = 0;
  time_FBE = 0;
  fval_FBE = 0;
  
  iter_DC = 0;
  time_DC = 0;
  fval_DC = 0;
  
  iter_major = 0;
  time_major = 0;
  fval_major = 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;
    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;
    
    count10 = count10 + count;
    
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter, time, fval);    
    
    fprintf(' ********** L1 - L2: FBS and majorization hybrid with line search**********\n')
    clear opts
    opts.x0 = zeros(n,1);
    opts.tol = 1e-4;
    opts.freq = inf;
    opts.maxiter = inf;
    opts.disp = 0;
    opts.L = nmA;
    opts.M = M2;
    
    [x2,iter2,time2] = FBS_major_ls_l1_l2(A,b,mu,lambda,opts);
    iter_DC = iter_DC + iter2;
    time_DC = time_DC + time2;
    fval2 = 0.5*norm(A*x2 - b)^2 + mu*norm(x2,1) - mu*lambda*norm(x2);
    fval_DC = fval_DC + fval2;
    
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter2, time2, fval2); 
    
    fprintf(' ********** L1 - L2: Majorization with line search**********\n')
    clear opts
    opts.x0 = zeros(n,1);
    opts.tol = 1e-4;
    opts.freq = inf;
    opts.maxiter = inf;
    opts.disp = 0;
    opts.L = nmA;
    opts.M = M3;
    
    [x3,iter3,time3] = major_ls_l1_l2(A,b,mu,lambda,opts);
    iter_major = iter_major + iter3;
    time_major = time_major + time3;
    fval3 = 0.5*norm(A*x3 - b)^2 + mu*norm(x3,1) - mu*lambda*norm(x3);
    fval_major = fval_major + fval3;
    
    fprintf(' iter = %5d  time = %5.1f  fval = %8.7e\n', iter3, time3, fval3); 
    
  end
  
  fprintf(fid,' %7.0f & %7.0f & %7.0f & %5.1f & %5.0f & %5.0f & %5.0f & %5.1f & %5.1f & %5.1f & %6.5e & %6.5e & %6.5e\n', m, n, k, time_lambda/repeats, iter_FBE/repeats, iter_DC/repeats, iter_major/repeats,...
      time_FBE/repeats, time_DC/repeats, time_major/repeats, fval_FBE/repeats, fval_DC/repeats, fval_major/repeats);  

end

fclose(fid);

save count10 count10