function [xo,Ot,nS]=rosembr(S,x0,ip,Lb,Ub,problem,tol,mxit,stp,amp,red)
%   Unconstrained optimization using Rosembrock.
%
%   [xo,Ot,nS]=rosembr(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 = 3)
%   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 22:20:15 $
%   Argimiro R. Secchi (arge@enq.ufrgs.br)

 if nargin < 2,
   error('rosembr 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=3;
 end
 if nargin < 11 | isempty(red) | red <= 0 | red >= 1,
   red=0.5;
 end

 x0=x0(:);
 yo=feval(S,x0)*problem;
 n=size(x0,1);
 D=eye(n);
 
 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'); h2=0;
   legend(h1,'start point');
   if ip > 0,
     disp('Pause: hit any key to continue...'); pause;
   end
 end
 
 xo=x0;
 Ot=yo;
 it=0;
 nS=1;
 
 while it < mxit,
                     % exploration
  delta=zeros(n,1);   
  idx=zeros(n,1);    % index of insuccess directions
  c=0;
                     % exit loop only if all directions (+/-) failed at least once,
                     % and at least one direction had success at least once
  while ~all(idx-(idx~=0)) | ~all(idx-c),
   for i=1:n,
     x=xo+stp(i)*D(:,i);
     y=feval(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 h2 == 0, legend(h1,'start point');
         else legend([h1,h2],'start point','new direction');
         end
         if ip > 0,
           disp('Pause: hit any key to continue...'); pause;
         end
       end

       xo=x;
       delta(i)=delta(i)+stp(i);
       stp(i)=stp(i)*amp;
     else
      stp(i)=-stp(i)*red;
      idx(i)=idx(i)+1;
     end
   end
   c=c+1;
  end

  it=it+1;
 
  if norm(x-x0) < tol*(0.1+norm(x0)) & abs(y-Ot) < tol*(0.1+abs(Ot)),
    break;
  end
                  % progression
  Ot=yo;
  x0=xo;

  if ip & n == 2,
    h2=plot(xo(1),xo(2),'kh');
    if it==1,
      legend([h1,h2],'start point','new direction');
    end
  end
  
  if any(delta-flipud(delta)),
    [dd,map]=sort(abs(delta));
    map=flipud(map);
    delta=delta(map);
    D=D(:,map');
    stp=stp(map);
  end
  
  L=triu(ones(n,1)*delta');
  A=D*L';
  m=min(find(abs(delta)<eps));
  if isempty(m),
    m=n;
  else
    m=m-1;
  end
  
  for i=1:m,
    a=A(:,i);
    b=a;
    for j=1:i-1
      b=b-a'*D(:,j)*D(:,j);
    end
    D(:,i)=b/norm(b);
  end
 end
 
 Ot=yo*problem;
 
 if it == mxit,
   disp('Warning Rosembr: reached maximum number of stages!');
 elseif ip & n == 2,
   h3=plot(xo(1),xo(2),'r*');
   legend([h1,h2,h3],'start point','new direction','optimum');
 end
 