-- Perform multidimensional function minimization.
~ The user must supply:

"f(parms)"  The function to be minimized where
            "parms" is an array of parameters.
"guess"     An array  of initial parameter values.
"deltas"    Initial step size for the parameters.
            If the parameters have different
            scale sizes, "deltas" should define an
            array with a step size for each,
            otherwise a scalar can be used.
"tol"       The termination criteria based on the
            amount of change in parameter values.
            A tol of 1.0e-6 gives about 3 digit
            accuracy.

The routine "minimize" steps until a minimum is found.
After minimize is run the following globals are set:

"minp"  The solution parameter array
"p"     The final simplex vertices
"y"     The function values at each vertex

~
------------ Example --------------------
-- fit a curve to data by minimizing the sum of the square of the differences between data and fit.

fit(a,x) = a[1]+a[2]*cos(a[3]*x)

data=read(xydata); xx[i]=data[i,1]; yy[i]=data[i,2] 

err(a) = sum((fit(a,xx[j])-yy[j])^2,j,1,count(data))

-- set up globals for minimize
f(a) = err(a)
guess = {80,-75,1}
deltas = .5
tol = 1e-5

plot data

minimize:;    minp:{80.2,-73.4,0.895}
plot fit(minp,X)

----------- minimization routines -------
-- These routines use the "downhill simplex method". A simplex is a N dimensional geometrical figure with N+1 vertices. The algorithim evaluates the function at each vertex and steps the simplex toward the minimum by: 1) reflecting a vertex away from the function's highest point. 2) reflection and expansion away from the highest point. 3) contraction away from the highest point. 4) contraction in all dimensions toward the lowest point.

-- step until parameters change by less than tol
minimize = init,
           step,
           step while pchange > tol,
           minp:=psum/m

pchange = sum((abs(p[lowest]-p[highest])/
           (abs(p[lowest])+abs(p[highest])))[ii],ii,
              1,n)

-- move simplex one step
step =  rank,
 try(-alpha),                -- reflect from highest
 try(gamma) when ytry  y[lowest], -- keep going
 (ysave:=y[highest],
  try(beta),                 -- contract from highest
  shrink when ytry  ysave) when ytry  y[high]

-- find lowest, highest and next highest vertex
rank = lowest:=1,high:=1,highest:=2,i:=1,
   (lowest:=i when y[i]<y[lowest],,
    (high:=highest,highest:=i) when y[i]>y[highest],,
    high:=i when y[i]>y[high] and i  highest,,
    i:=i+1) while im

-- move highest vertex through face by fac
try(fac) = ptry:=psum*(1-fac)/n-
                  p[highest]*((1-fac)/n-fac),
    ytry:=f(ptry),
    (y[highest]:=ytry,
     psum:=psum+ptry-p[highest],
     p:=replace(p,highest,ptry)) when ytry<y[highest]

-- shrink entire simplex
shrinkp(k)[j] = .5*(p[k,j]+p[lowest,j]) dim[n]
shrink = kk:=1,      
  ((psum:=shrinkp(kk),
    p:=replace(p,kk,psum)) when kklowest,,
    kk:=kk+1 ) while kkm,
  y:=fp

-- define finite arrays so they can be assigned
sump[j] = sum(p[ii,j],ii,1,m) dim[n]
initp[i,j]= guess[j]+(deltas[j] when i=j,0) dim[m,n]
fp[i]=f(p[i]) dim[m]

init = n:=count(guess), m:=n+1,
       p:=initp,
       psum:=sump,
       y:=fp

replace(array,index,newval)[i] = newval when i=index,
                                 array[i]

-- tuneable constants
alpha=1; beta=.5; gamma=2
