@=================================================================@
/*Translation of subroutine AMOEBA from Numerical Recipies into
   Gauss, by N.C. Mark.
   Input:  P: (ndim+1) x ndim.  Each row is a candidate parameter vector.
              The corresponding function value is found in the corresponding
              element of the vector, Y:
           Y: (ndim+1) x 1.
        FTOL: Fractional convergence tolerance to be achieved in the
              function values.

  Output:  P and Y are set to (ndim+1) new points, all within
           FTOL of a minimum function value.  Look at rows
           of P for parameter values, and corresponding elements of Y
           for the function values.

USAGE:
  Include the following proc to write out the results:

proc(0) = stmt(x);
 local k; k=rows(x)-2;
 output on;
"----------------- AMOEBA Parameters-----------------------";
 x[1:k,1];
 "Function value  :";;x[k+1,1];
 "Iterations:      ";;x[k+2,1];
 "Time in seconds: ";;(hsec-ts)/100;
 output off;
endp;

output off; ftol=1e-6; stp1=.1; maxits=1e+5;
 b=amoeba(&ofn,&strtval,ftol,stp1,maxits);
output on;
call stmt(b);

*/
@==================================================================@
@----------------------------------------------------------@
proc amoeba(&ofn,&startval,ftol,stp1,maxits) ;
local sv,np,alpha,beta,gam,iter,pr,prr,pbar,flo,ihi,inhi,i,j,
      e,rtol,ypr,yprr,x,mpts,ilo,p,y,ts ;
local ofn:proc,startval:proc   ;
@==========================================================@
/* Set up the initial simplex, an
   (np+1)x np dimensional matrix.  OUTPUT--P. The rows are the vertices
   of the starting simplex                    */
@========================================================@
@---------------------------------------------------------@
 sv=startval ;        rtol=99;      @ Initialize stuff @
 np=rows(sv) ;
 p=zeros(np+1,np) ;  ts=hsec ;
 e = zeros(np,np) ;
 y=zeros(np+1,1) ;
@---------------------------------------------------------@
 i=1 ;
 do while i <= np ;
  e[i,i] = 1 ;
  i=i+1 ;
 endo ;
 i=1 ;
 do while i <= np ;
  p[i,1:np] = sv[1:np,1]' + stp1*e[1:np,i]' ;
  i=i+1 ;
 endo ;
 p[np+1,1:np] = sv[1:np,1]' ;
@:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::@
@=====================================================================@
/* Evaluate the function at each of the vertices.
   The result is an np+1 dimensional vector, y.  The np+1's element
   corresponds to the unaltered values from startval.     */
