[xo,Ot,nS]=newton(S,x0,ip,G,H,Lb,Ub,problem,tol,mxit)

function [xo,Ot,nS]=newton(S,x0,ip,G,H,Lb,Ub,problem,tol,mxit)
%   Unconstrained optimization using Newton.
%
%   [xo,Ot,nS]=newton(S,x0,ip,G,H,Lb,Ub,problem,tol,mxit)
%
%   S: objective function
%   x0: initial point
%   ip: (0) no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
%   G: gradient vector function
%   H: Hessian matrix function
%   Lb, Ub: lower and upper bound vectors to plot (default = x0*(1+/-2))
%   problem: (-1): minimum (default), (1): maximum
%   tol: tolerance (default = 1e-4)
%   mxit: maximum number of iterations (default = 50*(1+4*~(ip>0)))
%   xo: optimal point
%   Ot: optimal value of S
%   nS: number of objective function evaluations

%   Copyright (c) 2001 by LASIM-DEQUI-UFRGS
%   $Revision: 1.0 $  $Date: 2001/07/16 08:10:12 $
%   Argimiro R. Secchi ([e-mail])

 if nargin < 2,
   error('newton requires two input arguments!');
 end
 if nargin < 3 | isempty(ip),
   ip=0;
 end
 if nargin < 4 | isempty(G),
   grad=0;
 else
   grad=1;
 end
 if nargin < 5 | isempty(H),
   hes=0;
   if ~grad,
     error('Newton: numerical gradient not allowed with numerical Hessian!');
   end
 else
   hes=1;
 end
 if nargin < 6 | isempty(Lb),
   Lb=-x0-~x0;
 end
 if nargin < 7 | isempty(Ub),
   Ub=2*x0+~x0;
 end
 if nargin < 8 | isempty(problem),
   problem=-1;
 end
 if nargin < 9 | isempty(tol),
   tol=1e-4;
 end
 if nargin < 10 | isempty(mxit),
   mxit=50*(1+4*~(ip>0));
 end

 x0=x0(:);
 y0=feval(S,x0)*problem;
 n=size(x0,1);
 
 if ~grad,
   d=1e-5*abs(x0)+1e-5*ones(n,1);
   CHG=0.01*d;
   gr=d;
 end
 
 if ip & n == 2,
   figure(abs(ip));
   [X1,X2]=meshgrid(Lb(1):(Ub(1)-Lb(1))/20:Ub(1),Lb(2):(Ub(2)-Lb(2))/20:Ub(2));
   [n1,n2]=size(X1);
   f=zeros(n1,n2);
   for i=1:n1,
    for j=1:n2,
      f(i,j)=feval(S,[X1(i,j);X2(i,j)]);
    end
   end
   mxf=max(max(f));
   mnf=min(min(f));
   df=mnf+(mxf-mnf)*(2.^(([0:10]/10).^2)-1);
   [v,h]=contour(X1,X2,f,df); hold on;
%   clabel(v,h);
   h1=plot(x0(1),x0(2),'ro');
   legend(h1,'start point');

   if ip > 0,
     disp('Pause: hit any key to continue...'); pause;
   end
 end
 
 xo=x0;
 yo=y0;
 it=0;
 nS=1;
 
 while it < mxit,
   if ~grad,
     y=yo;
     % Finite difference perturbation levels
     % First check perturbation level is not less than search direction.
     j=find(10*abs(CHG)>abs(d));
     CHG(j)=-0.1*d(j);
     % Ensure within user-defined limits
     for i=1:n
        xo(i)=xo(i)+CHG(i);
        yo=feval(S,xo)*problem;
        gr(i)=(yo-y)/CHG(i);
        if yo > y,
          y=yo;
        else
          xo(i)=xo(i)-CHG(i);
        end
     end
     % Try to set difference to 1e-8 for next iteration
     % Add eps for machines that can't handle divide by zero.
     CHG=1e-8./(gr+eps); 
     yo=y;
     nS=nS+n;
   else
     gr=feval(G,xo)*problem;
     gr=gr(:);
   end

  if ~hes,
    if problem == 1,
      fun={S,G};
    else
      fun={inline(['-feval(''' S ''',x)']),inline(['-feval(''' G ''',x)'])};
    end;
    
    Hstr=sparse(ones(n));        % dense matrix only
    group=color(Hstr,(n+1)*ones(n,1)-colmmd(Hstr)');
    he=full(sfd(xo,gr,Hstr,group,[],[{['fun_then_grad']},{['newton']},fun]));
  else
    he=feval(H,xo)*problem;
  end
  
  d=-he\gr;
  stepsize=1;
  x=xo+d;
  yo=feval(S,x)*problem;
  nS=nS+1;

  while yo < y0,
    stepsize=-0.5*stepsize;
    x=xo+stepsize*d;
    yo=feval(S,x)*problem;
    nS=nS+1;
  end

  xo=x;
  it=it+1;

  if norm(xo-x0) < tol*(0.1+norm(x0)) & abs(yo-y0) < tol*(0.1+abs(y0)),
    break;
  end
  
  if ip & n == 2,
    plot([x0(1) xo(1)],[x0(2) xo(2)],'r');
    if ip > 0,
      disp('Pause: hit any key to continue...'); pause;
    end
  end 

  x0=xo;
  y0=yo; 
 end
 
 Ot=yo*problem;
 
 if it == mxit,
   disp('Warning Newton: reached maximum number of iterations!');
 elseif ip & n == 2,
   h2=plot(xo(1),xo(2),'r*');
   legend([h1,h2],'start point','optimum');
 end