function [x,fval_1,iter,time] = NPG_sm(b,c1,c2,p,lambda,opts)

% This function computes a stationary point for the following
% nonconvex fused regularizer problem
%
% min 0.5*||x - b||^2 + c1*||x||_1 + c2*||Dx||_p^p, 0< p <=1
% where Dx = [x_2 - x_1,...,x_n - x_{n-1}]

% It is solved by NPG with smoothing
% min 0.5*||x - b||^2 + c1*||x||_1  + c2*||sqrt(Dx.^2 + lambda^2)||_p^p

% It uses subroutine lp_prox_mapping

% Input
%
% b           - m by 1 vector
% c1          - weight parameter for sparsity
% c2          - weight parameter for structured sparsity
% p           - 0< p <=1
% lambda      - smoothing parameter
% opts        - contains the following structure variables, default values are in brackets:
%               x0      - starting point [0]
%               tol     - tolerance [1e-5]
%               freq    - frequency for displaying output [1000]
%               maxiter - maximum number of iterations [inf]
%               disp    - display [1]: 1 - display per freq iterations
%                                      0 - no display
%               nonM       - memory in nonmonotone linesearch [5]
%
% Output
%
% x       - approximate stationary point
% fval_1  - 0.5*||x - b||^2 + c1*||x||_1
% iter    - number of iterations
% time    - cputime

% updated on April 22 to fix the bug in BB stepsize
% Further optimized on May 6, 2018

tstart = tic;
n = size(b,1); eta = 2; sigma = 1e-4; step_init = 1;

%% Parameters and defaults
if isfield(opts,'x0'),      x0      = opts.x0;      else x0 = zeros(n,1);end
if isfield(opts,'tol'),     tol     = opts.tol;     else tol = 1e-5;     end
if isfield(opts,'freq'),    freq    = opts.freq;    else freq = 1000;    end
if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = inf;  end
if isfield(opts,'disp'),    disp    = opts.disp;    else disp = 1;       end
if isfield(opts,'nonM'),    nonM    = opts.nonM;    else nonM = 5;       end

DThandle = @(z) [-z(1); -diff(z); z(end)];
x = x0;
%% Compute function value and gradient
Dx = diff(x);
fval_1 = 0.5*norm(x - b)^2 + c1*norm(x,1);
tmp = sqrt(Dx.*Dx + lambda^2);
fval_2 = c2*sum(tmp.^p);
fnew = fval_1 + fval_2;

grad = x - b + c2*p*DThandle((Dx.*(tmp.^(p-2))));

%% Set_up for nonmonotone line search
f1 = -inf(nonM,1);
f1(1) = fnew;
fmax = max(f1);

for iter = 1: maxiter
    fold = fnew;
    
    stepsize = step_init;
    xnew = x - stepsize*grad;
    xnew = lp_prox_mapping(xnew, c1*stepsize, 1); % proximal mapping
    % compute function value
    Dxnew = diff(xnew);
    fval_1 = 0.5*norm(xnew - b)^2 + c1*norm(xnew,1);
    tmp = sqrt(Dxnew.*Dxnew + lambda^2);
    fval_2 = c2*sum(tmp.^p);
    fnew = fval_1 + fval_2;
    count = 0;
    while fnew > fmax - sigma/2*norm(xnew - x)^2
        stepsize = stepsize/eta;
        xnew = x - stepsize*grad;
        xnew = lp_prox_mapping(xnew, c1*stepsize, 1); % proximal mapping
        % compute function value
        Dxnew = diff(xnew);
        fval_1 = 0.5*norm(xnew - b)^2 + c1*norm(xnew,1);
        tmp = sqrt(Dxnew.*Dxnew + lambda^2);
        fval_2 = c2*sum(tmp.^p);
        fnew = fval_1 + fval_2;
        count = count + 1;
    end
    
    %% Display
    if mod(iter,freq) == 0 && disp == 1
        fprintf(' iter = %5d, Fval_lam = %5.4e, ls_used = %3d, stepsize = %g, dsdy = %g\n',...
            iter, fnew, count, stepsize, dsdy)
    end
    
    %% Check for termination
    if norm(xnew - x)/max(1,norm(xnew)) < tol*stepsize || abs(fnew - fold)/max(1,abs(fnew)) < 1e-12
        if disp == 1
            fprintf(' TERMINATION \n')
            fprintf(' iter = %5d, Fval_lam = %5.4e, fval_1 = %5.4e\n', iter, fnew, fval_1)
            if abs(fnew - fold)/max(1,abs(fnew)) < 1e-12
                fprintf(' Terminate by lack of progress in function value\n')
            end
        end
        x = xnew;
        break
    end
    
    %% Barzilai-Borwein stepsize
    ds = xnew - x;
    gnew = xnew - b + c2*p*DThandle((Dxnew.*(tmp.^(p-2))));
    dy = gnew - grad;
    dsdy = ds(:)'*dy(:);
    if dsdy > 1e-12
        step_init = max(min(norm(ds(:))^2/(ds(:)'*dy(:)),1e8),1e-8);
    else
        step_init = max(min(2*stepsize,1e8),1e-8);
    end
    
    %% Update iterate, gradient and function value
    x = xnew;
    grad = gnew;
    f1(mod(iter,nonM)+1,1) = fnew;
    fmax = max(f1);   
end

time = toc(tstart);