/**************  Nelder-Mead Simplex search  *************************

author: J. Huston McCulloch         mcculloch.2@osu.edu
	Economics Dept.             voice (614) 292-0382
	Ohio State Univ.            FAX   (614) 292-3906
	1945 N. High St.
	Columbus, OH 43210
        URL: http://economics.sbs.ohio-state.edu/jhm/jhm.html 
date: 10 May 1999

This code is provided gratis without performance guarantees.


Author comments:
I wrote the following "blob2" to do a Nelder-Mead Simplex search, 
generally following Numerical Recipes.   
Note that it does an _uphillsearch, not a downhill search.  
Also, for some reason I now regret slightly, 
I made the vector "x" of parameters to be found 
a row vector rather than a column vector.  
Thus, $func(x) is expecting x to be a row vector.  
Blob2 requires an initial guess x0, 
plus a vector dx0 of perturbations to x0. 

The illustration sets $func 
to a negated quadratic form with a max at the origin, and then 
sets x0 equal to the unit 4-vector.  dx0 is also a unit vector, so that
the initial simplex consists of x0 = (1 1 1 1), plus (2 1 1 1),
(1 2 1 1), (1 1 2 1), and (1 1 1 2).

I also have versions that accept a general initial simplex matrix,
and/or allow you to restrict the initial element of x with a pair
of inequalities.  I find this version to be convenient, and the 
restricted version isn't very general, so this is the best for most 
purposes.

You can roughly impose complicated restrictions 
by setting func(x) = -1.e100 at excluded x values. 

If you are trying to maximize a likelihood function 
  func(x; data), 
the data must be snuck into func(x) 
with a global variable, since blob2 is
trying to maximize wrt all of the components 
of its stated argument.
*******************************************************************/

@ demonstrates blob2 @
new;
x0 = {1  1  1  1};
dx0 = {1  1  1  1};
{xhi, yhi, its} = blob2(x0, dx0, &quad, .0001, 1);
format /rdn 20,10;
print xhi yhi its;
end;

@ Illustrative $func for use with blob@ 
proc quad(x);
retp(-x*x');
endp;

proc(3) = blob2(x0, dx0, &func, tol, view);
@ Finds 1xn xhi that maximizes func(x), by uphill simplex method.
  Returns xhi, yhi, its.
  xhi = maximizing value of x
  yhi = maximized value of func = func(xhi)
  its = # of iterations required to meet convergence criterion. 
  tol = convergence criterion, eg .0001.  This is the maximum allowed
    dimension of the simplex, scaled by dx0, and also the maximum
    allowed component of the arc-gradient of func across the simplex.
  view = 1 for on-screen progress, 0 for none, 2 for x with pauses.
  Argument x0 is 1xn row vector of initial guesses.
           dx0 is 1xn row vector of perturbations to x0.
  func must expect a row vector
  converges when dx < tol*abs(dx0) _and_ max absolute component
    of the arc gradient of func is less than tol.
  Written by Hu McCulloch, Econ Dept, Ohio State U, 
  http://economics.sbs.ohio-state.edu/jhm/jhm.html
  Posted to Gaussians 5/10/Y2K-1
  @
local func:proc;
local n, its, i1, y, i, j, ylo, yhi, centroid;
local x, dx, xnew, ynew, maxgrad;
format /rdn 10,3;

n = cols(x0);
if cols (dx0) /= n; print "blob: x0 " n "dx0 " cols(dx0); end; endif;

x = x0|(ones(n,1)*x0+diagrv(zeros(n,n),dx0));
its = 0;
i1 = 1;
y = zeros(n+1,1);

start:
i = i1;
do until i > n+1;
y[i] = func(x[i,.]);
its = its + 1;
i = i + 1;
endo;

rank:
y = rev(sortc(y~x,1));
x = y[.,2:n+1];
y = y[.,1];
ylo = y[n+1];
yhi = y[1];

if view == 2; print x~y; wait; endif;

@ test @;
@
maxgrad = inv(x[2:n+1,.] - x[1,.]) * (y[2:n+1] - y[1]);
maxgrad = maxc(abs(maxgrad)); @

if maxc((maxc(x) - minc(x))./abs(dx0')) < tol;   @ AND maxgrad < tol; @
  retp(x[1,.], yhi, its);
  endif;

@ begin reflection @
centroid = sumc(x[1:n,.])'/n;
dx = centroid - x[n+1,.];
xnew = centroid + dx;
ynew = func(xnew);
if ynew <= ylo; goto halve; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "reflect " its xnew ynew yhi; endif;
if ynew <= yhi; goto rank; endif;

@ try doubling @
xnew = centroid + 2 * dx;
ynew = func(xnew);
if ynew <= y[n+1]; goto rank; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "double  " its xnew ynew yhi; endif;
goto rank;

halve:
xnew = centroid - .5 * dx;
ynew = func(xnew);
if ynew <= ylo; goto contract; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "halve   " its xnew ynew yhi; endif;
goto rank;

contract:
x = .5 * (x[1,.] + x);
i1 = 2;
if view; print "contract" its x[1,.] yhi yhi; endif;
goto start;

endp;

