/*-------------------------------------------------------------------
I and a student (Ekaterini Kyriazidou) wrote GAUSS versions of the
optimization routines in "Numerical Recipes". The are given below and
you are welcome to use and distribute them (although it would be nice
if you thanked us if you use them in a paper - the student was paid
from an NSF grant and it is always nice if one can  point to positive
externalities of ones work). 

You should know that you use the routines at your own risk. I started
writing the routines because the GAUSS minimizers crashed when I tried to
minimize x^2.

Please let me know if you get good response from your note. I to was
frustrated with GAUSS's routines, and I would like to know if anybody else
is distributing their own.

Sincerely

Bo Honore
Associate Professor
Department of Economics
Northwestern University



THE ROUTINES:
-------------------------------------------------------------------*/

 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


NCM Note: Modified to use numerical derivatives. In the main calling
program, provide  

1. proc fct(x) 
2. proc dfct(&fct,x); local fct:proc,dt; dt=gradp(&fct,x); retp(dt'); endp;


that returns the objective function evaluated at the parameter vector
x, and a proc that computes the gradient of fct at x.

------------------------------------------------------------------------*/
    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);@
    g=dfct(&fct,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 monit(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); @ @ unnecessary function evaluation @
    dg=g;
    fret=fct(p);
@   g=dfct(p);  MO, use this if separate function for gradient is available @
     g=dfct(&fct,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);  MO, use this if separate function for gradient is available @
    xi=dfct(&fct,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 monit(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);  MO, use this if separate function for gradient is available @
    gg=g'*g;
    @ dgg=xi'*xi;    @  @ Fletcher-Reeves @
    dgg=(xi+g)'*xi;     @ Polak-Ribiere   @
    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) = AMOEBA2(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,iter,tim); /* Add this line to observe each iteration */
      if rtol<ftol;
      call monit(p,y,iter,tim);
      RETP(p,y,iter,tim);
      endif;
      if (iter==maxit);
      call monit(p,y,iter,tim);
      "Maximum number of iterations exceeded";
      RETP(p,y,iter,tim);
      endif;
      if (tim .ge maxsec);
      call monit(p,y,iter,tim);
      "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.-gam)*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,iter,tim);
cls;
     format /rd 9,4;
/*     print "iteration: ";; ? iter;print "seconds elapsed: ";; ? tim;*/
     "iteration: ";;  iter; "seconds elapsed: ";;  tim;
     print "";
     format /ro 9,6;
     print "value of function at current simplex: ";
     ? y';
     print "";
     print "current simplex (transposed): ";
     ? p';
pause(3);
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;

    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 monit(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 xit ";
    ? xi;
    ? xit;
    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);
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 100); "linmib itry=100"; 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);  */
 PROC (2) = linmin(pii,xi,tol,&fct);  @ Mod 2/7/01 YK Moh @
/*
   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 @
   pii    (nx1)    vector of starting values    @ Mod 2/7/01 YK Moh @
   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,p,xd,x1,x2,d1,d2,fx1,fx2,itry,fct:proc;
    d1=0.618; d2=0.382;
    p1=pii;
    p2=pii+xi;
    f1=fct(p1);
    itry=0;
  be:
    itry=itry+1;
    if (itry.eq 100); "linmin itry=100"; stop; 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;
/* p1'; pm'; p2'; */
    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);
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;