%
% simplex - perform a simplex search
%
%   xend = simplex( func, xstart, istep, fplot, ...
%                        [ finc, fdec, [ minstep ]] )
%
% where  func    = name of objective function
%        xstart  = starting point vector
%        istep   = initial step size vector
%        fplot   = plot flag
%        finc    = step acceleration factor
%                  (default = 1, that is none)
%        fdec    = step decceleration factor
%                  (default = 1, that is none)
%        minstep = minimum step size (if deccelerating)
%
% global variable "stats" contains:
%        stats(1) = number of function evaluations
%        stats(2) = number of linear steps
%        stats(3) = number of linear reversals
%        stats(4) = number of linear searches
%
% Bob Newell, April 1996
%
%-----------------------------------------------------------
%
function xend = simplex( func, xstart, istep, ...
                            fplot, finc, fdec, minstep )
%
global stats
%
%  default parameters
%
if nargin < 6,
  finc = 1;
  fdec = 0.5;
end;
if nargin < 8,
  minstep = abs( istep ) / 50;
end;
if isempty(stats),	  
  stats = zeros( 1, 8 );
end;
%
% generate the starting simplex
%
n = length( xstart );
if n == 2,
  shape = [ 0.0 0.577 ; -0.5 -0.289 ; 0.5 -0.289 ];
else,
  disp( 'simplex.m: dimension not implemented' );
  return;
end;
step = istep;
xold = zeros( n+1, n );
fold = zeros( n+1, 1 );
for i = 1:n+1,
  xold(i,:) = xstart + shape(i,:) .* step;
  fold(i) = feval( func, xold(i,:) );
  stats(1) = stats(1) + 1;
end;
% plot it
if fold,
  for i = 1:n,
    plot( [ xold(i,1) xold(i+1,1) ], ...
	               [ xold(i,2) xold(i+1,2) ] );
  end;
  plot( [ xold(n+1,1) xold(1,1) ], ...
	               [ xold(n+1,2) xold(1,2) ] );
end;
kdiscard = 0;
%
while 1,
%
% search for worst corner
%
  k = 0;
  val = -1.0e30;
  for i = 1:n+1,
    if fold(i) > val,
	  k = i;
	  val = fold(i);
	end; 
  end;
  if k == 0,
    disp( 'simplex.m: worst not found' );
	return;
  end;
%
% if this is the one last discarded,
%
  if k == kdiscard,
    kdiscard = 0;
%
% check for convergence
%
    step = step * fdec;
	if step < minstep,
      disp( 'simplex.m: convergence' );
	  xend = xcent;
	  return;
	end;
%
% construct new simplex
%
    for i = 1:n+1,
      xold(i,:) = xcent + shape(i,:) .* step;
      fold(i) = feval( func, xold(i,:) );
      stats(1) = stats(1) + 1;
    end;
  else,
%
% replace with a reflection
%
    xdiscard = xold(k,:);
    fdiscard = fold(k);
    kdiscard = k;
%    [ kdiscard xdiscard fdiscard ]
% centroid
    xcent = zeros( 1, n );
    for i = 1:n+1,
      if i ~= k,
        xcent = xcent + xold(i,:);
      end;
    end;
    xcent = xcent / ( n );
%    [ xcent ]
% reflect
    xold(k,:) = 2 * xcent - xold(k,:);
    stats(4) = stats(4) + 1;
    fold(k) = feval( func, xold(k,:) );
    stats(1) = stats(1) + 1;
%    [ xold(k,:) fold(k) ]
  end;
% replot
  if fold,  
    for i = 1:n,
      plot( [ xold(i,1) xold(i+1,1) ], ...
	               [ xold(i,2) xold(i+1,2) ] );
    end;
    plot( [ xold(n+1,1) xold(1,1) ], ...
	               [ xold(n+1,2) xold(1,2) ] );
  end;
end;
%----------------------------------------------------------
%
