function [x,fval_1,iter,time] = SDCAM(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 SDCAM with subproblem
% min 0.5*||x - b||^2 + c1*||x||_1  + inf_y{1/(2*lambda)||y - Dx||^2 + c2*||y||_p^p}
% solved by proximal DCA with nonmonotone line search
% 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

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 = ones(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);
prox = lp_prox_mapping(Dx, c2*lambda, p); % proximal mapping
fval_1 = 0.5*norm(x - b)^2 + c1*norm(x,1);
fval_2 = norm(prox - Dx)^2/(2*lambda) + c2*norm(prox,p)^p;
fnew = fval_1 + fval_2;

g1 = x - b + DThandle(Dx)/lambda;
g2 = DThandle(prox)/lambda;
%% Set_up for nonmonotone line search
f1 = -inf(nonM,1);
f1(1) = fnew;
fmax = max(f1);

%% Main iteration
for iter = 1: maxiter
    fold = fnew;
    
    stepsize = step_init;
    xnew = x - stepsize*(g1 - g2);
    xnew = lp_prox_mapping(xnew, c1*stepsize, 1); % proximal mapping
    % Compute function value
    Dxnew = diff(xnew);
    prox = lp_prox_mapping(Dxnew, lambda*c2, p); % proximal mapping
    fval_1 = 0.5*norm(xnew - b)^2 + c1*norm(xnew,1);
    fval_2 = norm(prox - Dxnew)^2/(2*lambda) + c2*norm(prox,p)^p;
    fnew = fval_1 + fval_2;
    count = 0;
    while fnew > fmax - sigma/2*norm(xnew - x)^2
        stepsize = stepsize/eta;
        xnew = x - stepsize*(g1 - g2);
        xnew = lp_prox_mapping(xnew, c1*stepsize, 1); % proximal mapping
        % Compute function value
        Dxnew = diff(xnew);
        prox = lp_prox_mapping(Dxnew, lambda*c2, p); % proximal mapping
        fval_1 = 0.5*norm(xnew - b)^2 + c1*norm(xnew,1);
        fval_2 = norm(prox - Dxnew)^2/(2*lambda) + c2*norm(prox,p)^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\n', iter, fnew, count, stepsize)
    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;
    gnew1 = xnew - b + DThandle(Dxnew)/lambda;
    dy = gnew1 - g1;
    step_init = max(min(norm(ds(:))^2/(ds(:)'*dy(:)),1e8),1e-8);
    
    
    
    %% Update iterate, gradient and function value
    x = xnew;
    g1 = gnew1;
    g2 = DThandle(prox)/lambda;
    f1(mod(iter,nonM)+1,1) = fnew;
    fmax = max(f1);    
end
time = toc(tstart);