@=====================================================================@
i=1;
do while i <= np ;
 y[i,1] = ofn(p[i,1:np]') ;
 i=i+1 ;
endo;
y[np+1,1] = ofn(sv) ;
@==================================================================@
@ Begin the body of PROC AMOEBA                                    @
@==================================================================@
alpha=1 ;                   @ parameters that define expansions @
beta=0.5;                   @ and contractions.     @
gam=2.0 ;

pr=zeros(np,1) ; prr=zeros(np,1) ;
pbar=zeros(np,1); mpts=np+1 ;

iter=0 ;
pointer1:
iter=iter+1 ;
@================================================================@
@ STATUS REPORT                                                  @
@----------------------------------------------------------------@
"----------------- AMOEBA iteration ";;iter;;"-----------------";
"parameters"; p[1,1:np] ;
"function ";;y[1,1] ;
"relative tolerance ";; rtol;
"Iteration required ";;(hsec-ts)/100;;"  seconds.";
@timehsec(ts,date)/100;;" seconds.";@
"===============================================================
                                                                ";
@====================================================================@
@ Determine which point is highest (worst), next highest, and @
@ lowest (best) by looping over points in the simplex         @
@=====================================================================@
ilo=1;
if y[1,1] > y[2,1] ;
 ihi=1; inhi=2 ;
else ;
 ihi=2; inhi=1;
endif ;
i=1;
do while i <= mpts ;
 if y[i,1] < y[ilo,1] ; ilo=i ; endif ;
 if y[i,1] > y[ihi,1] ;
  inhi=ihi ; ihi=i;
 elseif y[i,1] > y[inhi,1] ;
  if i /= ihi ; inhi = i ; endif ;
 endif ;
 i=i+1 ;
endo ;
@====================================================================@
@ Compute fractional range from highest to lowest and return if      @
@ satisfactory.                                                      @
@====================================================================@
rtol = 2*abs(y[ihi,1]-y[ilo,1])/(abs(y[ihi,1])+abs(y[ilo,1])) ;
if rtol < ftol ;
 x=p[ilo,1:np]'|y[ilo,1]|iter ; goto pointer2 ;
endif ;
if iter == maxits ;
 x=p[ilo,1:np]'|y[ilo,1]|iter ; goto pointer2 ;
endif ;
@====================================================================@
@ Begin new iteration. Compute vector average of all points except   @
@ the highest. i.e., the center of the "face" of the simplex across  @
@ from the high point.  We will subsequently explore along the ray   @
@ from the high point through that center                            @
@====================================================================@
ts=hsec ;
pbar=zeros(np,1) ;
i=1;
do while i <= mpts ;
 if i /= ihi;
  pbar=pbar + p[i,.]' ;
 endif ;
 i=i+1 ;
endo ;
@=================================================================@
@ Extrapolate by a factor alpha through the face.  i.e., Reflect the@
@ simplex from the high point                                        @
@===================================================================@
pbar=pbar/np ;
pr[.,1]=(1+alpha)*pbar-alpha*p[ihi,.]' ;
@===============================================================@
@ Evaluate the function at the reflected point                  @
@===============================================================@
ypr = ofn(pr) ;
if ypr <= y[ilo,1] ;  /* gives result better than best pt. so try an
                         additional extrapolation by a factor gam    */
 prr=gam*pr+(1-gam)*pbar ;
 yprr = ofn(prr) ;    @ check out function there   @
 if yprr < y[ilo,1] ; @ additional extrapolation succeeds, replce high pt @
  p[ihi,.]=prr[.,1]' ;
  y[ihi,1] = yprr ;
 else ;             @ additional extrapolation failed-use reflected point @
 p[ihi,.]=pr[.,1]' ;
  y[ihi,1] = ypr ;
 endif;
elseif  ypr >= y[inhi,1] ;  @reflected point worse than 2nd highest @
 if ypr < y[ihi,1] ;        @ better than highest, so replace highest @
  p[ihi,.]=pr[.,1]' ;

  y[ihi,1] = ypr ;
 endif ;              /* Look for intermediate lower point.  Perform
                         contraction of simplex along one dimension
                         then evaluate      */
 prr=beta*p[ihi,.]' + (1-beta)*pbar ;
 yprr=ofn(prr) ;
 if yprr < y[ihi,1] ;   @contraction is an improvement @
  p[ihi,.]=prr[.,1]' ;
  y[ihi,1] = yprr ;
 else ;               @ Can't get rid of high point. Better contract @
  i=1 ;               @ about low point @
  do while i <= mpts ;
  if i /= ilo ;
   j=1 ;
   do while j <= np ;
    pr[j,1] = .5*(p[i,j]+p[ilo,j]) ;
    p[i,j] = pr[j,1] ;
    j=j+1 ;
   endo ;
   y[i,1] = ofn(pr) ;
  endif ;
  i=i+1 ;
  endo ;
 endif ;
else ;       /* Arrived here if original reflection gives a middling
                point.  Replace old high point and continue  */
 p[ihi,.]=pr[.,1]' ;
 y[ihi,1] = ypr ;
endif ;
goto pointer1 ;
@======================================================================@
pointer2:
retp(x) ;
endp;