function [xnext,Report] = proxLineSearch(xcurrent,xpast,costfun,proxop,cost0,J,e,k,fcount,opt)
%proxLineSearch Backtracking line search for proximal gradient descent.
% costfun does not include the sparsity penalty.
% proxop: handle to the prox operator with 3 input arguments:
%     1. parameter vector x
%     2. penalty gamma
%     3. weighting vector
% s: line shrinkage factor 0<s<1
% J, e: Jacobian matrix, error vector
% cost0: output of costfun(x0), optional.
% k: outer (main) loop iteration counter

% Reference: Proximal Gradient Descent, Ryan Tibshirani
% https://www.stat.cmu.edu/~ryantibs/convexopt/lectures/prox-grad.pdf

% Copyright 2021 The MathWorks, Inc.

% todo: parameter bounds
% consider: trust region + sparsity
% puzzle: starting point: unregularized answer or initialization result?
% puzzle: projection unification
%  PCA_Tolerance: eps
%  Sparsi

s = opt.StepReduction;
Accelerate = opt.Accelerate && ~isempty(xpast);
StepSize = opt.InitialStepSize;
PreCond = opt.Precondition;
Ns = opt.NumberOfSamples;
if nargin<6
   cost0 = Ns*feval(costfun, xcurrent);
   fcount = fcount+1;
else
   cost0 = Ns*cost0;
end
StepSize = StepSize/s;
iter = 0;
wt = opt.Weight;

x0 = xcurrent;

if Accelerate
   xcurrent = x0 + (k)/(k+3)*(xcurrent-xpast);
   [~,~,e,J] = feval(costfun, xcurrent);
   fcount = fcount+1;
end

NJ = norm(J,2);
grad = J'*e;
%grad = grad/norm(grad,2); % experimental..
if PreCond
   % gauss newton (called "preconditioning" with Hessian)
   Hess = J'*J/NJ^2;
   gdir = -pinv(J)*e;
else
   % oridinary gradient descent
   Hess = 1;
   gdir = -grad;
end
% normalize gradient
% https://jermwatt.github.io/machine_learning_refined/notes/3_First_order_methods/3_9_Normalized.html
gdir = gdir/norm(gdir,2);
n0 = norm(x0,2);

lstop = 0; % successful line search
Iterate = true;
Penalized = opt.PenalizedIndices;
while Iterate
   iter = iter+1;
   StepSize = StepSize*s;
   xnextNoProx = xcurrent + StepSize*gdir;
   p = proxop(xnextNoProx(Penalized),StepSize,wt);
   xnext = xnextNoProx;
   xnext(Penalized) = p;
   %grad(Penalized) = grad(Penalized).*wt(opt.gid);
   % Note: (xnextNoProx-xnext) = lambda*xcurrent or Inf for L1
   cost = Ns*feval(costfun, xnext);
   fcount = fcount+1;
   dx = xnext-xcurrent;
   RegulBuffer = grad'*dx + 0.5*(dx'*Hess*dx)/StepSize;
   TooSmall = StepSize < 1e-12*n0;
   Iterate = (cost-cost0 > RegulBuffer) && ~TooSmall;
   if TooSmall
      lstop = 3;
   end
end

if lstop==3
   disp('no improvement')
   xnext = xcurrent;
   StepSize = 1;
   cost = cost0;
else
   %RegulBuffer
   %cost_ = Ns*feval(costfun, xnextNoProx);
   %[cost0, cost_, cost]
   disp(sum(xnext==0))
end

Report = struct('Cost',cost,'StepSize',StepSize,'Iter',iter,'Stop',lstop,...
   'x_previous',x0,'x_current',xnext,'fcount',fcount);
