function [zbest,fbest,dval,dfeas,iter,time] = mFBS_logistic_flasso_uv(A,lambda1,lambda2,tol,maxiter,freq)

% This function uses the variant of Tseng's MFBS for solving
% fused lasso logistic regression.
%
% min sum(log (1 + exp(A*z))) + lambda1*norm(z(1:end-1),1) + lambda2*norm(diff(z(1:end-1)),1)
%
% where we define M*z to be [z(1:end-1);diff(-z(1:end-1))],
% with the dual problem given by
%
% max sum(nu*log(nu) + (1-nu)*log(1-nu)) 
% s.t. A'*nu + M'*y = 0, norm(y(1:n-1),'inf') <= lambda1,
%      norm(y(n:2*n-3),'inf')\le lambda2.
%
% We recommend using this code only when m <= 500 so that inverting A*A' is
% not prohibitively expensive
%
% Input
% A       - m by n matrix, with the first n-1 columns coming from data, the
%           last column is a {-1,1} vector representing the classifier.
% lambda1 - regularization for l1 penalty
% lambda2 - regularization for the fused lasso penalty
% tol     - termination tolerance
% maxiter - maximum number of iterations allowed
% freq    - frequency of display
%
% Output
% zbest   - approximate optimal solution at termination
% fbest   - approximate optimal value at z at termination
% dval    - approximate dual optimal value at termination
% dfeas   - dfeas infeasibility by the dual solutions (violation of the equality)
% iter    - number of iterations
% time    - total cpu time; measured by tic-toc
%
% Written by Ting Kei Pong
% Last modified on XXX

tstart = tic;
% Set parameters
AAT = A*A';
tmp = max(eig(AAT))/4; % use eigs if AAT is huge
L = 0.5*(tmp + sqrt(tmp^2 + 20))/0.95; 

% Initialization
n = size(A,2);
z = zeros(n,1);
y = zeros(2*n-3,1);
v = y;
iter = 0;

barz = zeros(n,1);
bary = zeros(2*n-3,1);
fbest = inf;
zbest = barz;

% main loop
while iter <= maxiter
  
  % compute v
  gradv = Mmap(-z,1);
  tmp = y - gradv/L;
  v(1:n-1) = max(min(tmp(1:n-1),lambda1),-lambda1);
  v(n:end) = max(min(tmp(n:end),lambda2),-lambda2);
  
  % compute u
  tmp = exp(A*z); 
  x = A'*(tmp./(1 + tmp)); % gradient of h at z
  gradu = x + Mmap(y,2);
  u = z - gradu/L;
  
  % update y
  y = v - (Mmap(-u,1) - gradv)/L;
  
  % update z 
  tmp = exp(A*u); 
  x = A'*(tmp./(1 + tmp)); % gradient of h at u
  z = u - ( (x + Mmap(v,2)) - gradu )/L;
  
  iter = iter + 1;  
  barz = u; % much slower with averaging
  bary = v;
  
  % compute function values and check for termination
  if mod(iter,freq) == 0
    
    fval = sum(log(1 + exp(A*barz))) + lambda1*norm(barz(1:end-1),1) + lambda2*norm(diff(barz(1:end-1)),1);
    if fval < fbest % find the smallest primal value up to the current iteration
      fbest = fval;
      zbest = barz;
    end
    
    tmp2 = Mmap(bary,2);
    testx = AAT\(-(A*tmp2)); % use y to solve for nu; can be expensive if AAT is huge.
    if min(testx) >= 0 && max(testx) <= 1 % if so, use this nu
      nu = testx;
    else
      tmp = exp(A*barz);
      barx = A'*(tmp./(1 + tmp)); % gradient of h at barz
      nu = AAT\(A*barx); % solve for nu for dual function evaluation; can be expensive if AAT is huge.
    end
    dval = sum(nu.*log(nu) + (1-nu).*log(1-nu));
    tmp1 = A'*nu; 
    dfeas = norm(tmp1 + tmp2)/max([norm(tmp1),norm(tmp2),1]);
    gap = abs(fbest + dval)/max(fbest,1);
  
  
%     fprintf(' iter = %d, fval = %g, dval = %g, dfeas = %g\n', iter, fval, dval, dfeas)
  
    if gap < tol && dfeas < tol/5
      break
    end
    
  end
end
time = toc(tstart);
fprintf(' Termination: iter = %d, fval = %g, dval = %g, dfeas = %g, time = %g\n', iter, fval, dval, dfeas, time)