/************* OPTIMIZATION ROUTINES  ******************************
The following are GAUSS versions of the optimization routines in
"Numerical Recipes" written by Bo Honore and Ekaterini Kyriazidou
of Northwestern University, with support from an NSF grant.
You are welcome to use and distribute them as long as you
include proper attribution to the authors.
You should know that you use the routines at your own risk.
THE ROUTINES FOLLOW:

***********************************************************/
/* corrected 9/9/96 */

PROC (4) = dfpmin(p,ftol,maxsec,maxit,&fct,&dfct,&linm);

/* This procedure minimizes a function called FCT using the
   Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.The procedure is tailored
   after DFPMIN in Numerical Recipes.

   INPUT

   p      (nx1)    starting vector
   ftol   (1x1)    fractional convergence tolerance in the function value
   maxsec (1x1)    maximum number of seconds allowed for running program
   maxit  (1x1)    maximum number of iterations allowed
   fct    proc     function we want to minimize
   dfct   proc     routine calculating  the gradient of the function
   linm   proc     routine performing line minimization

   OUTPUT

   p      (nx1)    minimizing vector
   fret   (1x1)    value of function at minimium,p
   iter   (1x1)    iterations taken
   tim    (1x1)    seconds elapsed
*/
    local tim,date1,fp,xi,g,hessin,hdg,dg,iter,eps,fret,n,tol,fac,fad,fae,
      fct:proc,dfct:proc,linm:proc;

    tim=0;
    date1=date;
    eps=0.0000000001;
    tol=0.00001;
    n=rows(p);
    fp=fct(p);
    g=dfct(p);
    xi=-g;
    hessin=eye(n);
    iter=1;
    do while iter<=maxit;
    {p,fret}=linm(p,xi,tol,&fct);
    tim=ethsec(date1,date)/100;
    call monitor(p,fret,iter,tim);
    if ((2*abs(fp-fret)).le ftol*(abs(fp)+abs(fret)+eps));
    RETP(p,fret,iter,tim);
    endif;
    if (tim .ge maxsec); "Maximum number of seconds exceeded";
    RETP(p,fret,iter,tim);
    endif;
    fp=fct(p);
    dg=g;
    fret=fct(p);
    g=dfct(p);
    dg=g-dg;
    hdg=hessin*dg;
    fac=dg'*xi;
    fae=dg'*hdg;
    if (fac==0); RETP(p,fret,iter,tim); else; fac=1/fac; endif;
    if (fae==0); RETP(p,fret,iter,tim); else; fad=1/fae; endif;
    dg=fac*xi-fad*hdg;
    hessin=hessin+fac*xi*xi'-fad*hdg*hdg'+fae*dg*dg';
    xi=-hessin*g;
    iter=iter+1;
    endo; "Maximum number of iterations exceeded";
    RETP(p,fret,iter,tim);
ENDP;
PROC (4) = frprmn(p,ftol,maxsec,maxit,&fct,&dfct,&linm);

/* This procedure minimizes a function called FCT using the conjugate
   gradients method of Fletcher-Reeves-Polak-Ribiere. The procedure is
   tailored after FRPRMN in Numerical Recipes.

   INPUT

   p      (nx1)    starting vector
   ftol   (1x1)    fractional convergence tolerance in the function value
   maxsec (1x1)    maximum number of seconds allowed for running program
   maxit  (1x1)    maximum number of iterations allowed
   fct    proc     function we want to minimize
   dfct   proc     routine calculating  the gradient of the function
   linm   proc     routine performing line minimization

   OUTPUT

   p      (nx1)    minimizing vector
   fret   (1x1)    value of function at minimium,p
   iter   (1x1)    iterations taken
   tim    (1x1)    seconds elapsed
*/

    local tim,date1,fp,xi,g,h,gg,dgg,gam,iter,eps,fret,n,tol,
       fct:proc,dfct:proc,linm:proc;

    tim=0;
    date1=date;
    eps=0.0000000001;
    tol=0.00001;
    n=rows(p);
    fp=fct(p);
    xi=dfct(p);
    g=-xi;
    h=g;
    xi=h;
    iter=1;
    do while iter<=maxit;
    {p,fret}=linm(p,xi,tol,&fct);
    tim=ethsec(date1,date)/100;
    if (tim .ge maxsec);
    "Maximum number of seconds exceeded";
    RETP(p,fret,iter,tim);
    endif;
    call monitor(p,fret,iter,tim);
    if ((2*abs(fp-fret)).le ftol*(abs(fp)+abs(fret)+eps));
    RETP(p,fret,iter,tim);
    endif;
    fp=fct(p);
    xi=dfct(p);
    gg=g'*g;
    dgg=(xi+g)'*xi;
    if (gg==0); RETP(p,fret,iter,tim); endif;
    gam=dgg/gg;
    g=-xi;
    h=g+gam*h;
    xi=h;
    iter=iter+1;
    endo;
    "Maximum number of iterations exceeded";
    RETP(p,fret,iter,tim);
