function [x_soln, z_soln, iter, fval_soln, time] = p_DCAe(A, b, k, r, lambda, mu, L, opts)
% This function solves the following optimization by proximal DCA with extrapolation:
% 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

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

% Initialization
x_old = x0;
x_new = x_old;
Ax_old = Ahandle(x_old);
Ax_new = Ax_old;
theta = 1;
theta0 = 1;
iter = 1;

% best fval
x_soln = x0;
% Projection
z = Proj(Ax_new - b,r);
z_soln = z;
[~, v] = sort(abs(x_new), 'descend');
fval_soln = 1/2*norm(Ax_new - b - z)^2 + lambda*norm(x_new, 1) - mu*lambda*norm(x_new(v(1:k)),1);

while iter < maxiter
   
    % Subgradient of mu*lambda*sum_{i=1}^k |x_[i]|
    tmp = zeros(n, 1);
    tmp(v(1:k)) = mu*lambda*sign(x_new(v(1:k)));
    
    % Extrapolation
    beta = theta0/theta - 1/theta;
    u = x_new + beta*(x_new - x_old);
    Au = Ax_new + beta*(Ax_new - Ax_old);
    
    % Proximal gradient
    u_tmp = u - (1/L)*(AThandle(Au - b - z) - tmp);
    x_old = x_new;
    Ax_old = Ax_new;
    x_new = soft_thresh(u_tmp,lambda/L);
    Ax_new = Ahandle(x_new);
    
    % update best fval
    z = Proj(Ax_new - b,r);
    [~, v] = sort(abs(x_new), 'descend');
    fval = 1/2*norm(Ax_new - b - z)^2 + lambda*norm(x_new, 1) - mu*lambda*norm(x_new(v(1:k)),1);
    if fval < fval_soln
        fval_soln = fval;
        x_soln = x_new;
        z_soln = z;
    end
    % Display
    if mod(iter, freq) == 0 && disp == 1
        fprintf('iter %d, fval %g, norm_suc_x %g, alpha %g\n', iter, fval, norm(x_new - x_old), beta);
    end
    
    % Check for termination
    if sqrt((sqrt(L)*norm(Ax_new - Au) + L*norm(x_new - u))^2 + norm(x_new - x_old)^2)/max(1,norm(x_new)) < tol
        iter = iter + 1;
        break
    end
    
    % Update theta
    theta0 = theta;
    theta = (1 + sqrt(1+4*theta^2))/2;
    
    if mod(iter,200) == 0
        theta0 = 1;
        theta = 1;
    end
    
    iter = iter + 1;
end
time = toc(tstart);

fprintf('Termination: iter = %g, time = %5.2f, fval = %8.7f\n', iter, time, fval);
end

