/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/* simplex algorithm for minimizing ff(b) call as

        {bend,vff,cnv} =simplex(&ff,b0,prt)
where
        b0: starting values
        ff: procedure returning value of fnction to be minimized
        prt: printing pointer,
                if prt = 1, intermediate steps will be printed
                if prt = 2, intermediate steps will be printed
                            with parameter estimates and gradients
        bend: best parameter vector found
        vff: value of function at bend
        cnv: if converged = 1
*/

/*::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
proc (3) = simplex(&ff,start,prt);
local n,fall,startall,i,nstart,nstart2,nstart3,faf,iter,c,fnew,fnew2,
fnew3,rtol,gam,db,maxiter,alpha,beta,ftol,cnv,kk,mask,fmt,pdummy;
local ff:proc;
output off;
maxiter = 500000;        /*maximum number of iterations*/
alpha = 1;            /*step size parameters            */
beta=2;               /*alpha > 0  beta > 1, 1 > gam > 0*/
gam = 0.5;
ftol = 0.0000001~0.0000001;
rtol = 0;
cnv = 0;
kk = 0;
let mask[1,3] = 1 1 1;
let fmt[3,3] = "*.*lG" 5 4 "*.*lf" 16 6 "*.*lf" 16 6;
format /rd 10,4;

nst:
n = rows(start);
fall = ff(start);   /* first function evaluation*/
startall = start;  /* startall vector of starting values*/
i = 1;
do until i > n;
    nstart = start;
    if start[i,.] .ne 0;nstart[i,.]=start[i,.]*1.1;
        else;nstart[i,.]=0.1;
    endif;
    fall = fall~ff(nstart);
    startall = startall~nstart;
    i = i + 1;
endo;
                        /*startall = nx(n+1) matrix of paramter vectors*/
              /*fall = 1x(n+1) corresponding vector of function values*/

iter = 1;
do until iter > maxiter;
    if (iter % 20000) == 0;
    if prt > 0;cls;locate 1,1;? "restarting";endif;
       start = startall[.,1];
       startall = start;
       fall = ff(start);
       i = 1;
       do until i > n;
           nstart = start;
           if start[i,.] .ne 0;nstart[i,.]=start[i,.]*1.1;
               else;nstart[i,.]=0.1;
           endif;
           fall = fall~ff(nstart);
           startall = startall~nstart;
           i = i + 1;
       endo;
    endif;
                                   /* sort according to function value */
    faf = (sortc(fall'~startall',1));
    fall = faf[.,1]'; startall = faf[.,2:cols(faf)]';
    /* smallest value fall[.,1], corresponding vector startall[.,1] */
    kk = key;
    if kk == 1120;goto eiter;endif;             /*ALT+1 stops program*/
    if kk == 1121;start=startall[.,1];goto nst;endif;/*ALT+2 rstarts*/

     /* check for convergence*/
     if minc(abs(startall[.,1])) > 0;
     rtol =maxc(abs((fall[.,1]-fall[.,2]))|
        abs((startall[.,1]-startall[.,2])/*./startall[.,1]*/));
     if rtol < ftol[.,2];cnv = 1;goto eiter;endif;
     else;
     rtol =abs((fall[.,1]-fall[.,2]));
     if rtol < ftol[.,1];cnv = 1;goto eiter;endif;
     endif;
/*

     rtol =abs((fall[.,1]-fall[.,2]));
     if rtol < ftol[.,1];cnv = 1;goto eiter;endif;
*/
    if prt==1;locate 2,1;
        ? "iter      rtol      function";
        pdummy = printfm(iter~rtol~fall[.,1],mask,fmt);
    endif;
    if prt==2;locate 2,1;
        ? "iter      rtol      function";
        pdummy = printfm(iter~rtol~fall[.,1],mask,fmt);
        locate 5,1;
        if minc(abs(startall[.,1])) > 0;
        ? " parameter estimates and psuedo gradients ";
        seqa(1,1,rows(startall))~startall[.,1]~((startall[.,1]-startall[.,2])./startall[.,1]);
        else;
        "parameter estimates";
        seqa(1,1,rows(startall))~startall[.,1];
        endif;
    endif;
    c = meanc(startall[.,1:n]');
    nstart = c + (c-startall[.,n+1])*alpha;
    fnew = ff(nstart);

    if (fnew > fall[.,1]) .and (fnew < fall[.,n]);
    if prt > 0;locate 1,1;? "condition 1";endif;
        startall = nstart~startall[.,1:n];
        fall = fnew~fall[.,1:n];
        goto niter;
    endif;

    if fnew <= fall[.,1];
    if prt > 0;locate 1,1;? "condition 2";endif;

        nstart2 = c + (nstart-c)*beta;
        fnew2 = ff(nstart2);
            if fnew2 < fnew;
                startall = nstart2~startall[.,1:n];
                fall = fnew2~fall[.,1:n];
            else;
                startall = nstart~startall[.,1:n];
                fall = fnew~fall[.,1:n];
            endif;
        goto niter;
    endif;

    if fnew >= fall[.,n];
    if prt > 0;locate 1,1;? "condition 3";endif;
        if fnew >= fall[.,n+1];
            nstart3=(1-gam)*startall[.,1]+gam*startall[.,n+1];
        else;
            nstart3=(1-gam)*startall[.,1]+gam*nstart;
        endif;
        fnew3 = ff(nstart3);
        if fnew3 <= minc(fnew|fall[.,n+1]);
            startall = nstart3~startall[.,1:n];
            fall = fnew3~fall[.,1:n];
        else;
            i = 2;
            do until i > n;
               nstart3 = (startall[.,i]+startall[.,1])*0.5;
               fall[.,i]=ff(nstart3);
               startall[.,i] = nstart3;
               i = i + 1;
            endo;
        endif;
        goto niter;
    endif;

    niter:
iter = iter + 1;
endo;
eiter:
cls;
output on;
retp(startall[.,1],fall[.,1],cnv);
endp;

