function [x,time,iter] = proj_ls(A,b,r,k,epsilon,type,tol,x0,maxiter,disp,freq)
% This code solves the following problem using grad proj with nonmonotone
% linesearch:
%
% min sum f( dist^2(x,Di)) s.t. ||x||_0 <= k, ||x||_infty <= r
% 
% The Di takes the form {x : ai'*x <= bi}
%
% Inputs
% A        - its rows are ai'
% b        - collection of bi
% r        - bounds on ||x||_infty
% k        - bounds on the cardinality of x
% epsilon  - parameters in f
% type     - choice of f: 'frac', 'log' or '', corresponding to fraction
%            penalty, log penalty and f(t) = t
% tol      - termination tolerance
% x0       - random initial point. Does not have to be feasible
% maxiter  - maximum number of iterations
% disp     - {0} no display; {1} display iterations
% freq     - the frequency of display
%
% Outputs
% x        - approximate solution
% time     - run-time
% iter     - number of iterations

tstart = tic;
[m,n] = size(A);
M = 10; % nonmonotone line search parameters
c = 1e-4; % nonmonotone line search parameters

if strcmp(type,'frac')
    f = @(x) x./(x + epsilon) + epsilon*x; % sum f gives the function value
    g = @(x) epsilon./((x + epsilon).^2) + epsilon;
elseif strcmp(type,'log')
    f = @(x) log(x/epsilon + 1);% - log(epsilon);
    g = @(x) 1./(x + epsilon);
else
    f = @(x) x/m;
    g = @(x) ones(m,1)/m;
end

normrA2 = sum(A.^2,2); % vector containing the norm squares of rows of A

tmp = x0; %%%%%%%%%%%%%% Initial point
tmp0 = max(min(tmp,r),-r);
tmp2 = tmp.^2 - (tmp0 - tmp).^2;
[~,I] = sort(tmp2,'descend');
x = zeros(n,1);
x(I(1:k)) = tmp0(I(1:k));
% function value
temp = (A'*spdiags(max(A*x - b,0)./normrA2,0,m,m)); % x - proj x
t = sum(temp.^2,1);
fnew = sum(f(t)); 
% compute grad
utemp = g(t');
grad = 2*(temp*utemp);
% initial stepsize
step_init = 1;

iter = 1;
f1 = -inf(M,1);
f1(1) = fnew;
fmax = max(f1);
neg_cur = 0; % keep track of negative curvature

while iter <= maxiter
    
    % line search
    stepsize = step_init;
    
    tmp = x - stepsize*grad;
    tmp0 = max(min(tmp,r),-r);
    tmp2 = tmp.^2 - (tmp0 - tmp).^2;
    [~,I] = sort(tmp2,'descend');
    xnew = zeros(n,1);
    xnew(I(1:k)) = tmp0(I(1:k));
    temp = (A'*spdiags(max(A*xnew - b,0)./normrA2,0,m,m)); % x - proj x
    t = sum(temp.^2,1);
    fnew = sum(f(t));
    while fnew > fmax - c/2*norm(x - xnew)^2
        stepsize = stepsize/2;
        tmp = x - stepsize*grad;
        tmp0 = max(min(tmp,r),-r);
        tmp2 = tmp.^2 - (tmp0 - tmp).^2;
        [~,I] = sort(tmp2,'descend');
        xnew = zeros(n,1);
        xnew(I(1:k)) = tmp0(I(1:k));
        temp = (A'*spdiags(max(A*xnew - b,0)./normrA2,0,m,m)); % x - proj x
        t = sum(temp.^2,1);
        fnew = sum(f(t));
    end
    % check for termination
    if norm(x - xnew) <= tol*max(norm(xnew),1) % heuristic to prevent premature convergence in one step
        x = xnew;
        break
    end
    
    if disp == 1 && mod(iter,freq) == 0
        fprintf(' iter %d fval %g err %g\n', iter, sum(f(t)), norm(x - xnew))
    end
    
    % BB stepsize
    utemp = g(t');
    gnew = 2*(temp*utemp); 
    dg = gnew - grad;
    dx = xnew - x;
    inner = dx'*dg;
    if inner > 1e-12
        step_init = max(min(norm(dx)^2/inner,1e10),1e-10);
    else
        step_init = max(min(2*stepsize,1e10),1e-10);
%         fprintf(' Warning: negative curvature. Iter %d \n', iter)
        neg_cur = neg_cur + 1;
    end
    
    % update
    x = xnew;
    grad = gnew;
    f1(mod(iter,M)+1,1) = fnew;
    fmax = max(f1);
    
    iter = iter + 1;
   
end
fprintf(' Termination: iter %d fval %g neg_cur %g\n', iter, fnew, neg_cur)
time = toc(tstart);

