function [w] = lrtrain

global x
global y
global C
global P
global eps

max_iter = 1000;

% parameters of trust region algorithm
eta0 = 1e-3 ; eta1 = 0.25 ; eta2 = 0.75 ;
sigma1 = 0.25 ; sigma2 = 0.5; sigma3 = 4.0 ;

w = zeros(size(x,2),1);
f = fun(w); g = grad(w);
init_gnorm = norm(g);
delta = norm(g) ;
   
iter = 1;
while iter<=max_iter,
  [s, cg_iter] = trcg(w, delta, g);
  s  = s .* P;
  w_new = w + s;
  fnew = fun(w_new) ;
  actred =  f - fnew ;  
      
  % Line search
  gs = g'*s ;
  theta=1;
  nu = 0.01;
  while(actred<nu*theta*gs)
      theta = theta * 0.5;
      w_new = w - theta * s;
      fnew = fun(w_new);
      actred = f - fnew;
  end;
 
  liblr_eps = 0.01*min(length(y(y>0)), length(y(y<0))) / length(y);
  if (norm(g) < liblr_eps * init_gnorm)
      stopped = 1;
  else
      stopped = 0;
  end
  fprintf(1, 'iter %2d f %5.16e max|g| %5.5e CG %3d stopped %1d \n', iter, f, norm(g,inf), cg_iter, stopped);

  iter = iter + 1 ; 
  w = w_new ; 
  f = fnew;
  g = grad(w) ;
  if (max(abs(g)) < eps)
    fprintf(1, 'iter %2d f %5.16e max|g| %5.5e CG %3d stopped %1d \n', iter, f, norm(g,inf), cg_iter, stopped);
    break ;
end
end

function [f] = fun(w)

global x
global y
global C

f = w'*w/2 + C*sum(log(1+exp(-y.*(x*w))));

function [g] = grad(w)

global x
global y
global C
global D
global P

z=1./(1+exp(-y.*(x*w)));
g = w + C*(((z-1).*y)'*x)'; 
D = z.*(1-z);
P = 1;


function [Hv] =  Hv(s)

global x
global D
global C
global P

Hv = s + C*((D.*(x*s))'*x)';

function [s, cg_iter] = trcg(w, delta, g)

global y
global x
global D
global P

s = zeros(size(g));
g_hat = P.*g;
cgtol = 0.1*norm(g_hat);


r = -g_hat;
d = r;

cg_iter = 0;
rTr = r'*r;
while 1,
  if norm(r) <= cgtol,
    break
  end
  cg_iter = cg_iter + 1;
  Hd = Hv(d.*P).*P;
  alpha = rTr/(d'*Hd);
  s = s + alpha*d;
  r = r - alpha*Hd;
  rnewTrnew = r'*r;
  beta = rnewTrnew/(rTr);
  d = r + beta*d;
  rTr = rnewTrnew;
end