ENDP;
PROC (4) = AMOEBA(p,ftol,maxsec,maxit,&fct);

/* This procedure minimizes a function called FCT using the simplex
   method. The procedure is tailored after AMOEBA in Numerical Recipes.
   We are minimizing over an NDIM-dimensional vector of parameters.
   Input is a matrix P, whose NDIM+1 rows are NDIM-dimensional vectors
   which are the vertices of the starting simplex.
   FTOL is the fractional convergence tolerance to be achieved in the function
   value.
   ALP, BET, GAM below, are parameters which define the expansions and
   contractions.

   INPUT

   p      ((ndim+1) x ndim)   starting simplex
   maxsec (1x1)               maximum number of seconds allowed
   maxit  (1x1)               maximum number of iteraions allowed
   ftol   (1x1)               fractional convergence tolerance
   fct    proc                function we want to minimize

   OUTPUT

   p      ((ndim+1) x ndim)   final simplex; first row is the minimizing
                  vector
   y      ((ndim+1) x 1)      vector of values of fct at final simplex;
                              first number is the value of the function at
                              the minimum
   iter   (1x1)               number of iterations taken
   tim    (1x1)               number of seconds of running time

   While running the following are printed on the screen:
   the number of the current iteration and of seconds of running time,
   the value of the function at the current simplex (transposed), and
   the current simplex (transposed).
   In the end the following are printed on the screen:
   the number of the final iteration, the number of seconds taken,
   the value of the function at the final simplex (transposed); the
   first number is the value of the function at the minimum
   the final simplex (transposed); first column is the minimizing vector.
*/

     local y,j,date1,tim,ndim,npts,pr,prr,pbar,ind,ihi,inhi,ilo,rtol,
       alp,bet,gam,ypr,yprr,i,iter,fct:proc   ;

     tim=0;
     date1=date;
     alp=1.; bet=0.5; gam=2.0;
     ndim=cols(p);
     npts=ndim+1;
     y=zeros(ndim+1,1);
     j=1;
       do while j<=npts;
       y[j,1]=fct(p[j,.]');
       j=j+1;
       endo;


     iter=0;
     begy:
      tim=ethsec(date1,date)/100;
      ind=sortind(y);
      ihi=ind[npts,1];
      inhi=ind[npts-1,1];
      ilo=ind[1,1];
      rtol=2.*abs(y[ihi,1]-y[ilo,1])/abs(y[ihi,1]+y[ilo,1]);
     call monit(p,y,tim,iter);  /* Add this line to observe each iteration */
      if rtol<ftol;
      call monit(p,y,tim,iter);
      RETP(p,y,iter,tim);
      endif;
      if (iter==maxit);
      call monit(p,y,tim,iter);
      "Maximum number of iterations exceeded";
      RETP(p,y,iter,tim);
      endif;
      if (tim .ge maxsec);
      call monit(p,y,tim,iter);
      "Maximum number of seconds exceeded";
      RETP(p,y,iter,tim);
      endif;
      iter=iter+1;
      pbar=(sumc(p)-p[ihi,.]')/ndim;
      pr=(1.+alp)*pbar-alp*p[ihi,.]';
      ypr=fct(pr);
      if ypr <=y[ilo,1];
     prr=gam*pr+(1.-gam)*pbar;
     yprr=fct(prr);
     if yprr < y[ilo,1];
    p[ihi,.]=prr';
    y[ihi,1]=yprr;
     else;
    p[ihi,.]=pr';
    y[ihi,1]=ypr;
     endif;
      elseif ypr >= y[inhi,1];
     if ypr < y[ihi,1];
    p[ihi,.]=pr';
    y[ihi,1]=ypr;
     endif;
     prr=bet*p[ihi,.]'+(1.-bet)*pbar;
     yprr=fct(prr);
     if yprr < y[ihi,1];
    p[ihi,.]=prr';
    y[ihi,1]=yprr;
     else;
    p=.5*(p+p[ilo,.]);
    i=1;
    do while i<=npts;
       y[i,1]=fct(p[i,.]');
       i=i+1;
    endo;
     endif;
      else;
     p[ihi,.]=pr';
     y[ihi,1]=ypr;
      endif;
      goto begy;

ENDP;

PROC (0) = monit(p,y,tim,iter);

     format /rd 9,2;
     print "iteration: "; ? iter;print "seconds elapsed: "; ? tim;
     print "";
     format /rd 11,3;
     print "value of function at current simplex: ";
     ? y';
     print "";
     print "current simplex (transposed): ";
     ? p';

ENDP;
PROC (4) = powell(p,ftol,maxsec,maxit,&fct,&linm);

/* This procedure minimizes a function called FCT using Powell's
   method. The procedure is tailored after POWELL in Numerical Recipes.

   INPUT

   p      (nx1)      starting value
   ftol   (1x1)      fractional tolerance in the function value
   maxsec (1x1)      maximum number of seconds allowed for running program
   maxit  (1x1)      maximum number of iterations allowed
   fct    proc       function we want to minimize
   linm   proc       line minimization routine

   OUTPUT

   p     (nx1)      minimizing vector
   fret  (1x1)      value of function at minimum,p
   iter  (1x1)      number of iterations taken
   tim   (1x1)      number of elapsed seconds
*/

    local  date1,tim,fret,pt,iter,n,i,xi,xit,tol,ptt,fptt,t,ibig,del,fp,
           fct:proc,linm:proc;

" starting";
    tim=0;
    date1=date;
    tol=0.00001;
    n=rows(p);
    xi=eye(n);
    fret=fct(p);
    pt=p;
    iter=0;
begy:
    tim=ethsec(date1,date)/100;
    iter=iter+1;
    call monitor(p,fret,iter,tim);
    fp=fret;
    ibig=0;
    del=0;
    i=1;
    do while i<=n;
    xit=xi[.,i];
    {p,fret}=linm(p,xit,tol,&fct);
    if (abs(fp-fret).gt del);
    del=abs(fp-fret);
    ibig=i;
    endif;
    i=i+1;
    endo;
    if ((2*abs(fp-fret)).le ftol*(abs(fp)+abs(fret)));
    RETP(p,fret,iter,tim);
    endif;
    if (iter.eq maxit);
    "Maximum number of iterations exceeded";
    RETP(p,fret,iter,tim);
    endif;
    if  (tim.ge maxsec);
    "Maximum number of seconds exceeded";
    RETP(p,fret,iter,tim);
    endif;
    ptt=2*p-pt;
    xit=p-pt;
    pt=p;
    fptt=fct(ptt);
    if (fptt.ge fp); goto begy; endif;
    t=2*(fp-2*fret+fptt)*(fp-fret-del)^2-del*(fp-fptt)^2;
    if (t.ge 0); goto begy; endif;
    {p,fret}=linm(p,xit,tol,&fct);
    xi[.,ibig]=xit;
    goto begy;
ENDP;
PROC(4)=genalg(bst,sig,m2,sig2,ga,maxit,maxsec,param,&fct,&calp);

/*
This procedure minimizes a function called FCT using a "Generic Algorithm".

    INPUT

    bst     (kx1)   vector of starting values
    sig     (1x1)   variance of initial draws
    m2      (1x1)   half the number of draws
    sig2    (1x1)   variance of change
    ga      (1x1)   fraction changed
    maxit   (1x1)   maximum number of iterations allowed
    maxsec  (1x1)   maximum number of seconds allowed
    param   (2x1)   vector indicating minimization options:
                    param[1,1]=0:  cross-over by changing coordinates
                              =1:  cross-over by random averaging
                    param[2,1]=0:  re-use best
                              =1:  best not automatically re-used.
    fct     proc    function to be minimized
    calp    proc    calculates probabilities from f.
                    Two alternative procs are provided in the end of this
                    program: calp_fu and calp_ra.

   OUTPUT

   b        (kx1)   minimizing vector
   f        (1x1)   function of value at minimum
   iter     (1x1)   number of iterations taken
   tim      (1x1)   seconds elapsed

  While running, the following are printed in each row on the screen:
  number of iteration;seconds elapsed;current value of function;current
  parameter vector.
*/


local date1,iii,g,gg,f,ii,aa,p,ij,i,m,k,ii1,ii2,mix,iter,tim,
      fbest,gbest,m1,t2,fct:proc,calp:proc;

"START";
date1=date;
k=rows(bst);
m1=m2*2;
if (param[2].eq 0); m=m1+1; else; m=m1; endif;
iii=seqa(1,1,m);
g=bst+sig*rndn(k,m-1);
iter=1;
gbest=bst; fbest=fct(gbest);
g=bst~g;
xxlbl:

f=fct(g);
ii=minindc(f);
tim=ethsec(date1,date)/100;
format /rd 5,0;
iter;;
format /rd 7,2;
tim;;
format /rd 9,4;
fbest~gbest';
if (fbest.ge f[ii]); fbest=f[ii]; gbest=g[.,ii]; endif;
if (iter.ge maxit);
    "Maximum number of iterations exceeded";
    RETP(gbest,fbest,iter,tim);
endif;
if (tim.ge maxsec);
    "Maximum number of seconds exceeded";
    RETP(gbest,fbest,iter,tim);
endif;
gg=g[.,ii];
p=calp(f);
p=cumsumc(p);
i=1;
if (param[1].eq 0);
ij=int(rndu(m2,1)*(k-1))+1;
else;
mix=rndu(k,m2);
endif;
do while i<=m2;
ii1=rndu(1,1);
ii1=counts(ii1,p)'*iii;
ii2=rndu(1,1);
ii2=counts(ii2,p)'*iii;
if (param[1].eq 0);
g=g~(g[1:ij[i],ii1]|g[ij[i]+1:k,ii2])
   ~(g[1:ij[i],ii2]|g[ij[i]+1:k,ii1]);
else;
g=g~(g[.,ii1].*mix[.,i]+g[.,ii2].*(1-mix[.,i]))
   ~(g[.,ii2].*mix[.,i]+g[.,ii1].*(1-mix[.,i]));
endif;
i=i+1;
endo;
g=g[.,m+1:m+2*m2];
aa=rndu(k,2*m2);
aa=(aa.lt ga);
g=g+aa.*(sig2*rndn(k,2*m2));
if (param[2].eq 0); g=gbest~g; endif;
iter=iter+1;
goto xxlbl;
ENDP;

/* CALP PROCS */

PROC(1)=calp_ra(f);
local aa;
aa=(rows(f)+1-rankindx(f,1));
RETP(aa/sumc(aa));
ENDP;

PROC(1)=calp_fu(f);
local aa;
aa=-f+maxc(f)+0.0001;
RETP(aa/sumc(aa));
ENDP;
PROC (2) = linmib(pst,xi,tol,&fct);

/* This routine performs line minimization of a (n-dimensional) function
   called FCT, using Brent's one-dimensional minimization method without
   derivatives. It is tailored after BRENT in Numerical Recipes, starting
   off with an ad hoc initial bracketing of the minimum. Given as inputs
   a vector P and a direction vector XI, it minimizes FCT(P+X*XI) w/r/t
   the scalar X, isolating X (the best abscissa) at the midpoint of two
   scalars, a & b, that are 2*X*TOL apart, i.e. at a fractional precision
   of about TOL. The maximum allowed number of iterations is 100.

   INPUT

   pst    (nx1)     starting value
   xi     (nx1)     direction vector
   tol    (1x1)     fractional tolerance
   fct    proc      function we want to minimize

   OUTPUT

   pst+x*xi    (nx1)     minimizing vector along direction xi; new vector
                         direction
   fx          (1x1)     value of fct at minimum
*/

    local maxit,cgold,zeps,a1,a2,ad,am,f1,f2,fm,itry,
          v,u,w,fu,fv,fw,x,fx,tol1,tol2,e,etemp,iter,q,r,a,b,xm,d,p,fct:proc;

    cgold=0.3819660;  maxit=100; zeps=.0000000001;
    a1=0;
    a2=0.01;
    f1=fct(pst+a1*xi);
    itry=0;
  be:
    itry=itry+1;
    if (itry.eq 10); "linmib itry=10"; stop; endif;
    f2=fct(pst+a2*xi);
    if(f1.eq f2); a2=a2+1; goto be; endif;
    if (f1.le f2);
    am=a1; fm=f1;
    f1=f2; a1=a2;
    else;
    am=a2; fm=f2;
    endif;
    ad=5*(am-a1);
    a2=am+ad;

    f2=fct(pst+a2*xi);
    do until (f2.gt fm);
    ad=2*ad;
    a2=am+ad;
    f2=fct(pst+a2*xi);
    endo;
    if (a1.<a2);
    a=a1; b=a2;
    else;
    a=a2; b=a1;
    endif;
    v=am; x=v; w=v;
    e=0;

    fx=fct(pst+x*xi); fv=fx; fw=fx;


    iter=1;
    do while (iter <=maxit);
    xm=(a+b)*0.5;
    tol1=tol*abs(x)+zeps;
    tol2=2*tol1;
    if (abs(x-xm).le (tol2-0.5*(b-a))); goto lb3; endif;
    if (abs(e).gt tol1);
        r=(x-w)*(fx-fv);
        q=(x-v)*(fx-fw);
        p=(x-v)*q-(x-w)*r;
        q=2*(q-r);
        if (q.gt 0); p=-p; endif;
        q=abs(q);
        etemp=e;
        e=d;
        if ((abs(p).ge abs(0.5*q*etemp)).or (p.le q*(a-x))
        .or (p.ge (q*(b-x)))); goto lb1; endif;
        d=p/q;
        u=x+d;
        if (((u-a).lt tol2).or ((b-u).lt tol2));
        d=abs(tol1)*((xm.>x)-(xm.<x));
        endif;
        goto lb2;
    endif;
      lb1:
    if (x.ge xm);
        e=a-x;
    else;
        e=b-x;
    endif;
    d=cgold*e;
      lb2:
    if (abs(d).ge tol1);
        u=x+d;
    else;
        u=x+abs(tol1)*((d.>0)-(d.<0));
    endif;
    fu=fct(pst+u*xi);
    if (fu.le fx);
        if (u.ge x); a=x; else; b=x; endif;
        v=w; fv=fw; w=x; fw=fx; x=u; fx=fu;
    else;
        if (u.lt x); a=u; else; b=u; endif;
        if ((fu.le fw).or (w.eq x));
        v=w; fv=fw; w=u; fw=fu;
        elseif ((fu.le fv).or (v.eq x).or (v.eq w));
        v=u; fv=fu;
        endif;
    endif;
    iter=iter+1;
    endo;
    "BRENT EXCEEDED MAXIMUM NUMBER OF ITERATIONS";
  lb3:

RETP(pst+x*xi,fx);
ENDP;
PROC (2) = linmin(p,xi,tol,&fct);

/*
   This routine performs an ad hoc n-dimensional Golden Section Search for
   the minimum of a function called FCT.

   INPUT

   pi    (nx1)    vector of starting values
   xi    (nx1)    direction vector
   tol   (1x1)    tolerance
   fct   proc     function we want to minimize

   OUTPUT

   p     (nx1)    minimizing vector
   f     (nx1)    value of fct at minimum,p
*/

    local f1,f2,fm,f,p1,p2,pm,xd,x1,x2,d1,d2,fx1,fx2,itry,fct:proc;
    d1=0.618; d2=0.382;
    p1=p;
    p2=p+xi;
    f1=fct(p1);
    itry=0;
  be:
    itry=itry+1;
    if (itry.eq 10); "linmin itry=10 (probably a problem)";
       RETP(p1,f1); endif;
    f2=fct(p2);
    if(f1.eq f2); p2=p2+xi; goto be; endif;
    if (f1.le f2);
    pm=p1; fm=f1;
    f1=f2; p1=p2;
    else;
    pm=p2; fm=f2;
    endif;
    xd=pm-p1;
    p2=pm+xd;
    f2=fct(p2);
    do until (f2.gt fm);
    xd=2*xd;
    p2=pm+xd;
    f2=fct(p2);
    endo;
    stit:
    x1=d1*p1+d2*p2;
    x2=d2*p1+d1*p2;
    fx1=fct(x1);
    fx2=fct(x2);
    if (fx1.le fx2);
    p2=x2;
    else;
    p1=x1;
    endif;
    if (sqrt(sumc((p1-p2).^2)).ge tol); goto stit; endif;
    if (fx1.lt fx2);
    p=x1;
    else;
    p=x2;
    endif;
    f=fct(p);
    if (fm.lt f); f=fm; p=pm; endif;
RETP(p,f);
ENDP;
PROC (4)=lnsrch(p1,p2,n,&fct);

/*
This routine computes the value of a function called FCT at N, equally
spaced, points on the line connecting two vectors P1 & P2.
It returns the minimizing vector (the endpoints included), the value of the
function at that point, as well as those vectors that yield a value equal
or less to the value of the function at max(p1,p2), along with their
respective values.

INPUT

p1  (rx1)   endpoint
p2  (rx1)   endpoint
n   (1x1)   number of equally spaced points
fct proc    function we are minimizing

OUTPUT

PPPPP'   (rx?)  minimizing vector(s?)
z        (1x1)  value of function at minmium
pppp'    (rx?)  vectors with function values <= fct(max(p1,p2))
x        (?x1)  vector of function values at pppp'
*/


    local i,r,x,y,z,m,p,pp,ppp,pppp,ppppp,ff,f1,f2,f,s,ss,fct:proc;

    f1=fct(p1);f2=fct(p2);
    r=rows(p1);
    p=zeros(r,n+1);
    f=zeros(n+1,1);
    i=1;
    pp=zeros(r,n-1);
    ff=zeros(n-1,1);
    do while i<n;
    pp[.,i]=p1+i*(p2-p1)/n;
    ff[i,1]=fct(pp[.,i]);
    i=i+1;
    endo;
    p=p1~pp~p2;
    f=f1|ff|f2;
    y=(f[.,1] .>f1).and (f[.,1] .>f2);
    x=delif(f,y);
    ppp=p';
    s=sumc(y);
    pppp=zeros(n+1-s,r);
    pppp=delif(ppp,y);

    z=minc(f);
    m=f[.,1] .==z;
    ss=sumc(m);
    ppppp=zeros(ss,r);
    ppppp=selif(ppp,m);

RETP(ppppp',z,pppp',x);
ENDP;


PROC (0) = monitor(p,y,iter,tim);
     format /rd 9,4;
     print "iteration: ";? iter;
     print "";
     format /rd 16,6;
     print "current value of function: ";? y;
     print "";
     format /rd 10,3;
     print "current parameter vector: " ;? p';
     print "";
     print "seconds elapsed:" ;? tim;
ENDP;