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 (arge@enq.ufrgs.br)

 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=evstr(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)=evstr(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);
   contour(X1(1,:),X2(:,1),f',df); set(gca(),"auto_clear","off");
//   clabel(v,h);
   plot(x0(1),x0(2),'ro'); h1=gca();
   legend(h1,'start point');

   if ip > 0,
     disp('Pause: hit any key to continue...'); halt;
   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=evstr(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=evstr(G+"(xo)")*problem;
     gr=gr(:);
   end

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

  while yo < y0,
    stepsize=-0.5*stepsize;
    x=xo+stepsize*d;
    yo=evstr(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...'); halt;
    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,
   plot(xo(1),xo(2),'r*'); h2=gca();
   legend([h1,h2],'start point','optimum');
 end
endfunction
