function [x, fnew, iter, time] = NPG_BB(A, b, r, s, L, opts)
% This function solves the following optimization problem by nonmonotone proximal gradient with majorization,
% min 1/2d^2(Ax,D),
% s.t. ||x||_0<=s, ||x||_infty<=le8,
% D = {y: ||y||_0 <=r}+b.
 
% Input 
%
% A - m by n matrix (n>m)
% b - m by 1 vector
% r - positive integer
% s - positive integer
% L - Lipschitz contant, used for termination
% opts    - contains the following structure variables, default values are in brackets:
%           x0      - starting point [0]
%           tol     - tolerance [1e-5]
%           freq    - frequency for displaying output [5000]
%           maxiter - maximum number of iterations [inf]
%           disp    - display [1]: 1 - display per freq iterations
%                                  0 - no display
%           M       - memory in nonmonotone linesearch [5]


% Output
% 
% x    - approximate solution at termination
% iter - number of iterations
% fnew - function value at termination
% time - CPU time

tstart = tic;

% Parameters
[~,n] = size(A);

if isfield(opts,'x0'), x0 = opts.x0; else x0 = zeros(n,1); end 
if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = inf; end 
if isfield(opts,'tol'), tol = opts.tol; else tol = le-5; end 
if isfield(opts,'freq'), freq = opts.freq; else freq =5000; end 
if isfield(opts,'disp'), disp = opts.disp; else disp = 1; end 
if isfield(opts,'M'), M = opts.M; else M = 5; end

% Set matrix handles
Ahandle = @(x) A*x;
AThandle = @(x) A'*x;

% Initialization
eta = 2;
sigma = 1e-4; % line search parameter
step_init = 1; % line search parameter

x = x0;
iter = 1;

% Compute function value
tmp = Ahandle(x) - b;
tmp_1 = tmp - ProjD(tmp,r);
fnew = 0.5*norm(tmp_1)^2;

% compute gradient
g = AThandle(tmp_1);

% Set up for nonmonotone line search
f1 = -inf(M,1);
f1(1) = fnew;
fmax = max(f1);

%% Main interation
while iter <= maxiter
    
    tmpold = tmp;
    stepsize = step_init;
    z_2 = x - stepsize*g;
    xnew = ProjC(z_2,s);
    
    tmp = Ahandle(xnew) - b;
    tmp_1 = tmp - ProjD(tmp,r);
    fnew = 0.5*norm(tmp_1)^2;
    count = 0;
    while fnew > fmax - sigma/2*norm(xnew - x)^2
        % Line search
        stepsize = stepsize/eta;
        z_2 = x - stepsize*g;
        xnew = ProjC(z_2,s);
        tmp = Ahandle(xnew) - b;
        tmp_1 = tmp - ProjD(tmp,r);
        fnew = 0.5*norm(tmp_1)^2;
        count = count+1;
    end
    
    %% Display
    if mod(iter,freq) == 0 && disp == 1
        fprintf('iter = %5d, fval = %5.4e, ls_used = %3d, stepsize = %g\n', iter, fnew, count, stepsize);
    end
    
    % Brazilai-Borwein stepsize
    ds = xnew - x;
    gnew = AThandle(tmp_1);
    dy = gnew - g;
    dsdy = ds'*dy;
    if  dsdy >= 1e-12
        step_init = max(min(norm(ds)^2/dsdy,1e8),1e-8);
    else
        step_init = max(min(2*stepsize, 1e8), 1e-8); 
    end
    
    %% Check for termination
    nrm = norm(tmp - tmpold);
    if sqrt((sqrt(L)*nrm + 1/stepsize*norm(ds))^2 + norm(ds)^2)/max(1,norm(xnew)) < tol
        x = xnew;
        iter = iter + 1;
        break
    end
       
    %% Update iterate, gradient and function value
    x = xnew;
    g = gnew;
    f1(mod(iter,M)+1,1) = fnew;
    fmax = max(f1);
    iter = iter + 1;
end
time = toc(tstart);
fprintf('iter = %d, fval = %g, time = %g\n', iter, fnew, time)
