function [x_soln, z_soln, iter, fval_soln, time] = NPG_DC(A, b, k, r, lambda, mu, L, opts)
% This function solves the following optimization by nonmonotone proximal gradient?
% min 0.5*||Ax - b - z||^2 + lambda*||x||_1 - mu*lambda*sum_{i=1}^k |x_[i]|,
% s.t.||z||_0 <= r,
% where x_[i] is the i-th largest entry in x in magnitude.


% Input
%
% A  - m by n matrix (n > m)
% b  - m by 1 vector measurement
% k  - positive integer <= n
% r  - positive integer
% lambda  - real number > 0
% mu  - real number in (0, 1)
% L  - Lipschitz constant



% Ouput
% x_soln,z_soln - solution with the best function value
% iter          - number of iterations
% fval_soln     - best function value
% 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 = 1e-4;     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;
step_init = 1;

x = x0;
iter = 1;

% Compute function value
tmp = Ahandle(x) - b;
z = Proj(tmp,r);
[~, v] = sort(abs(x), 'descend');
fval = 1/2*norm(tmp - z)^2 + lambda*norm(x, 1) - mu*lambda*norm(x(v(1:k)),1);
g = AThandle(tmp - z);

% save best fval
x_soln = x0;
fval_soln = fval;
z_soln = z;

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

while iter < maxiter
    tmp_old = tmp;
    stepsize = step_init;
    % Subgradient of mu*lambda*sum_{i=1}^k |x_[i]|
    v_tmp = zeros(n, 1);
    v_tmp(v(1:k)) = mu*lambda*sign(x(v(1:k)));
    dir = g - v_tmp;
    
    u_tmp = x - stepsize*dir;
    x_new = soft_thresh(u_tmp, lambda*stepsize);
    
    tmp = Ahandle(x_new) - b;
    z_new = Proj(tmp,r);
    [~, v] = sort(abs(x_new), 'descend');
    fnew = 1/2*norm(tmp - z_new)^2 + lambda*norm(x_new, 1) - mu*lambda*norm(x_new(v(1: k)),1);
    count = 0;
    while fnew > fmax - sigma/2*norm(x_new - x)^2
        % Line search
        stepsize = stepsize/eta;
        u_tmp = x - stepsize*dir;
        x_new = soft_thresh(u_tmp, lambda*stepsize);
        
        tmp = Ahandle(x_new) - b;
        z_new = Proj(tmp,r);
        [~, v] = sort(abs(x_new), 'descend');
        fnew = 1/2*norm(tmp - z_new)^2 + lambda*norm(x_new, 1) - mu*lambda*norm(x_new(v(1: k)),1);
        count = count + 1;
    end
    
    % Display
    if mod(iter, freq) == 0 && disp == 1
        fprintf('iter %d, fval %g, norm_suc_x %g, ls_used %g, step_init %g\n', iter, fnew, norm(x_new - x), count, step_init);
    end
    
    % Barzilai-Borwein stepsize
    s = x_new - x;
    gnew = AThandle(tmp - z_new);
    y = gnew - g;
    dsdy = s'*y;
    if  dsdy >= 1e-12
        step_init = max(min(norm(s)^2/dsdy,1e8),1e-8);
    else
        step_init = max(min(2*stepsize, 1e8), 1e-8); % code changed for r=r*1.1, BB fails happens
%         fprintf('Warning: BB fails, iter = %d\n', iter);
    end
    
    % update best fval
    if fnew < fval_soln
        x_soln = x_new;
        fval_soln = fnew;
        z_soln = z_new;
    end
    
    
    % Check for termination
    nrm = norm(tmp - tmp_old);
    if sqrt((sqrt(L)*nrm + 1/stepsize*norm(s))^2 + norm(s)^2)/max(1,norm(x_new)) < tol
        iter = iter + 1;
        break
    end
    
    % Update maximal previous function value
    f1(mod(iter,M)+1,1) = fnew;
    x = x_new;
    g = gnew;
    fmax = max(f1);
    iter = iter + 1;
    
end
time = toc(tstart);
fprintf('Termination: iter = %g, time = %5.2f, fval = %8.7f\n', iter, time, fnew);
end