function [x, fnew, iter, time] = NPG_constant(A, b, r, s, L, opts)
% This function solves the following optimization problem by 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 and constructing the constant
%     stepsize
% 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
%

% x    - soulution with the best function value
% iter - number of iterations
% fnew - best function value
% time - CPU time

tstart = tic;
L_0 = L + 1e-4; % reciprocal of the constant stepsize

% 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-4; end 
if isfield(opts,'freq'), freq = opts.freq; else freq =5000; end 
if isfield(opts,'disp'), disp = opts.disp; else disp = 1; end 


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

% Initialization
x = x0;
iter = 1;

% Compute function value
tmp = Ahandle(x) - b;
tmp_1 = tmp - ProjD(tmp,r);

% compute gradient
g = AThandle(tmp_1);


%% Main interation
while iter <= maxiter
    
    xold = x;
    tmpold = tmp;
    z_2 = x - 1/L_0*g;
    x = ProjC(z_2,s);
    
    tmp = Ahandle(x) - b;
    tmp_1 = tmp - ProjD(tmp,r);
    
    
    %% Display
    if mod(iter,freq) == 0 && disp == 1
        fnew = 0.5*norm(tmp_1)^2;
        fprintf('iter = %5d, fval = %5.4e\n', iter, fnew);
    end
    
    ds = x - xold;  
    
    %% Check for termination
    nrm = norm(tmp - tmpold);
    if sqrt((sqrt(L)*nrm + L_0*norm(ds))^2 + norm(ds)^2)/max(1,norm(x)) < tol
        iter = iter + 1;
        break
    end    
    
    %% Update iterate and gradient    
    g = AThandle(tmp_1); 
    iter = iter + 1;
    
  
end
fnew = 0.5*norm(tmp_1)^2;
time = toc(tstart);
fprintf('iter = %d, fval = %g, time = %g\n', iter, fnew, time)
