function [xo,Ot,nS] = complex_(S,x0,ip,Lb,Ub,problem,tol,mxit,edge,r,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10)
//   Unconstrained optimization using Extended Flexible Polyhedron.
//
// [xo,So,nS] = complex(S,x0,ip,Lb,Ub,problem,tol,mxit,edge,r,P1,...,P10)
//
//   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)))
//   edge: initial size of the edges (default = max(0.01*norm(x0),0.1))
//   r: extra number of vertices (default = 0)
//   xo: optimal point
//   Ot: optimal value of S
//   nS: number of objective function evaluations
//   P1,...,P10: variable arguments

//   Copyright (c) 1997-2001 by LASIM-DEQUI-UFRGS
//   $Revision: 1.0 $  $Date: 1997/05/29 10:09:20 $
//   $Revision: 2.0 $  $Date: 2001/07/05 22:20:15 $
//   Argimiro R. Secchi (arge@enq.ufrgs.br)

 if nargin < 2,
   error('complex 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(edge),
   edge=max(0.01*norm(x0),0.1);
 end
 if nargin < 10 || isempty(r),
   r=0;
 end

 x0=x0(:);
 n=size(x0,1);

// design parameters
 alfa=1;
 beta=0.5;
 gama=2;

// function conversion
 f=S;
 //if ~or(S<48)
   f=f+"(x";
   for i=1:nargin-10
      f=f+",P"+int2str(i);
   end
   f=f+")";
 //end
 f=f+"*(-problem)";
 
// polyhedron dimension
 p=n+1+r;
 N=p-1;
 u=ones(1,p);

// polyhedron building
 X=[x0 eye(n,n) rand(n,r)];
 OLDX=x0;
 x=x0;
 S=evstr(f);
 for i=2:p
    X(:,i)=x0+edge*X(:,i);
    OLDX=X(:,i-1);
    x=X(:,i);
    S(i)=evstr(f);
 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);
   ff=zeros(n1,n2);
   for i=1:n1,
    for j=1:n2,
      x=[X1(i,j);X2(i,j)];
      ff(i,j)=evstr(f);
    end
   end
   mxf=max(max(ff));
   mnf=min(min(ff));
   df=mnf+(mxf-mnf)*(2.^(([0:10]/10).^2)-1);
   contour(X1(1,:),X2(:,1),ff',df); set(gca(),"auto_clear","off");
//   clabel(v,h);
   plot(x0(1),x0(2),'ro'); h1=gca();
   plot([X(1,:) X(1,1)],[X(2,:) X(2,1)],'r');
   legend(h1,'start point');
   if ip > 0,
     disp('Pause: hit any key to continue...'); halt;
   end
 end

// start loop search
 nS=p;
 k=0;

 while k < mxit,

   k=k+1;

// find worst (h=p) and best vertices (L=1)
   [S I]=gsort(S,"g","i");
   X=X(:,I);

// centroid evaluation
   x0=sum(X(:,1:N),2)/N;

// convergence criterion
   if max(abs(S-S(1)*u')) < tol*(1+abs(S(1))) && max(max(abs(X-X(:,1)*u))) < tol*(1+norm(X(:,1))),
     break;
   end

// reflexion
   xr=x0+alfa*(x0-X(:,p));
   OLDX=x0;
   x=xr;
   Sr=evstr(f);
   nS=nS+1;

// expansion
   if Sr <= S(1)
     xe=x0+gama*(xr-x0);
     OLDX=xr;
     x=xe;
     Se=evstr(f);
     nS=nS+1;
     
     if Se < Sr
       X(:,p)=xe;
       S(p)=Se;
     else
       X(:,p)=xr;
       S(p)=Sr;
     end

// contraction
   elseif Sr > S(N)
         xc=x0+beta*(X(:,p)-x0);
         OLDX=xr;
         x=xc;
         Sc=evstr(f);
         nS=nS+1;
         
         if Sc < S(p)
           X(:,p)=xc;
           S(p)=Sc;
         else
// reduction
           for i=2:p
              X(:,i)=X(:,1)+0.5*(X(:,i)-X(:,1));
              OLDX=X(:,1);
              x=X(:,i);
              S(i)=evstr(f);
           end
             nS=nS+N;
         end
   else
     X(:,p)=xr;
     S(p)=Sr;
   end
     
   if ip && n==2,
     plot([X(1,:) X(1,1)],[X(2,:) X(2,1)],'r');
     if ip > 0,
       disp('Pause: hit any key to continue...'); halt;
     end
   end
 end

 xo=X(:,1)';
 Ot=S(1);

 if k == mxit
   disp('Warning Complex: reached maximum number of iteration');
 elseif ip && n == 2,
   plot(xo(1),xo(2),'r*'); h2=gca();
   legend([h1,h2],'start point','optimum');
 end
endfunction
 
