function [xo,Ot,nS]=hkjeeves(S,x0,ip,Lb,Ub,problem,tol,mxit,stp,amp,red)
//   Unconstrained optimization using Hooke & Jeeves.
//
//   [xo,Ot,nS]=hkjeeves(S,x0,ip,Lb,Ub,problem,tol,mxit,stp,amp,red)
//
//   S: objective function
//   x0: initial point
//   ip: (0): no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
//   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 stages (default = 50*(1+4*~(ip>0)))
//   stp: stepsize vector for the independent variables (default = max(0.01*abs(x0+~x0),0.1))
//   amp: stepsize enlargement factor (1,oo) (default = 1.5)
//   red: stepsize reduction factor (0,1) (default = 0.5)
//   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/05 21:10:15 $
//   Argimiro R. Secchi (arge@enq.ufrgs.br)

//   direction = - 1 --> reverse
//                 1 --> direct
//   idx: variables index vector with enlarged stepsize
//   top: end of idx list
//   bottom: beggin of idx list
//   next: index of enlarged variable
 
 if nargin < 2,
   error('hkjeeves requires two input arguments');
 end
 if nargin < 3 || isempty(ip),
   ip=0;
 end
 if nargin < 4 || isempty(Lb),
   Lb=-x0-~x0;
 end
 if nargin < 5 || isempty(Ub),
   Ub=2*x0+~x0;
 end
 if nargin < 6 || isempty(problem),
   problem=-1;
 end
 if nargin < 7 || isempty(tol),
   tol=1e-4;
 end
 if nargin < 8 || isempty(mxit),
   mxit=50*(1+4*~(ip>0));
 end
 if nargin < 9 || isempty(stp),
   stp=max(0.01*abs(x0+~x0),0.1);
 else
   stp=abs(stp(:));
 end
 if nargin < 10 || isempty(amp) || amp <= 1,
   amp=1.5;
 end
 if nargin < 11 || isempty(red) || red <= 0 || red >= 1,
   red=0.5;
 end

 x0=x0(:);
 yo=evstr(S+"(x0)")*problem;
 n=size(x0,1);

 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
 
 x=x0;
 xo=x0;
 Ot=yo;
 it=0;
 nS=1;
 top=0;
 bottom=0;
 idx=zeros(n+1,1);
 idx(bottom+1)=top;
 
 while it < mxit,
  next=bottom;
  norma=0;
                     // exploration
  for i=1:n,
   stp_i = stp(i);
   
   for direction=[1 -1],
     x(i)=xo(i)+stp_i*direction;
     y=evstr(S+"(x)")*problem;
     nS=nS+1;
     
     if y > yo,     // success
       yo=y;
       if ip && n == 2,
         plot([x(1) xo(1)],[x(2) xo(2)],'r');
         if ip > 0,
           disp('Pause: hit any key to continue...'); halt;
         end
       end

       xo(i)=x(i);
       idx(next+1)=i;
       next=i;
       break;
     end
   end  
  
   x(i)=xo(i);
   norma=norma+stp_i*stp_i/(x(i)*x(i)+(abs(x(i))<tol));
  end
   
  it=it+1;

  if sqrt(norma) < tol && abs(yo-Ot) < tol*(0.1+abs(Ot)),
    break;
  end
                  // progression
  if next==bottom,
    stp=stp*red;
  else
    good=1;
    idx(next+1)=top;

    while good,
      Ot=yo;
    
      next=idx(bottom+1);
      while next ~= top,
        x(next)=x(next)+amp*(x(next)-x0(next));
        x0(next)=xo(next);
        next=idx(next+1);
      end
    		  
      if idx(bottom+1) ~= top,
        y=evstr(S+"(x)")*problem;
        nS=nS+1;
      
        if y > yo,
          yo=y;
          good=1;
        
          if ip && n == 2,
            plot([x(1) xo(1)],[x(2) xo(2)],'r');
            if ip > 0,
              disp('Pause: hit any key to continue...'); halt;
            end
          end
        else
          good=0;
        end

        next=idx(bottom+1);
        while next ~= top,
          if good,
            xo(next)=x(next);
          else
            x(next)=xo(next);
          end
          next=idx(next+1);
        end
      end
    end
  end
 end

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

endfunction
