
/**********************************************
  GRBL2_MLE.INC  -- maxlik alternative (procedures for maximum likelihood estimation)

Procedures here are:  maxlikc & compcvc & figstuff.

( & Other procs:
    mloptsgt figmaxst makestep goldenc  askmlalg grblbtch
    btchbnds askbound invpd2 setflu mlinvpd max_prc1 max_prc2 max_prc3
    max_prc4 gt_tech mk_techs do_tmph1 figfreev dofreev figfree0
    do_base2 binmxst  helptech mle_constraints askconst shwconst)


  Maximum likelihood estimation, ala DMH.

  Features:
        1) Will do constrained optimization, given a set of linear constraints
           Constraints consists of a set (of any size) of either equality or
           inequality contraints, with inequality constraints incorporating any
           number of parameters.
        2) Will switch between a variety hill climbing algorithims, if 'stuck'.
           Allows a large amount of control of when to use different
           algorithims (e.g.; at specific iterations, at every nth iteration,
           only between certain iterations), and which to use first at any
           given iteration.

Many of the bells and whistles of the MAXLIK procedures are NOT
supported -- although a certain amount of "changing strategy on the fly" is
permitted.

 ********* The user is responsible for displaying the results.  *********


 To detect exiting via A :
        if likelihood=missing, then A was hit.

Currently, to detect exiting via C, check value of _forcecvg (or check
the gradient vector stored in _MLEGRAD to tell whether an optimum
has been reached).


Basic algorithim is, where the -loglik is being minimized:
              Note: B2= B1 - s*Inv(D)*G
                     s = step
                     H = Direction matrix
                     G = Gradient

--------- MAXLIKC   ------
    {li,beta,freei}=
         maxlikc(b0,cnstlist,hessalg,dogolden,prnt,header,&model,parnames);

Outputs:
   li = Log Likelihood at maximum. Missing if "aborted"   (scalar)
   beta  = coefficient vector      (Kx1)
   freei = Points into beta -- the coefficients that are NOT constrained.

GLOBAL OUTPUTS:
    _mlegrad = Gradient (summed) at maximum (a Kx1 vector).
    _mlehess = Hessian (BHHH or NR, or quasi-Newton) at maximum (KxK)
    _mliter = Actual # of iterations required.
    _mlebuse = Current estimate of beta coefficients (useful if you want to examine current values)

    _CtCall = A three element counter of how many times the MODEL has been called (the model may be called
	      several times in a single iteration). 
	       _CtCall[1]= total number of calls, across all iterations (including this one)
		      you can use this to initialize stuff at the beginning of the  1st iteration.
	       _CtCall[2]= number of call in this iteration (including current one)
		      you can use this to initialize stuff at the beginning of each iterations.
	       _CtCall[3]= total number of iterations (including the current one)

   
Inputs:
 b0 =  A Kx1 vector of starting values.  BETA will
         have the same number of rows as b0.
 cnstlist = A list containing constraint information.  See mle_constraints for a
           handy routine to create this list from an ascii string containing
           constraint info.
      NOTE: The constraint list may also contain equality constraints.
            Equality constrained variables are NOT maximized over
            (no coefficient or std. errors computed)

 hessalg  = Algorithim to use for 'direction' matrix.
             1= Use 2nd derivs  = H
             2= Use inner product of gradient = G'G
             3 = Steepest descent
             4= DFP
             5= BFGS
             6= Reserved for User-method (e.g. -- GridSearch)
             7= Random
 dogolden  = Set to 1 if GOLDEN step search desired.
  prnt     = Set to 1 to print intermediate results
  header  =  Header displayed during estimation (a string). If 0,
                nothing will be displayed.
 model     = Name of proc that will compute gradient, likelihood, and
             'direction' matrix for entire dataset.  This proc should be of
              form: {li,grad,dirmtx}=model(b,arg1)
                        b=beta
                        arg1 == Compute:
                           0) Li,G and DirMtx=G'G
                           1) Li, G and DirMtx=Hessian
                           2)           DirMtx=Hessian only
                           3) Likelihood only
                           4)           DirMtx=G'G only
                           5) G and Li only

                 li and grad will be (1x1) and (Kx1) respectively
                 DirMtx, if requested, will be KxK.
                 Missing is returned for components not requested; e.g.
                        for the DirMtx when arg1=3.

          !!!   MODEL is assumed to read the entire dataset  !!!
                MODEL returns + (true) loglik, etc.; and NOT -loglik !!

   parnames = Names corresponding to B0 (used to display constraints?)

GLOBAL INPUTS:

 _maxiter =  Maximum # of iterations
 __rows  = # of observations per read (0= program set automatically)
 _GRBLCVG = Global used to force convergence (set by user supplied routines)
            If set to 1, then force convergence.
            If set to 2, then force abort.
 _dsp_Frc = If user forced convergence, if _dsp_frc= 1, then output gradient
 _mlmnstp = Minimum step size in several places.
 _mlepsb,_mlepsg = Stopping criteria (b elasticity and grad elasticity)
 _mlatmax = How many stopping criteria to satisfy (1,2, or 3)
            Third criteria is dl.
 _mlvrylo = Alternate convergence criteria. If =1.0 or mlatmax=1, ignored;
            otherwise, if dl is  < _mlepsb/_mlvryol, and other terms are
           < _mlepsb*_mlvrylo, then converge.  In other words, if
            dl is getting tiny, the other's don't have to be as stringent.
            A typical value for _mlvrylo would be 500.

 _mlbhhh   = Controls use of BHHH as backup technique. 1=Yes, 0=No.
 _mlsteep  = Controls steepest descent as backup technique.
             If -0, don't do it.
             If 1, use all coefficients.
             If 2, then first use all coeffs, then try 1-at-a time
             (starting with largest absolute value) if no improvement
 _mlmarq   = Use Marquardt algorithim if direction matrix is not
             positive definite:: H = (H + t*I), t>0.
              0=No, 1=Yes
 _mldornd  = Controls use of random search as a backup technique.
             The value sets how many alternate á vectors to try. If 0,
             don't do it.
 _mldorn2  = Used with _mldornd.  How many successes to look for.
 _mlbkup   =  1=Try moving backwards if stuck (e.g.; at an inflection point)
 _mlrndvc  =  Controls H in N(á,H) selection of random á (in RANDOM
             algorithim).  1=Use Hessian.  2=Use BHHH.
                           Values < 1.0 mean use a diagonal matrix equal to
                            _mlrndvc * á.
 _mlelist  = List of techniques. These are done first.
             Set to 0 to suppress.
             Otherwise, use this to specify a fairly intricate search
             path.  _mlelist should be Nx5, with the 5 columns:
                  arg~step~freq~min~max
                        arg == see hessalg above
                        step == -1 = Backward
                                0  = One
                                1  = Golden
                        freq ==  every nth iteration
         (values below 0 have special, archaic meanings -- no longer supported -
          superseded by the use of min/max range

                           Less frequent methods should precede more
                           frequent methods (the first match is the one
                           that's done).

                        min  == First iteration to consider this on
                        max  == Last iteration to consider this on
  _no_cls     = 0 = clear screen between iterations
		1 = do NOT clear screen (useful if you write status messages to lower portion of screen)
  
  _useallb   = 0 = If first unbound direction yields no improvement, then do
                   next technique (e.g.; give up after only 1 unbound
                   direction).
               1 = Keep trying unbound directions before doing next technique
                   (e.g.; keep trying unbound directions until either none
                   left to try, or success).
               2 = Try all unbound directions for each technique, &
                   choose best.
                   (e.g.; choose the best unbound direction, only if none offer
                    an improvment then attempt next technique)>

 _error_out  = 0  Do not output error messages to output file
               1  Output error messages to output file
 
Basic idea:
      A "technique list" is constructed from _mlelist AND from hessalg,
dogolden and the _mlxxxx parameters.  Starting from the top of the _mlelist,
techniques are attempted.  If a technique fails to yield an improvement (and
gradient is non-zero), then the next technique (that is permitted on this
iteration) is attempted.  If none of the techniques on the _mlelist work (or if
_mlelist=0, in which case there are none!) then the strategies implied by the
set of hessalg, dogolden and _mlxxx parameters is attempted.

 By 'not work' I mean that either the inverse does not exist,
  or no improvment (regardless of stepsize) is possible.


Example.  For TOBIT procedure, b0=b_ols|variance.
{li,b,vc,vcused}=
       maxlikc(b0,                              @ starting values @
               miss(1,1),                       @ no constraints @
               1,                               @ use NR as initial algorithim@
               0,                               @ do not use GOLDEN @
               0,                          @ do not print intermediate results@
               " TOBIT MLE ",              @ displayed on screen @
              &TOBIT,                    @ name of TOBIT likelihood,etc. proc@
              parnames) ;          @list of parameter names @


------ COMPVC ------

        {vc,vcused}=compcvc(b0,vctype,f0,&model);

INPUTS:
 b0     =  A Kx1 vector of starting values.
 vctype = Algorithim to use for variance-covariance matrix:
            1= -inv(H)                  : hessian
            2= inv(G'G)                 : BHHH matrix
            3= inv(H)*(G'G)*inv(H) -- 'white' heteroscedastic consistent matrix
 f0 = Jx1, (J <= K) vector 'free variables'.  This is an index, into b, of
      variables that were maximized over.
      Eg: if K=4, and f0=1 3 4, then b[2] is treated as a constant.
      If f0=miss(1,1); then all variables are considered free to vary.
      FIGFREEV can be used to find f0, given a cnstlist.

 model     = Name of proc that will compute model.
             This proc should be of form:
               {li,grad,dirmtx}=model(b,narg)
                        b=beta
                        narg=
                           2) Compute Hessian
                           4) Compute G'G
                (see description under MAXLIKC)

OUTPUTS:
  vc    = A JxJ  (NOT KxK) covariance matrix, corresponding to the J 'free
          to vary' variables.
  vcused  = Which algorithim is actually used.  If 1 or 3 are requested, but
            are incomputable, then 2 is used.  If VC matrix can not
            be computed, an identity matrix is returned in VC, and vcused=4.


---- DO_TSTAT  --
 Compute t-stats and standard errors, given Beta, VC matrix, and "free index")
 The Free INdex points to non-constrained variables (if = ., then all
 variables are unconstrained).  Constrained variables will have their
 respective elements of stderr and tstat set to missing).

  call as : {stderr,tstat}=do_tstat(beta,vc,freeindx)
  Note: freeindx=0 means "no constrained elements"

------ FIGSTUFF ------
  Compute and output elasticites and E(Y) for each observation.
  Output either to default output file, or to user defined GAUSS dataset.

     Call as:
        call figstuff(&preds,beta,f1,xns,dataset)

  Inputs:
        Preds: a proc that computes elasticities and E(Y). See below.
        Beta: Current coefficient vector. Passed to Preds.
        F1 :  Current dataset. Assumed known by Preds.
        XNS: Names of variables involved in model. Assumed know by Preds.
        Dataset: Name of dataset associated with F1.

  Note that all output is to datasets or output files.

        Preds is a user written procedure (eg; FIGPRED) with the structure:

       {elasts,yhat,idvals}=preds(beta,idindx);

 Inputs:
                Beta: as defined above.
                Idindx : Index into variables in DATASET.
 Outputs:
                elasts: Elasticities for all variables in XNS.
                yhat: E(Y)=Yhat=Predicted value of Y
                idvals: Values of variable indexed by idindx. If idindx =0,
                  miss(1,1) is returned.

                Note that E(Y) and ALL the elasticities should be computed.

      FIGSTUFF will decide which of these values to output; eg. only 1 of 3
              elasticities.



**************/

declare _mlegrad = 0 ;         @ gradient at maximum @
declare _mlehess = 0 ;         @ hessian at maximum @
declare _mlebuse = 0 ;
declare _mliter = 0 ;
declare _tchlist = 0 ;      @ used to observe techniques used this run @
declare __rows = 0 ;         @ # of obs per read, 0=program sets @
declare matrix _grblcvg = 0 ;@ set by user supplied routine (1=cnvrg,2=abort)@
declare _mlepsb = 0.00001 ;
declare _mlepsg = 0.00001 ;
declare _mlvrylo = 1.0 ;    @ if dl very small convergence tester @
declare _mlatmax = 3 ;   @ most stringent (1 2 or 3 ) @
declare _mlbhhh = 1 ;   @ 1=use BHHH as backup techniqe, 0=no @
declare _mlsteep = 0 ;      @ 0=no, 1=all , 2=one at a time @
declare _mlmnstp =0.0001 ;
declare matrix _cowstep = 0 ;
declare matrix _maxiter = 50 ;        @ max iters, a global @
declare matrix _dsp_frc = 1 ;   @ output gradient if user forced convergence@
declare _mlmarq = 1;      @ 1=use marquardt algorithim if necessary @
declare _mldornd = 0 ;    @ 0=don't do it @
declare _mldorn2 = 5 ;    @ 5 successes @
declare _mlrndvc = 2 ;    @ use BHHH in random algorithim @
declare _mlbkup = 1 ;    @ 1=do backup if needed, 0=don't @
declare _mlelist = 0 ;
declare _useallb = 2 ;  @ 0 =don't try unbound direction @
                        @ 1=try all unbound directions before giving up@
                        @ 2=select best unbound direction @

declare _CtCall=0 ;    @ counter of number of times MODEL callled @

declare _no_cls=0 ;   @ 0=do CLS before writing Iter=1 status message at top of screen, 1= do NOT cls @
declare _error_out=1 ;  @ write error messages to otuptu file @

@ used to hold messages for SET in ALG_SET @
declare mlmess01=0;declare mlmess02=0;declare mlmess03=0;declare mlmess04=0;
declare mlmess05=0;declare mlmess06=0;
declare mlmess07=0;declare mlmess08=0;declare mlmess09=0;
declare mlmess10=0;declare mlmess11=0;declare mlmess12=0;declare mlmess13=0;
declare mlmess14=0;declare mlmess15=0;declare mlmess16=0; mlmess17=0;
/********************************/

proc (3)=maxlikc(b0,cnstlst0,hessalg,dogolden,prnt,header,&model,parlist);

  local model:proc ;
  local b,li ;
  local q0,q1,q2,g1,c1,c2,h1,dq1,tmp,tmp2,mname2;
  local relgrad,relb ;
  local forcec,iter,mname,kb,stp,prnt2,foo,mtype0,mm,methchng,ok,ar1,ac1,amess;
  local dvs1,dvs2,invh1,maxstep,gtmax,ordst,minstep,attempt1 ;
  local ifree,freevars,freevar0,freevar,freeindx,freeind0;
  local g0,db,TMPINVH1,stepname;            @ DFP and BFGS need these @
  local techlist,attempt,mtype,alglist,steptype,priorinv,do_1ata ;
  local bestli,savefree,savestp,savedvs,cnstlist,notok,prc3_ind ;
  local user_li,user_b;

    let mname=NR BHHH STEEPEST DFP BFGS OTHER RANDOM STEEPONE;
    let mname2=NR BHHH STEEPEST DFP BFGS USER RANDOM ;
    let alglist=NR BHHH Steepest DFP BFGS RANDOM ;
    let stepname=BACKWARD NOT_PD ONE GOLDEN ;

@ REMINDER:  Basic algorithim: B2= b1 - s*Inv(D)*G @
@ (Except for RANDOM ) @

   if hessalg==0 ;
     call grbl2_error(" ERROR: 0 is not an allowed matrix type in maxlikc ",1);
     stop;
  endif;

@ see if hessalg is string; if so, convert to number @
    tmp=0 $+ hessalg ;
    tmp2=indcv(tmp,mname2);
    if not(scalmiss(tmp2)) ; hessalg=tmp2 ; endif;

    ifree=miss(1,1);

@ initialize some stuff @

@ Find set of coefficients that are NOT equality contrained.
  And set equality constrained coefficients to their constraints,
  and remove equality constraints from list

  freevar0 is KBx1 -- is index of non-equality constrained variables @


@ see if hessalg is string; if so, convert to number @
   {b0,freevar0,cnstlist,notok}=figfree0(cnstlst0,b0);
   if not(scalmiss(notok)) ;
    if _error_out ne 0 ; output on ; endif;
     ? ;  call boldon ;
" Starting values are IMPERMISSIBLE.
  The following INEQUALITY constraints are violated:";
       call mle_show_constraints(cnstlist[notok,.],parlist);
       call normalon ;       output off ;
       call g_pause(10,-1); 
       _mlebuse=b0;
       retp(miss(1,1),b0,miss(1,1));
   endif ;

   freeindx=packr(miss(freevar0,0));

@ is there anything to do ? @
      if scalmiss(freeindx) ;
       if _error_out ne 0 ; output on ; endif;
         ? ; ?;
         " All parmeters are bound at start !  No estimation possible !  ";
          ? ; ?;
	  output off ; 
         _mlebuse=b0 ;
         retp(miss(1,1),b0,freeindx);
      endif;

@ yup, there is .... @
      _mlegrad=miss(1,1); _mlehess=_mlegrad ;
      priorinv=_Mlegrad;              @ none available @
      iter=0; kb=rows(b0);
      q0=0; q1=0;  q2=0; li=0;b=0;
      g0=0; g1=0;
      c1=b0; h1=0;
      db=0; ordst=0;
      user_li=miss(1,1); user_b=b0;

@ make an instruction list, of techniques, that tells which techniques to
  use on each ATTEMPT of an iteration @

     techlist=mk_techs(hessalg,dogolden,kb) ;
     _tchlist=techlist ;
     do_1ata=not(scalmiss(indnv(8,techlist[.,1])));

@ print some stuff at top of screen @
   if _no_cls==0 ;cls ; endif;                    @ clear screen .. @
   {ar1,ac1}=max_prc2(header,0,0,0,0,0,0,miss(1,1),0,0);  @ just print header @
   locate 4,1 ;

@ Each technique/attempt combo has different requirements and outputs @

   _Ctcall=zeros(3,1);   @ col1= # of iterations, col2=number of times modelproc called in THIS iteration @
			 

iter1: ;                    @ -------  start iteration here --------- @
   _Ctcall[3]=_ctcall[3]+1; @ iteration counter @
   _ctcall[2]=0;	   @calls within iteration counter @
   attempt=0 ;              @ start at bottom @
   attempt1=0 ;           @ count of actual attempts tried @
   scroll ar1|1|ar1+2|80|0|7 ;  @ clear some space @


   iter=iter+1;

   if iter>_maxiter ;                  @ too many iterationS @
         if _error_out ne 0 ; output on ; endif;
        " MAXIMUM # of iterations exceeded = " _maxiter ;
        " Current coefficient values will be used. ";
	  output off ;

          goto done0 ;
    endif;

     methchng=0;
     {forcec,prnt,mm}=mloptsgt(prnt,ar1);        @ allow for user options @
     if mm==1; methchng=1 ; endif;
     if forcec>0 ;
          goto forcecvg ;
     endif;



iter2: ;        @ ------------------- come here when an an attempt fails @


@ try first, or next, technique in techlist. Note that attempt is the
  "last attempted" technique, NOT the "next" (so, next =attempt+1)
@

    {attempt,mtype,steptype,minstep}=gt_tech(attempt,techlist,iter);

    attempt1=attempt1+1 ;

@ mtype is zero either specfically (i.e. -- as coded for in _mlelist)
  or if at end of attempt list @
    if mtype==0 ;               @ signals that there is no more @
       goto done0 ;             @ +grad but done @
     endif;

@ compute li,grad, and desired 'hess' for dataset @

@ WE either compute hessian (h1), or we directly compute inv(hessian)
  (Tmpinvh1).  Note that invh1 =inv(h1), is used below by NR and BHHH,
  while invh1=tmpinvh1 is used below by DFP and BFGS.  Note that all
  gradients or hessian'ettes returned by MODEL should be
  KBxKB -- for all vars, including constrained vars.  Constraints are
  taken care of in max_prc3, or at iter3 @


  if attempt1==1 ;     @ then we are using the default technique for this iter@
@ get loglik,etc . TMPINVH1 is used only by DFP and BFGS @

         {Q1,G1,H1,TMPINVH1}=
                  max_prc1(&model,mtype,c1,priorinv,iter,g0,db);
           forcec=_grblcvg ; _grblcvg=0;
           if forcec > 0 ; goto forcecvg; endif ;

@ Uncomputable likelihood ? @
         if scalmiss(q1) or scalmiss(g1);
          " q1 " q1 ; " g1 " g1 ;
           q1=miss(1,1) ;
          if _error_out ne 0 ; output on ; endif;
           format /lz 10,5 ;
   ?;?;        " Can't compute likelihood at these coefficients:"; c1' ;?;?;
           output off ;
          if _error_out ne 0 ;   call g_pause(5,-1); endif;
         forcec=1;
           goto done1 ;
        endif ;

        dq1=q0-q1;                        @ displayed as dL @
        mtype0=mtype ;                  @ needed if Random called later @

@ We need this if we ever do 1-ata-time steepest?@
        if do_1ata==1 ;
           ordst=rev(sortind(abs(g1))) ;
        endif;

  else    ;                @ ******** trying backup techniques **** @


    if mtype==1 and mtype0 ne 1;                    @ nr as backup @
         _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; 
        {foo,foo,h1}=model(c1,2);    @ just get hessian matrix @
       mtype0=1;
    elseif mtype==2 and mtype0 ne 2 ;                    @ bhhh as backup @
         _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; 
        {foo,foo,h1}=model(c1,4);    @ just get gg' matrix @
        mtype0=2;
    elseif mtype==4 or mtype==5 ;    @ get tmpinvh1, for DFP or BFGS @
       if mtype ne mtype0 ;
          tmpinvh1=do_tmph1(mtype,g0,g1,priorinv,db,kb);
          mtype0=mtype;
       endif ;
    elseif mtype==3 or mtype==8 ;     @ all, or 1-ata, steepest as backup @
         h1=eye(kb);
    elseif mtype==7 ;                   @ may need to set random VC mtx @
        foo=varget("_mlrndvc");
        if foo==1 and mtype0 ne 1;
           _CtCall[2]=_CtCall[2]+1 ;    _CtCall[1]=_CtCall[1]+1 ; 
           {foo,foo,h1}=model(c1,2);
           if not(scalmiss(h1));
               h1=-h1 ;                    @ used as VC of multi-var normal @
               mtype0=1 ;                           @ h1 is now from hess @
           endif;
        elseif foo==2 and mtype0 ne 2 ;
            _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; 
            {foo,foo,h1}=model(c1,4);     @ bhhh @
            mtype0=2 ;                   @ h1 now from g'g @
        endif;                           @ use fraction @
        if (scalmiss(h1) or foo lt 1) and (mtype0 ne 8) ;
            if (foo le 0) or (foo ge 1); foo=0.1; endif;
@note: cv=inv(h1), so use reciprocal of desired 'cv' @
            h1=eye(rows(c1))*(1/foo);   @ use fraction of beta @
            mtype0=8;
        endif;

        forcec=_grblcvg ; _grblcvg=0;
        if forcec > 0 ; goto forcecvg; endif ;

    endif ;                     @ mtype # if @

  endif;                        @ attempt1 # if @


@ display user info, etc. @
   if iter==1;foo=miss(1,1); else; foo=dq1 ; endif;

   {ar1,ac1}=max_prc2(header,prnt,iter,attempt1,
                       q1,foo,g1,mtype,steptype,ordst);


@ Note that binding constraints are nested within techniques ! @
@ furthermore, hessian and gradient elements are computed as
 if there are no binding constraints @

@ ... start with "equality only" constraints @
     freevar=freevar0 ; freevars=miss(1,1); bestli=1e25;
     ifree=0 ;                          @ signals 0th constraint set @

@ ---------- Hop here if a constraint(s) is binding  ----
  Note: random search never jumps here
 Freevar and Freeindx are KBx1 -- they are a list of candidate
  "free to vary "variables.  Freevar is either equal to freevar0, or
  was pulled off of the freevarS matrix
-------------------------------------------@

iter3:  ;
     freeindx=packr(miss(freevar,0));
     {forcec,prnt,mm}=mloptsgt(prnt,ar1);
     if mm==1; methchng=1 ; endif;
     if forcec>0 ;
          goto forcecvg ;
     endif;

@ NOte: for 1-ata-time, we can use standard constraint 'enforcement' @
@  for other random search, max_prc4 does it all @
    if mtype==8 ;               @ 1-ata-time @
       freeindx=ordst[steptype] ;               @ only 1 coeff free ata-time @
       if freevar0[freeindx] == 0 ;        @ an equality constrained var@
           goto iter2 ;                  @ so get next "technique" @
       else ;
          freevar=zeros(kb,1); freevar[freeindx]=freeindx ;
       endif;
       invh1=1 ;
       prc3_ind=freeindx  ;    @ Used to set priorinv @

    else;                       @ not 1-ata-time steepest @

@ Compute inverse, freevarXfreevar. possibly marquardt, possibly quasi-newton @
@                                   possibly backwards @

       invh1=max_prc3(h1,freeindx,tmpinvh1,mtype,steptype,kb);
       prc3_ind=freeindx  ;    @ Used to set priorinv @

@ Note: invh1 has dimension equal to the # of free to vary variables -- it is
        'defined' only for non-0 freevar variables @
    endif ;

    if scalerr(invh1) or ismiss(invh1);
        gosub unusab("inv(h1)") ;
        goto iter2 ;            @ and try next method @
     endif;

@ Either compute step length, or start random searches. If random searches,
  then call max_prc4 -- note that max_prc4 does all constraint enforecment @

     if mtype==7 ;                        @ RANDOM SEARCH @
        G1=0;                       @ cancel it @
        {c2,q2,ok,forcec,prnt,mm}=
           max_prc4(c1,invh1,-q1,&model,freevar0,
                    cnstlist,prnt,ar1,minstep);

        if forcec> 0 ; goto forcecvg ; endif;
        if mm==1; methchng=1 ; endif;

        forcec=_grblcvg ; _grblcvg=0;
        if forcec > 0 ; goto forcecvg; endif ;

        if ok==0 ;                                 @ no improvement @
            goto done0 ;

@ Note that random search requires improvements in (li and beta)--
  not elasticities thereof @
        else ;               @ got an improvement -- do next iteration? @
            q2=-q2 ;            @ make negative (minimization @
            foo=0;
            if (abs(c2-c1) le _mlepsb) ; foo=foo+1; endif;
            if (abs(q2-q1) le _mlepsg) ; foo=foo+1; endif;
            if foo> minc(2|_mlatmax);
                  goto done1 ;
            else      ;                 @ keep on trucking @
                stp=1; dvs2=c2-c1 ;relb=miss(1,1);relgrad=relb;  @ used below @
                goto iterend ;              @ do some end'o iteration stuff @
            endif;
        endif ;                         @ improvement @
    endif ;                             @ random method @

@ If here, did not do random ..... @
@ compute step direction. Notice use of - sign, and step only using
 unconstrained vars (invh1 incorporates these constraints, see max_prc3   @
     dvs1=-invh1*g1[freeindx] ;     @ direction vector, of unconstrained vars@

     dvs2=zeros(kb,1);
     dvs2[freeindx]=dvs1 ;         @ note: 0 step for constrained vars @

@ find maximum step size, this freevar and this direction @
@ (check only using non-equality cosntraints. Maxstep is 1x1 @
     maxstep=figmaxst(c1,dvs2,cnstlist,freevar,steptype==1,minstep) ;

@ for this set of freevar, a constraint is binding? @
     if maxstep le minstep ;

@ Get possible freevars, given current freevars and direction, etc.
  Note that we add to old set of possible freevars, and we do not add
  duplicates.
  Point: FreevarS is only non-missing if a constraint is binding !  @

           if mtype==4 or mtype==5 ;
             if _error_out ne 0 ; output on ; endif;
               ?; format /ld 3,0 ;
" ! Warning: DFP/BFGS used, but constraints are binding (@ " iter attempt ")";
                  ?; output off;
            else ;
                  prnt2=prnt; prnt=0;
                   gosub unusab("Free_var") ;
                  prnt=prnt2 ;
                  if prnt==1 ;
                      format /ld 4,0 ; 
                       output on ; 
                     " ... Unusable set of free_vars = " freeindx';
                     output off ;
                  endif;
            endif;

@ augment freevars, given current freevar and direction etc.
 Reminder: freevar is a list of variables that failed the constraint test!
Figfreev looks for subsets of freevar that pass the constraint test, at
current direction.  Freevars contains prior set of variables that pass
constraint test, but at different directions (associated with prior
realizations of freevar).  Thus, before using a vector of off freevars,
it needs to be tested against its OWN direction .
Note: only subsets of freevar that are not already on freevars will be added
to freevars (no point in duplicating)
@

          freevars=figfreev(cnstlist,c1,dvs2,steptype==1,
                                          minstep,freevar,freevars);

/*" after figfreev, freevars " freevars;call wait ;*/

          ifree = ifree+1;
@ get next freevar candidate @
          if (ifree le rows(freevars)) and (not(scalmiss(freevars)));
               freevar=freevars[ifree,.]' ;
               goto iter3 ;                    @ try with this freevar @
          else ;
              if prnt==1 ;   output on ; endif;
   ? ; " Can not find unbound set of parameters for this direction. "; ?;
               output off;
               goto iter2 ;            @ get next attempt @
          endif;

     endif ;

@ when here, got a candidate step that will NOT go out of bounds @
@ So .. take a step ! @
     {q2,stp,c2}=makestep(-q1,c1,dvs2,maxstep,steptype==1,&model,minstep);
     forcec=_grblcvg ; _grblcvg=0;
     if forcec > 0 ; goto forcecvg; endif ;
     if scalmiss(q2);   @ couldn't find an improvement.@
        gosub unusab("Lik_got") ;

        if _useallb gt 0 and ifree lt rows(freevars) ;
           ifree=ifree+1 ;
           freevar=freevars[ifree,.]' ;
           goto iter3 ;
        else ;                     @ give up if constrained, try new direction@
            if _useallb == 2 and bestli < 1e25 ;   @ got an earlier improve@
               goto gotimp ;
            else ;
                goto iter2  ;    @ try next method @
            endif;
        endif;
     else ;
        q2=-q2 ;          @ makestep returns + loglik @
     endif;

     if scalmiss(stp)  ;        @grad +, but no improvement... @
          gosub unusab("Step_dir") ;
          goto iter2 ;          @ try next method @
     endif;

@ an improvement, do we want to try other constraint sets @

gotimp: ;
     if _useallb==1 and not(scalmiss(freevars)) and prnt==1 ;
          format /ld 4,0 ;
           output on ;
         " Using free_vars: " freeindx';
          output off ;

     elseif _useallb==2 and not(scalmiss(freevars));

        if q2<bestli;                           @ minimizing -loglik @
            bestli=q2 ; savedvs=dvs2 ;
            savefree=freevar ;    savestp=stp ;
        endif;
        if ifree lt rows(freevars) ;
           if prnt==1 ; format /ld 4,0 ; 
                  output on ;
                 " Looking for better free_vars then: " freeindx';
                output off ;
           endif;
           ifree=ifree+1 ;
           freevar=freevars[ifree,.]';  @nb: freevar is KBx1 ! @
                           @ freevars is M*xKB (M* changes over time @
           goto iter3 ;
        else ;
           freevar=savefree; stp=savestp ;
           dvs2=savedvs ; q2=bestli ;
           freeindx=packr(miss(freevar,0));
           if prnt==1 ; format /ld 4,0 ; output on ;
                 " Using free_vars:"  freeindx';
                output off ;
           endif;
        endif;
     endif;

@ ARE WE MAXED OUT -- check 'elasticity' wrt coefficients and dL ? @

@ _mlatmax varies between 1 and 3:
   1 = Done if 1 critera. 2= .. if 2, 3= .. if 3
@
     foo=0;
     relgrad = abs(g1).*(maxc(abs(c1)|1)/maxc(abs(q1)|1));
     relb = abs(stp*dvs2)./maxc(abs(c1)|1);


     tmp=minc(_mlepsb|_mlepsg);
     if abs(dq1)<tmp ; foo=foo+1 ; endif;
     if abs(relb[freeindx]) < _mlepsb  ; foo=foo+1; endif ;
     if abs(relgrad[freeindx]) < _mlepsg ; foo=foo+1; endif ;

@format /ld 8,5 ;dq1 _mlepsb _mlepsg _mlatmax foo ; waitc;@

     if foo ge minc(3|_mlatmax) ; goto done1 ; endif ;

@ also check _mlvrylo condition @
     if _mlvrylo > 1 and _mlatmax>1 ;
        if abs(dq1)<(tmp/_mlvrylo) ;
          if abs(relb[freeindx]) < (_mlvrylo*_mlepsb)  ; goto done1; endif ;
          if abs(relgrad[freeindx]) < (_mlvrylo*_mlepsg) ; goto done1 ; endif ;
        endif;
    endif;

@ is the step size too small to continue to bother with (but +grad)?@
    if  (stp<minstep)  ;
         gosub unusab("TinyStep") ;
         goto iter2 ;                 @ try next method @
    endif;



@ Ready to take next iteration... @
iterend: ;                                      @ random jumps here @
     c1=c2 ;  q0=q1 ; q1=q2 ;
     db=stp*dvs2 ;                  @ db is used by DFP and BFGS@

@ set priorinv -- for possible use by DPF or BFGS @
     if mtype==7 or mtype==8 or mtype==3 ;
        priorinv=miss(1,1);             @ can't trust prior inverse(h1) @
     else ;
        priorinv=eye(kb);
        priorinv[prc3_ind,prc3_ind]=invh1 ;
     endif ;

@Note: DFP and BFGS are not recommended if/when model hits constraint.  The
       above setting of priorinv is just a heuristic solution @

@ print some current info @
     if prnt==1;
         scroll (ar1+3)|1|24|80|0|0;
         locate ar1+3,1 ;
          output on ; format /ld 3,0;
         ?; "Iter= " iter " , tech#= " attempt ;;format /lz 9,4 ;
          ", Step size= " stp ;; format /lz 12,6;
          ", LogLik= " q1 ;
          format /lz 10,5 ;
         " Coeffs   =" c1';
         " Direction=" dvs2';
         if mtype ne 7 ;
             " RelB     =" relb';
             " Gradient =" g1';
             tmp=zeros(rows(relgrad),1);
             tmp[freeindx]=relgrad[freeindx] ;
            "Relg w/con=" tmp';
         endif;

         output off;
      endif;

forcecvg: ;              @ if converge/ abort selected above@
      if forcec==1;
          if _error_out ne 0 ; output on ; endif;
           ?;?;     "  Convergence forced by user ! ";
           if _dsp_frc==1 ;
              " -Gradient: ";; format /lz 10,5 ; g1' ; ? ; ?;
           endif ;
           output off ;
           goto done1 ;
      elseif forcec==2 ;                     @ immediate abort @
          _mliter=iter ;
          gosub savestuf ;          @ save some globals @
         _mlebuse=c1;
          retp(miss(1,1),c1,freeindx);
      endif;


@ if change of method signaled @
       if methchng==1 ;
          $chrs(7);
           locate ar1+1,1 ;
           dogolden=g_yesno(" ? Use Golden ",dogolden);
           foo=g_askopt(" ? Algorithim ",alglist,hessalg,0,2,2);
           if not(scalmiss(foo));
              hessalg=indcv(foo,alglist);
           endif;
           if hessalg==6; hessalg=7; endif;
           techlist=mk_techs(hessalg,dogolden,kb) ;
           foo=csrlin ;
           scroll (ar1+1)|1|foo|80|0|0 ;
         endif;


@ save stuff that may be needed if CONVERGE hit @
   if not(badmtx(q1|c1)) ;
       user_li=q1 ; user_b=c1 ;
   endif;

@ Do next iteration @
   goto iter1 ;

/********************  JUmp Past here when done, for several reasons ...  **/

@ ending up with bad gradient @
done0: ;
          if _error_out ne 0 ; output on ; endif;
           ?;"  ! Warning: non-zero gradient at maximum !";
           format /lz 10,5;    "  -Gradient: " g1'; ?;
           output off ;

@ c'est tout, @
done1: ;

    foo=kb-rows(freeindx);

    if _error_out ne 0 ; output on ; endif;


@ is their a binding constraint that can't be beat @
    if (not scalmiss(ifree)) ;

@"Ifree " ifree';
"freevars " freevars;
" freeindx " freeindx' ;
" parlist " $parlist'; @

      if  ifree gt rows(freevars);
        ?;?; " Warning: ALL variables constrained at maximum! ";?;?;
        pause(0.5);
      elseif foo ne 0 ; format /ld 4,0 ;
        ?;?; " Warning: " foo " variables constrained at maximum: ";
           parlist[freeindx]=zeros(rows(freeindx),1);
           parlist=packr(miss(parlist,0));
            format 1,8;
           " : " $parlist';
            ?;?;
            pause(0.5);
      endif ;
     endif;
    output off ;

    gosub savestuf ;
    _mlebuse=b;
    retp(missrv(li,missrv(q1,0)),b,freeindx) ;


@ Subroutine to save some stuff ... @
savestuf: ;
           if forcec==1 ;
             li=user_li ; b=user_b ;
           else;
               b=c1 ;  li=q1 ;
           endif;

           if not(scalmiss(li));  li=-li ; endif;      @ ready to save it@
           _mliter=iter ;                           @ some background stuff @
          _mlegrad=-missrv(g1,0) ;
          _mlehess=-missrv(h1,0) ;     @ ... into permanent storage @
           _mlebuse=c1;

  return;

@ unusable hessian, etc @
unusab:;
  pop amess ;
       foo=csrlin ;
       locate ar1,ac1 ;
       call boldon ;
       if prnt==1; output on ; endif; format 1,8 ;
           $amess " unusable, "  $mname[mtype] ", "  $stepname[3+steptype]
           " Trying next method .. ";
        output off;
        pause(0.4); call normalon ;
        locate foo,1 ;
  return;

endp;
@ retp @
@model @
/**********************************/
@ create a technique list.  Two elemets:
  The hessian type -- 1,2,3,4,5,7,8   (NR BHHH STEEPEST DFP BFGS RANDOM 1-ata
  The steptype     -- -1,0,1        (Backstep, Stepsize of 1, GOLDEN)
  If 8 == then steptype refers is in 1..KB

Uses globals and hessarg,steparg @

  proc mk_techs(hessarg,steparg,kb);

 local tlist,foo1,foo2,amlelist,foo ;
  amlelist=_Mlelist ;

 if cols(amlelist)==5 ;                 @ use these instructions first @
    if hessarg==4 or hessarg==5 ;       @ use BHHH to start quasi-newtons @
        tlist=amlelist|(2~0~0~1~1)|(hessarg~steparg~1~1~_maxiter);
    else ;
        tlist=amlelist|(hessarg~steparg~1~1~_maxiter) ;
    endif;
 else ;                            @ else, use tlist as is @
    tlist=hessarg~steparg~1~1~_maxiter ;   @ first step is the one requested @
    if hessarg==4 or hessarg==5 ;       @ use BHHH to start quasi-newtons @
        tlist=(2~1~0~1~_maxiter)|tlist ;
    endif;
 endif;

@ bhhh backup@
 if (hessarg==1 or hessarg==4 or hessarg==5) and _mlbhhh==1;
     tlist=tlist|(2~steparg~1~1~_maxiter) ;
 endif;

@ is backstepping desired@
 if _mlbkup==1 ;

@ backup hessarg @
   if (hessarg==1 or hessarg==4 or hessarg==5 or hessarg==3);
      tlist=tlist|(hessarg~-2~1~1~_maxiter);    @ use -inv(h1) @
      if _mlmarq==1 and hessarg ne 3 ;   @ try non-marq augmented direction@
          tlist=tlist|(hessarg~-1~1~1~_maxiter);   @ use inv(h1), not invpd @
       endif;
    endif;
@ backup bhhh or bhhh itself @
   if  (hessarg==2) or
       ( (hessarg==1 or hessarg==4 or hessarg==5) and _mlbhhh==1) ;
         tlist=tlist|(2~-2~1~1~_maxiter);
         if _mlmarq==1 ;           @ try non-marq augmented direction@
             tlist=tlist|(2~-1~1~1~_maxiter);
          endif;
   endif;

endif;

@ next steps, bad NR,DFP,BFGS or BHHH -- is steepest descent desired @
        if (hessarg==1 or hessarg==2 or hessarg==4 or hessarg==5)
                 and _mlsteep ne 0 ;         @ either all or one at a time @
            tlist=tlist|(3~0~1~1~_maxiter) ;         @ give up on golden @
            if _mlbkup==1 ;              @ backstep here, too @
               tlist=tlist|(3~-2~1~1~_maxiter);
            endif;
        endif;

@ a desperation move: one-at-a-time (but only if RANDOM not chosen @
@ Note: one at a time is never user selected @
       if hessarg le 5 and _mlsteep==2 ;
            foo1=ones(kb,1) ;
            foo2=seqa(1,1,kb);
            tlist=tlist|((foo1*8)~foo2~foo1~foo1~(foo1.*_maxiter)) ;
       endif;

@ a desperation move, if 1,2,3,4, or 5 failed, then try random? @
         if hessarg le 5 and _mldornd>0 ;
           tlist=tlist|(7~0~1~1~_maxiter) ;    @note that steptype is ignored @
          endif;


@ tlist now contains instructions as to the order of techniques to be
  used.  Note that main program will need to behave in an optimal fashion,
  given this ordering @
     foo=zeros(rows(tlist),2);
     retp(tlist~foo);      @ add 2 columns for debug info @
 endp ;


/*****************************************/
@ pull info off of techlist @

 proc (4)=gt_tech(attempt,techlist,iter);
     local freq,mtype,steptype,iat,minstep,miniter,maxiter ;
     iat=attempt ;

foo1: ;           @ hop here if this technique not allowed on this iteration@
     iat=iat+1;
     if iat>rows(techlist)  ;  @ no more attempts available @
        retp(attempt,0,0,0) ;
     endif;

     mtype=techlist[iat,1];
     steptype=techlist[iat,2] ;      @ step type @
     freq=techlist[iat,3];
     miniter=techlist[iat,4];
     maxiter=techlist[iat,5];

@ note: since we may have a hard time starting, first time only techniques
        are allowed as backups @

     if iter>maxiter or iter<miniter ;   @ if not in permitted range @
        goto foo1 ;
     endif ;

     if iter>1 and freq==0 ;       @ first time techniques only first time @
        goto foo1 ;
     endif ;

     if freq< 0 ;                 @ iteration specific technique ? @
        if iter ne -freq ;
           goto foo1 ;
        endif;
     endif ;

     if freq>1 ;             @ is this every freq.th iteration technique @
       if (iter % freq) ne 0 ;    @ not time yet @
          goto foo1 ;
       endif;
     endif;

@ set up minimum step size to use (before giving up).  Note that for
 one-at-a-time, or steepest as backup technique, then use small steps @
      minstep=_mlmnstp;
      if mtype==8 or (attempt>1 and mtype==3 ) ;
         minstep=_mlmnstp^2;
         if steptype == 1 ;
            steptype=0 ;                   @and don't use golden @
         endif;
      endif;
_tchlist[iat,6]=iter ; @ some debug info @
_tchlist[iat,7]=iat ;
      retp(iat,mtype,steptype,minstep);
 endp;



/************************************/
@ invpd, with rescaling for nearly collinear matrices @
proc invpd2(mm);
    local ii,amult,am,ii10;
    if badmtx(mm) ; retp(miss(1,1)); endif;

    if _mlmarq ne 0;
        ii10=10;
    else;
        ii10=1;
    endif;

    ii=0;
    trap 1 ;
    do until ii ge ii10 ;          @ rescale for nearly collinear matrices @
           amult=1 + (ii*100) ;
           am=invpd(amult*mm);
           if scalerr(am) ;
              ii=ii+1;
            else;
              am=am*amult ;   @ inv(a*x)= inv(x)/a -- so a*inv(a*x)=inv(x) @
              ii=10;
            endif;
    endo;
   trap 0 ;
   retp(am);
endp;

/*************************/
@ compute g1,h1,etc. NOTE:: we always compute at least q1 and g1 @
@ (this routine is typically called as first attempt in an iteration @

  proc (4)= max_prc1(&model,hessalg,c1,preh1,iter,g0,db) ;

     local g1,h1,q1,tmpinvh1,foo,gnorm,kb ;
     local model:proc ;

    kb=rows(c1); tmpinvh1=preh1 ;
    if hessalg==1 ;                      @ NR @
            _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
           {q1,g1,h1}=model(c1,1);            @ compute grad and hessian @
           h1=-h1 ;
    elseif hessalg==2 ;                     @ BHHH @

           _CtCall[2]=_CtCall[2]+1 ;
           _CtCall[1]=_CtCall[1]+1 ; 


           {q1,g1,h1}=model(c1,0);
    elseif hessalg==3 or hessalg==8 ;                     @ Steepest @
            _CtCall[2]=_CtCall[2]+1 ; _CtCall[1]=_CtCall[1]+1 ; 
           {q1,g1,h1}=model(c1,5);
            if not(scalmiss(g1));
               gnorm=sqrt(g1'g1);
               h1=eye(kb)*gnorm;           @ will be inverted @
            endif;

@ for above, tmpinvh1 remains unset  @

    elseif hessalg eq 4 or hessalg eq 5 ;        @ DFP  or BFGS @
          _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
          {q1,g1,foo}=model(c1,5);
           tmpinvh1=do_tmph1(hessalg,g0,g1,preh1,db,kb);
           h1=0;                        @ we didnt compute hessian @

   elseif hessalg eq 7 ;
        q1=0;   foo=varget("_mlrndvc"); h1=miss(1,1);
        if foo==1 ;
           _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
           {q1,g1,h1}=model(c1,1);
           if not(scalmiss(h1));
               h1=-h1 ;                     @ used as VC of multi-var normal @
           endif;
        elseif foo==2 ;
            _CtCall[2]=_CtCall[2]+1 ;    _CtCall[1]=_CtCall[1]+1 ; 
            {q1,g1,h1}=model(c1,0);     @ bhhh @
        endif;
        if scalmiss(h1) or foo<1;                          @ use fraction @
            if (foo le 0) or (foo ge 1); foo=0.1; endif;
@note: cv=inv(h1), so use reciprocal of desired 'cv' @
            h1=eye(rows(c1))*(1/foo);   @ use fraction of beta @
            _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; 
            {q1,g1,foo}=model(c1,5);    @ likelihood needed, grad maybe later@
         endif;

   endif;

   if ismiss(q1) or ismiss(g1) ;
      retp(miss(1,1),miss(1,1),miss(1,1),miss(1,1));
   else;
      retp(-Q1,-G1,H1,TMPINVH1) ;    @ return -log,grad (minimizing) @
   endif;

   endp;


@-------------------@
@ compute DFP or BFGS inv(h1) @

 proc do_tmph1(hessalg,g0,g1,tmpinvh1,db,kb);

    local vv,shy,yshy,ys,shys,dg ;

@ if bad info, then just return identity matrix @
    if ismiss(g0) or ismiss(g1) or ismiss(tmpinvh1) or ismiss(db);
         retp(eye(kb));
    endif;

@ else, use 'em @
    dg=-(g1+g0) ;
    if hessalg==4 ;       @DFP from APtech and Judge et al pg 959@
       vv=tmpinvh1*dg ;
       tmpinvh1= tmpinvh1 + ((db*db')/(db'dg)) - ((vv*vv')/(dg'vv)) ;
    ELSEIF hessalg==5 ;               @ BFGS from APtech @
       shy=db-(tmpinvh1*dg); shys=shy*db';
       yshy=dg'shy; ys=dg'db;
       tmpinvh1=tmpinvh1 + ((shys + shys')   - (yshy*(db*db'))/ys)/ys;
    endif;
    retp(tmpinvh1);
 endp;

/************************/
@ display user info, etc. @

   proc (2)=max_prc2(header,prnt,iter,attempt,q1,dq1,g1,mtype,steptype,ordst);

    local mname,stepname,r1,c1,r2,c2 ;

    let mname=NR BHHH STEEPEST  DFP BFGS other RANDOM STEEPONE;
    let stepname=BACKWARD NOT_PD ONE GOLDEN ;

     c1=csrcol ;r1=csrlin ;

     output off ;
     locate 1,1 ;

     scroll 1|1|2|80|0|7 ;

     if not(scalmiss(header));    @  print header message in reverse video@
        call hiliteon;
        format 1,8 ; $header ;
        call normalon ;
     endif;

     format /ld 4,0; "Iter=" iter ;; format 3,0; ", " attempt ;;

     if not(scalmiss(mtype)) ;
      if mtype ne 8 ;
        format /lz 9,4; "-Loglik=" q1 ;;
        call boldon ;
        format /lz 6,2 ; ", dL=" dq1 ;;
        call normalon ;
        format 1,8;  " :" $mname[mtype] ", " $stepname[steptype+3];;
      else ;
       if not(ismiss(g1)) and not(ismiss(ordst)) and not(scalmiss(steptype)) ;
        format /ld 2,0;
          "  :: Steepest, grad[" ;; ordst[steptype] ;; "] =" ;;
           format /lz 9,4 ; -g1[ordst[steptype]]  ;;
       endif;
     endif;
    endif;
     "(A,C,M,P,?";

     c2=csrcol ;r2=csrlin ;

     locate r1,c1 ;

     retp(r2,c2);
  endp;



/*************/
@ return inv(h1), but only for currently free-to-vary variables @
@ Note: given non-free to vary variables are treated as constants, we should
  NOT include then in computation of direction matrix ! @

proc max_prc3(h1,freeidx,tmpinvh1,mdo,steptype,kb);

  local ok,invh1,mname,ordst,foo;
  ok=1;
  foo=missrv(freeidx,0);
  if mdo< 4 ;               @ check agains h1 @
     if (minc(foo) le 0) or (maxc(foo) gt rows(h1) ) ;
      if _error_out ne 0 ; output on ; endif;
        " % Inv(h1) ERROR, dirmtx type= ";; MDO ;
        " Free variable Index is " freeidx' ;
        " Size of h1 = " rows(h1);
        OUTPUT OFF;
         retp(miss(1,1));
     endif;
  elseif mdo==4 or mdo==5 ;
     if (minc(foo) le 0) or (maxc(foo) gt rows(tmpinvh1) ) ;
      if _error_out ne 0 ; output on ; endif;
        " % Inv(h1) ERROR, dirmtx type= ";; MDO ;
        " Free variable Index is " freeidx' ;
        " Size of h1 = " rows(tmpinvh1);
        OUTPUT OFF;
         retp(miss(1,1));
     endif;
  endif;

  if mdo eq 7;                        @ random @
         invh1=mlinvpd(h1[freeidx,freeidx]);
         if scalmiss(invh1);                    @if uninvertible @
            h1=eye(rows(freeidx))*(0.1);   @ use fraction of beta @
         endif;

  else;                         @ some gradient related step @
     if mdo < 4   ;                @ if bhhh or nr or steepest@
         if steptype ge 0 ;     @ not backwards @
@ better check that gradient matrix is correct size @
             invh1=mlinvpd(h1[freeidx,freeidx]);
         else ;                         @ back step request @
            trap 1 ;
            invh1=inv(h1[freeidx,freeidx]);
            if not(scalmiss(invh1)) ;
             if steptype eq -2 ;                  @ make it negative @
                invh1=-invh1 ;
             endif;
            endif;
            trap 0 ;
         endif ;

     else;                                    @ dfp,bfgs directly compute@

         invh1=tmpinvh1[freeidx,freeidx];      @ inv hessian @

    endif;

 endif ;                @ not random @
 retp(invh1);

 endp;


/********************/
@ do random search, using globals. Note that cv is the covariance matrix ONLY
  for the non-equality constrained variables.
  Note that   random search handles constraints
  internally (not by reiterating in main routine).
Stop trying after 5 sucesses or _mlrndo attempts @

  proc (6)=max_prc4(c1,invh1,q1,&model,freevar0,
                    cnstlist,prnt,ar1,minst);

     local model:proc ;
     local ii,ok,ido,brnd,bcand,qcand,q2,foo,tol,c2,cdiff,g1,relb,maxstep,dq ;
     local cc,mm,mm2,lin1,col1,ioutof,cvtmp,freeind,gotchol ;

      lin1=csrlin;col1=csrcol ;
      freeind=packr(miss(freevar0,0));       @ not equality constrained @

@      cvtmp=invh1[freeind,freeind];@
     cvtmp=invh1 ;

      cc=0;   mm2=0; gotchol=0;
      tol=varget("_mlepsb");
     ii=0;   ok=0; bcand=c1; qcand=q1 ; g1=0; relb=0;
     ido=maxc(varget("_mldornd")|75);
     do until ii==ido or ok ge maxc(varget("_mldorn2")|1);

        ioutof=0;
        {cc,prnt,mm}=mloptsgt(prnt,ar1);
        if mm> 0; mm2=1; endif;
        if cc>0 ;
           retp(bcand,qcand,ok,cc,prnt,mm2);
        endif;


        ii=ii+1;
        c2=c1 ;                          @ used to set equality constraints@

@ Note: if it were possible to evaluate at more than 1 beta, this would be
  quicker.  But, I assume that this is not done, so 1 at a time. @

@ get a random vector, mean of c1[freeind] and cv of cvtmp. ONLY for
  unconstrainted coeffs.  Note hackamatic time saving trick with cvtmp. @
doit1: ;                        @jump here if out of bounds beta selected @
       if gotchol==0 ;
          {brnd,cvtmp}=rndn2(c1[freeind],cvtmp,1)  ;
          gotchol=1;
       else;
          {brnd,cvtmp}=rndn2(c1[freeind],-cvtmp,1)  ;
       endif;
       c2[freeind]=brnd';           @ note brnd is a row vector @
       cdiff=c2-c1 ;

@ is it in bounds? @

     maxstep=figmaxst(c2,miss(1,1),cnstlist,freevar0,0,minst);

@ for this set of variables (freevar), is a constraint binding? @
     if maxstep le minst ;
         ioutof=ioutof+1 ;
         if ioutof < 10 ;        @ this direction is getting hopeless @
              goto doit1 ;
         else ;
            ii=ii+1;            @ so hopeless, we pretend we got 1 sucess@
            continue ;
         endif;
       endif;
       _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
       {q2,foo,foo}=model(c2,3);       @ get loglik @

       if not(scalmiss(q2));
       if q2 > q1  ;                      @ loglikelihood improved? @
          ok=ok+1;
          dq=q2-QCAND;            @ BETTER THAN PRIOR BEST GUESS @
          if dq > tol ;
             bcand=c2;
             qcand=q2 ;
             if prnt==1 ;
                locate lin1,col1;
                format /ld 6,0; " Random attempt " ii " = " ;;
                format /lz 9,5; qcand "  @ " bcand';
             endif;
          endif;
      endif;
      endif;

    endo;

@ if ok=0, bcand and qcand are == original values @
   retp(bcand,qcand,ok>0,0,prnt,mm2);

  endp;


/******************/

@ get inverse, using marquardt method if need be @
    proc mlinvpd(mtx);

     local ainv,i,ieye;
     i=1;

     if ismiss(Mtx); retp(mtx) ; endif ;      @ if missing @

     trap 1 ;
getinv: ainv=invpd(mtx);
     if scalmiss(ainv);
        if _mlmarq gt 0;
          mtx=mtx+eye(rows(mtx))*_mlmarq*(i^1.8) ;   @ speed up with each iter@
           i=i+1;
           if i>50 ;
               ainv=inv(mtx);
               trap 0;
               retp(ainv); endif;        @ forget it, or use non pd inv ! @
           goto getinv ;
        else;
               ainv=inv(mtx);
              trap 0;
              retp(ainv);         @ forget it, or use non pd inv ! @
        endif;
    else;
        trap 0;
        retp(ainv);
    endif;
endp;

/***********************/
@ return run time options. ALso, check globals _GRBLCVG
  where _GRBLCVG is set to 1 if convergence was forced
   and  _GRBLCVG is set to 2 if abort was forced.
@
     proc (3)=mloptsgt(prnt,ar1);
       local cc,ky,meth,ab,foo,tmp ;
       cc=0; meth=0; ab=0;

@ check _GRBLCVG @
        tmp=_grblcvg ;
        _grblcvg=0;
        if tmp==1;
          retp(1,prnt,meth);
        elseif tmp==2 ;
          retp(2,prnt,meth);
        endif;

aa:    ky=key ;
       if ky==0;
          if ab==1; cc=2; endif;
          retp(cc,prnt,meth);
       endif;

       if (ar1+9) gt 25 ; goto aa ; endif;   @ avoid trouble @

       if ky== ky==vals("H") or ky==vals("h") or  ky==vals("?");

          scroll (ar1+3)|1|(ar1+9)|80|0|7 ;
          locate ar1+3,1 ;
            "     Run Time Options: ";
            "          A = Abort estimation ";
            "          C = Force convergence ";
            "          M = Change algorithim and StepSize ";
            "          P = Output intermediate results ";
            "          D = Temporary exit to DOS (estimation resumes on EXIT) ";
            " Hit ESC now to select an option .. ";
            if g_pause(1.5,27);
                foo=g_askopt(" Select ","A"|"C"|"M"|"P"|"D"|"continue",
                            "continue",0,2,2);
                if scalmiss(foo) ; ky=0 ;
                  elseif foo$=="continue"; ky=0;
                  else ; ky=vals(""$+foo);
                endif;
             endif;
             scroll (ar1+3)|1|(ar1+9)|80|0|7 ;


        endif;

       if ky==vals("C") or ky==vals("c");
          cc=1 ;
       endif;
       if ky==vals("a") or ky==vals("A");
          ab=1 ;
       endif;
       if ky==vals("P") or ky==vals("p");
          prnt=1-prnt ;
       endif;
       if ky==vals("M") or ky==vals("m");
          meth=1 ;
       endif;
       if ky==vals("D") or ky==vals("d");
           foo=coreleft ;
           if foo<100000 ; format /ld 9,0;
             " Sorry, not enough free memory to open a DOS session ( " foo ;
              call g_pause(4,0);
           else;
              " Type EXIT to restart estimation .... ";
              dos ;
              cls ; " Restarting estimation ... ";;
            endif;
       endif;

       goto aa    ;

  endp;

/****************************************/
@  Given:
    cnstlist -- list of constraints
    b0  -- Current parameter values

   Set
       Equality constrained coefficients, list of free-to-vary
       coefficients, and list of non-equality constraints.

constraints:     1= lt   (le is converted to this )
                 2= eq

Retp(curvals_reset,freevar_nonequality,cnstlist_nonequality)

@

   proc (4)=figfree0(cnstlist,curvals);

   local nk,m,afreevar,amtx,cnstvals,amiss,im,ith,notok,st,minstep;

   nk=rows(curvals);
   m=rows(cnstlist);

   afreevar=seqa(1,1,nk);        @ assume all are free.. @

   

   if scalmiss(cnstlist) or missrv(cnstlist,0)==0;       @ no constraints ! @
      retp(curvals,afreevar,cnstlist,miss(1,1));
   endif;

@ Else, let's check. First, pull some info out of 'storage' @
   amtx=cnstlist[.,1:nk];
   cnstvals=cnstlist[.,nk+2];
   amiss=miss(1,1);

@ remove equality constraints.  We assume that syntax is correct (only
  one coefficient per each equality constraint @

  im=1;
  do until im> m ;
    if cnstlist[im,nk+1]==2 ;            @ an equality constraint @
       ith=maxindc(amtx[im,.]');     @ the single non-zero value @
       afreevar[ith]=0;              @ it aint' free babe @
       curvals[ith]=cnstvals[im] ;       @ enforce it @
       cnstlist[im,1]=amiss;
    endif;
    im=im+1;
  endo;

   cnstlist=packr(cnstLIST);      @ remove equality constraints @


@ are the starting values (adjusted for equality constraints) permissible ? @

    notok=miss(1,1); minstep=0.001 ;

    if not(scalmiss(cnstlist));
       im=1;
       amtx=cnstlist[.,1:nk];
       do until im>rows(cnstlist);
           st=binmxst(amtx[im,.],curvals,miss(1,1),cnstlist[im,nk+2],
                      minstep,1.0,afreevar);
           if st le minstep ;
               notok=appnd_it(notok,im) ;
           endif;
           im=im+1;
       endo;
    endif ;

    retp(curvals,afreevar,cnstlist,notok);

  endp;



/****************************************/
@  Given:
    cnstlist -- list of non-equality constraints
    wasfree  -- A candidate list of non-equality constrained variables
               that FAILED.
    oldfree  -- a matrix containing several lists of untried
                candidate "unconstrained" parameters.
                Unconstrained means "max step size in tested
                direction, using these parameters (setting other parameters
                at current values) is greater than 0 "
    current parameter values and current direction ---

   Find
       List(s) of new candidate "unconstrained" parameters (and add to
       oldfree).

Idea:
  A binding constraint may be due to any combinaton of the several
  variables involved in a constraint.  Thus, there may be several subsets
  of these variables that are unconstrained.  Try all of them!

 Note: 2  binding constraints may generate similar lists, so duplicates are
       removed.

straints:     1= lt   (le is converted to this )
              2= eq


Retp(oldfree|new_freevars)

Note that oldfree|new_freevars is MxK -- K=# of parameters, M=# of
candidate sets of variables that can be unconstrained.
Typically, this will be a 1xK containing no zeros (no binding constraints)


@

   proc figfreev(cnstlist,curvals,direc0,stype,minstep,wasfree,oldfree);

   local nk,amtx1,avals,m,amiss,im,ith,tmp,ii,i ;
   local badvars,st,maxsts,maxst0,direc;
   local amtx,cnstvals,maxst,this1,mm,imm,freevar,wasindx ;

    AMISS=MISS(1,1);
    if stype==1 ;             @ golden @
      maxst=2^10 ;
    else;
      maxst=1;
    endif;
    maxst0=maxst ;

@ only  1 free var, then no possible partial steps  @
   if sumc(wasfree./=0) le 1 ;
      freevar=miss(1,1);
      goto returnit ;
   endif ;

   nk=rows(curvals);
   m=rows(cnstlist);

   if scalmiss(cnstlist);       @ no constraints ! @
      retp(miss(1,1));
   endif;

@ pull some info out of 'storage' @
   amtx=cnstlist[.,1:nk];
   cnstvals=cnstlist[.,nk+2];

   if (maxc(cnstlist[.,nk+1]) ne 1 ) or (minc(cnstlist[.,nk+1]) ne 1 ) ;
      if _error_out ne 0 ; output on ; endif;
         call grbl2_error(" % Equality constraint passed to FIGFREEV ",1);
         output off ;
          stop ;
   endif;

  wasindx=packr(miss(wasfree,0)) ;

@ Any variable not included in WASFREE is assumed to be bounded, either
  as an equality constraint or as the result of a previously
  binding non-equality constraint. So set direc of non-WASFREE parameters to
  zero @

  direc=zeros(nk,1);
  direc[wasindx]=direc0[wasindx];

@ Given WASFREE, examine each constraint and find the maximum step we can take
  Presumably, one of these maxsteps will be 0 (the constraining constraints)

Note: freevar is being built -- it is the list of subset(s) of wasfree that
       are potentially unconstrained (at current direction).  Note that these
       will have to be tested against their own direction.
@

 badvars=amiss ; freevar=amiss ;

@ check wasfree (current set of variables) against M constraints @

  im=1;
  do until im>m;
     amtx1=amtx[im,.];
     st=binmxst(amtx[im,.],curvals,direc,cnstvals[im],minstep,maxst,wasfree);

@Note: maxst will shrink as prior constraints impose their maximum steps!@
@ Each row of Badvars contains the parameters in a binding constraint @

     if st le minstep  ;                   @ constrained in this direction ! @
        tmp=seqa(1,1,nk)'.*(amtx1./=0) ;
        badvars=appnd_it(badvars,tmp);
     else ;
        maxst=minc(maxst|st);        @ minimum of permitted steps @
     endif;
     im=im+1;
  endo;      @ get next constraint @


@ Each row of badvars refers to a binding constraints, and it lists
  the parameters of this constraint that are active (non-zero)
  Take badvars and create a list of freevar
Note: the base list is wasfree == any badvar not in wasfree MUST be
  excluded from freevar list
@

   if scalmiss(badvars) ;       @ no binding non-equality constraints @
        retp(oldfree);
   endif;

@ !!!! oh yeah, there are binding NON-EQUALITY constraints @

   maxst=0 ;              @ signal that current wasfree is bound @

@ however, there may be subsets of this wasfree that are not bound, given
  current direction vector ...@


@ Remove candidate constraints in which all variables are definitely bounded
  For those that contain variables not definitely bounded, create all
  possible freevar and add to freevars@

   m=rows(badvars);       @ m now refers to binding non-equality constraints@
   im=0;
   do until im ge m ;

      im=im+1;
      this1=badvars[im,.]';
      this1=packr(miss(this1,0));

@ This1 is a list of coefficients that appear in a binding constraint @

      mm=rows(this1);  imm=1;

@ check which vars in this are 'already' constrained @
      do until imm>mm ;
          if wasfree[this1[imm]]==0;
              this1[imm]=0;
          endif;
          imm=imm+1;
       endo;

@ see if all are fixed.  If so, this constraint has no impact, so don't
  bother with it @

       this1=packr(miss(this1,0)); @ removes "elsewhere constrained" vars @
       if scalmiss(this1);           @ all variables already fixed @
          continue ;                  @ so skip it @
       endif;



@ if here, one or more remaining variables to be fixed.  We do a combinatronic
  removal of variables that are not fixed (the variables identified in this1)
  and add these to the freevar matrix.  Note that wasfree
  contains info on the variables that are definitely fixed (e.g.; equality
  constraints, or single valued binding inequality constraints).
 Note that wasfree is 1xKB, while freevar is NNxKB (NN varies )
@

       freevar=dofreev(wasfree,this1,freevar);

   endo ;                         @ next constraint @

/*
format /m1 ;" after all calls to dofreev, freevar " freevar ; call wait ;
*/

@ we now have a new set of candidate free-to-vary lists, stored in
 freevar.  Go through them and find
  the max step size associated with each one.  If < minstep, remove it.
  THus, we will have only those freevars that one can actually step in
  the direction of.  The point is to remove constraints that are
  impossible, say due to one of the elements dominating a constraint @

  maxsts=zeros(rows(freevar),1);
  m=rows(amtx);
  ii=1 ;
  do until ii>rows(freevar);    @ for each candidate freevar@
    im=1;
    tmp=freevar[ii,.]';
    st=maxst0 ;        @ prior non-binding constraint may be indirectly killed@

@ check against each constraint that is a non-equality @
    do until im>m or st le minstep ;
@ the maxstep for this freevar is ... @
       st=binmxst(amtx[im,.],curvals,direc,cnstvals[im],minstep,st,tmp);
       im=im+1;
/*" in loop st,im, tmp' " st im tmp' ;*/
    endo;

    if st le minstep  ;                 @ can't use this candidate freevar @
       freevar[ii,1]=amiss ; maxsts[ii]=amiss;
    else ;
       maxsts[ii]=st ;
    endif;
    ii=ii+1;
  endo;            @ get next freevar @


  freevar=packr(freevar); maxsts=packr(maxsts);

/*" after removal of st=0s, freevar"; freevar; call wait ;*/

@ sort in order of maxsts @
  tmp=rev(sortind(maxsts));
  freevar=freevar[tmp,.];

@ jump here when ready to return freevar.. but first remove duplicatses @
returnit: ;

@ add freevar to oldfree ..............................................@
  if scalmiss(freevar);
    if scalmiss(oldfree);
      retp(miss(1,1));
    else ;
       retp(Oldfree);
    endif;
   endif;

@ remove duplicates............................. @
   if scalmiss(oldfree);
     oldfree=freevar[1,.];
     i=2;
   else ;
     i=1;
   endif;

@ compare new freevars to old list, add if new @
  do until i>rows(freevar);
     tmp=oldfree.==freevar[i,.] ;              @ element by element compare @
     tmp=sumc(tmp');
     if tmp<cols(oldfree);       @ there are no exact matches @
         oldfree=oldfree|freevar[i,.];
     endif ;
     i=i+1;
  endo;

  retp(oldfree);


endp;



/****************************************/
@ create freevar list from @
@ Freevar0 -- list of parameters that are not equality constrained,
              and are not subject to a binding, single variable,
              inequality constraint .
 this1     -- list of parameters that appear in a binding constraint.
               We selectively remove from freevar0 variables that
               appear in this1
 oldlist   --  A storage matrix.  After using this1 to selectively remove,
               place results at end of oldlist
@
  proc dofreev(freevar0,this1,oldlist);

      local nn,todo,tmp,i,foo,a,freevar,tmp2;
     nn=rows(this1);
     todo=(2^nn)-1;
     tmp=zeros(todo,nn);


     freevar=miss(1,1);

/*
"dofreev: freevar0'" freevar0';
" this1' " this1';
" oldlist " oldlist ;
call wait ;
*/


@ use binary expansion of nn.   @
     i=1 ;
     do until i>todo;
        foo=do_base2(i,nn);
        tmp[i,.]=(foo.*this1)';
        i=i+1;
     endo;
     foo=tmp./=0 ;
     a=sumc(foo');
     foo=sortc(tmp~a,nn+1); foo=foo[.,1:nn];

@ foo lists all possible variable combinations to constrain.
  Create freevar by zeroing out these combinations, from freevar0, one
  at a time @
  i=1;
  do until i>todo ;
     tmp=packr(missrv(foo[i,.],0));
     tmp=packr(indnv(tmp',freevar0));
     tmp2=freevar0;
     tmp2[tmp]=zeros(rows(tmp),1);
     freevar=appnd_It(freevar,tmp2);
     i=i+1;
  endo ;

@ we now have new list. Add to new ones to oldlist @
  if scalmiss(oldlist) ;                 @ nothing in old list @
     retp(freevar) ;
  endif;

@ else, stuff to check @
  i=1 ;
  do until i>todo ;
     tmp=oldlist.==freevar[i,.] ;              @ element by element compare @
     tmp=sumc(tmp');
     if not(tmp<cols(oldlist));       @ there is at least one exact match @
         freevar[i,1]=miss(1,1);
     endif ;
     i=i+1;
  endo;
  retp(packr(oldlist|freevar));

endp;

/**************/
@ convert number to base2, return in matrix @
  proc do_base2(i,ni);
   local foo,a,num;
   foo=zeros(32,1);
   num=trunc(i);

  i=0 ;
  do until num==0;
    i=i+1;
    a=num%2 ;
    if a==0 ; foo[i]=0; else ; foo[i]=1; endif;
    num=trunc(num/2);
 endo;
 retp(foo[1:ni]);
endp;

/***************************/
@ given freevar list, and constraint list, return maximum allowable step@
     proc figmaxst(curvals,direc,cnstlist,freevar,stept,minstep);
     local im,maxstep,st,kb,amtx ;


     if stept==1 ; maxstep=1024000 ; else ; maxstep=1; endif;

     if scalmiss(cnstlist); retp(maxstep); endif ;

     kb=cols(cnstlist)-2 ;
     im=1;
     amtx=cnstlist[.,1:kb];
     st=maxstep ;
     do until im>rows(cnstlist);
      if cnstlist[im,kb+1] ne 1 ;       @ an equality constraint? @
         continue ;                     @ then ignore it here @
      endif ;
      st=binmxst(amtx[im,.],curvals,direc,cnstlist[im,kb+2],
                minstep,st,freevar);
      if st le minstep ;
         retp(0);
      endif;
      im=im+1;
     endo;
     retp(st);
    endp;

/********************************/
@ return maximum allowable step size @
@ Note: recompute freeindx here as a favor to figfreev @

@ Note: E*Be + V*Bv < C
   E =Equality constrained vars
   V =Other vars
   Be, Bv - coeffieicnts
   C = Constraints

Thus,
    C2 = C - E*Be
is the portion of the constraint that is relevant to our free-to-vary
variables (the V variables).
@


proc binmxst(amtx0,baseval0,direc0,cnstval0,minstep,maxstep,freevar);

    local oldfrac,frac,maxf,minf,diff,oldmax,tmp,amtx1,freeind ;
    local val0,val1,basevals,direc,cnstval,tmp1;
@ used when checking for validity of "starting" Or random values @
    if scalmiss(direc0);
       direc0=zeros(rows(baseval0),1);
    endif ;

@ look for maximum step using ONLY non-constrained variables@
    freeind=packr(miss(freevar,0));

    if scalmiss(freeind) ; retp(0) ; endif ;

    amtx1=amtx0[freeind] ;
    direc=direc0[freeind] ;
    basevals=baseval0[freeind] ;

@ we want to compare to free variables in constraint @
@ val1-val0 is attributable to constrained vars @
    val0=amtx0*baseval0 ;      @ actual value of to compare to @
    val1=amtx1*basevals ;      @ non-fixed portion of this compare value @
    cnstval=cnstval0-(val0-val1);    @non-fixed portion of constraint@
    frac=minc(maxstep|1) ; oldmax=frac ;
    oldfrac=frac;

@ is max <  minstep ?  @
    do until frac < minstep  ;
        tmp1=(amtx1*(basevals+(frac*direc))) ;
        if tmp1 gt cnstval ;           @ le condition @
           oldfrac=frac ;
           frac=frac*0.5 ;
        else;
            break ;
        endif;
    endo;


    if frac<minstep;            @ binding!@

       retp(0);
    elseif frac ge maxstep ;       @ ok at current binding step @
       retp(maxstep) ;
    endif;

@ is it > 1 @
     if frac==oldmax ;              @ perhaps we can have a big step @
      oldfrac=oldmax ;
      frac=frac*2;
      do until frac ge maxstep ;
        tmp=amtx1*(basevals+(frac*direc))  ;
        if tmp < cnstval ;
           oldfrac=frac ;
           frac=frac*2 ;
        else;
            break ;
        endif;
     endo;
   endif;


@ when here, we have bounds (oldfrac and frac @
  minf=minc(frac|oldfrac);
  maxf=minc(maxc(frac|oldfrac)|maxstep);
  diff=maxf-minf ;
  do until diff<minstep;
    frac=minf+diff/2 ;
    if (amtx1*(basevals+(frac*direc))) gt cnstval ;  @ le condition @
           maxf=minf+diff/2 ;
    else;
           minf=minf+diff/2 ;
    endif;
    diff=maxf-minf ;
  endo;

  retp(minf);

endp;




/****************************************/
@ Compute stepsize, likelihood and coeff vector using this step size @
@ Return +loglikelihood,stepsize, coefficients @
@ Assume we are maximizing +loglik @

  proc (3)=makestep(q1,x,di,maxstep,dogolden,&model,minstep);

     local stp,c1,c2,q2,maxst,ok,foo,FORCEC ;
     local model:proc ;

      c1=x ;

strt: ;
     stp=minc(1.0|maxstep) ;
     _cowstep=stp;

     if dogolden==1 ;              @  use golden algorithim @
        stp=minc(1.0|maxstep) ;
        _cowstep=stp;
        maxst=minc(maxstep);
        {stp,ok,q2}=goldenc(x,di,stp,maxst,&model);
        forcec=_grblcvg ; @ STOP OR ABORT ? @
        if forcec > 0 ;
            retp(q1,miss(1,1),c1);
        ENDIF;
        if ok==0 or q1 > q2;           @ golden did not work @
           dogolden=0;
           goto strt ;
        else ;
           c2=x+(stp*di);

           retp(q2,stp,c2);
        endif;

@ stepsize of 1 (if possible ) @
      else;                             @ use 1.0, if possible ..@


stp1: ;
        c2=x+(stp*di) ;                @ try a step @
       _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 

        {q2,ok,foo}=model(c2,3) ;       @ and see if its an improvement @
        forcec=_grblcvg ;
        if forcec > 0 ;   @ STOP OR ABORT ? @
            retp(q1,miss(1,1),c1);
        ENDIF;

        if not(scalmiss(q2));          @ some kind of bad value @
          if q2> q1 ;                     @ yes @

             retp(q2,stp,c2);
          endif;
       endif;

@ if not, try a smaller step ? @

format /ld 9,4;

        stp=stp/5.0;
        if stp < minstep ;
@          "  Warning: Gradient positive at max. " ;@
            retp(q1,miss(1,1),c1);           @ signal + gradient at max @
        else;
            goto stp1 ;
        endif;

    endif;

  endp ;

/*********************************/
@ golden algorithim to maximize function.
        {stp,ok,q2}=goldenc(x,di,stp,maxst,&model);

            {stepsize,ok,val}=goldenc(b0,di,t,maxst,model);

         INPUTS
           b0   -- starting b
           di  -- direction vector
           t  --  starting stepsize
           maxst -- maximum stepsize permissible
           MODEL  - function to compute likelihood of dataset at b:
                        {li,arg2,arg3}=model(b,3)
                        (for arg2 and arg3, see above. They are not used
                           in this proc).
         OUTPUTS
            stepsize = step size
            ok   = Flag, 1 if good, otherwise 0
           val  = log likelihood at stepsize*di

                Eg: {stp,ok,val}

Note: proc badmtx detects missings, infs, and nans.  If these are found,
      something wierd has happened in the calculation of the likelihood function

@

proc (3)=goldenc(b0,di,t,maxst,&model);

local model:proc ;
local maxtries,  v1, v2, cd, cs, vr0, vr1, vr2, a0, b00, g0, g1, g2,
l, i, rd, rs, del, sig, s;
local ktht,diff0,diff1,eps,foo;


/* --------------- initializations ------------------- */
   maxtries=100;    @ Will stop trying to find min after this many tries @
   v1=(3-sqrt(5))/2; v2=(sqrt(5)-1)/2;
   cd=1; cs=1;    @ if 1, compute rd, rs, respectively, in seven below @
   eps=0.001;     @min step size tolerated @

/* --------------- test to see if [0,t] brackets minimum ---------------- */
step1: ;
        _cowstep=0;
     _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
    {vr0,foo,foo}=model(b0,3); vr0=-vr0;
    if badmtx(Vr0) OR (_GRBLCVG NE 0);        @ shouldn't happen, but .. @
        retp(0,0,0);
    endif;
step2: ;
        _cowstep=t;
       _CtCall[2]=_CtCall[2]+1 ; _CtCall[1]=_CtCall[1]+1 ; 
    {vr1,foo,foo}=model(b0+(t*di),3); vr1=-vr1;
    if badmtx(vr1) OR (_GRBLCVG NE 0);           @try smaller initial step @
       t=t/20;
       if t< eps ;
          retp(0,0,0);
       else ;
          goto step2;
       endif;
    endif;

if vr0 <= vr1; a0=0; b00=t; diff0=abs(vr1-vr0); goto seven; endif;

g0=0; g1=t; ktht=1;

/* -- if not, expand interval by increments of t until min is bracketed -- */
four:  g2=g1+(t*(1.5^ktht));     @bigger expansions at each iteration @
    ktht=ktht+1;

    if g2>maxst ;                    @ maxed out at maxstep @
        retp(g1,1,-vr1);
    endif;

     _cowstep=g2;
     _CtCall[2]=_CtCall[2]+1 ; _CtCall[1]=_CtCall[1]+1 ; 

    {vr2,foo,foo}=model(b0+(di*g2),3); vr2=-vr2;

    if badmtx(vr2)  OR (_GRBLCVG NE 0);  @ cant bracket minimum, so return VR1@ 
          retp(g1,1,-vr1);
    endif ;

    if vr2 >= vr1; a0=g0; b00=g2; diff0=abs(vr2-vr0);
        goto seven;
    else; vr0=vr1 ; vr1=vr2; g0=g1; g1=g2; goto four;
    endif;

/* ---------------- search to isolate minimum ---------------------------- */
seven: i=1;

eps=abs(b00-a0)/12;
diff1=diff0;                    @hack to skip first condition @
if eps<1e-05; eps=1e-05; endif;         @ just in case of wierdness@

do until i > maxtries;
    IF (_GRBLCVG NE 0) ; BREAK ; ENDIF;
    l=b00-a0;
    if abs(l) < eps; goto fin; endif;      @step size small? enuf @
    if diff1< (diff0/50) ;  goto fin ;endif; @ improvement small? enuf @
    if cd; del=a0+v1*l;
        _cowstep=del;
       _CtCall[2]=_CtCall[2]+1 ;  _CtCall[1]=_CtCall[1]+1 ; 
       {RD,foo,foo}=model(b0+(DI*DEL),3); RD=-RD;
    endif;
    if cs; sig=a0+v2*l;
        _cowstep=sig;
       _CtCall[2]=_CtCall[2]+1 ;    _CtCall[1]=_CtCall[1]+1 ; 
       {RS,foo,foo}=model(b0+(DI*SIG),3); RS=-RS;
    endif;
    if badmtx(rd|rs);            @if bad values, give up via a hack @
        b00=a0;
        goto fin;
    endif;

    diff1=abs(rd-rs) ;                          @dmh addition @
    if rd < rs; b00=sig; rs=rd; sig=del; cs=0; cd=1;
          else; a0=del; rd=rs; del=sig; cs=1; cd=0;
    endif;
i=i+1;
endo;

fin:
s=(a0+b00)/2;
 _cowstep=s;
 _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; 
{RS,foo,foo}=model(b0+(DI*S),3);


retp(s,1,rs);

endp;


/*******************/
@
  Generate predictions/elasticies.  Requires a routine that will
  predict E(Y;X,á) and dY/dX * x/y (the elasticity).
  This routine will only output chosen elastcities.

   Call as: call figstuff(&preds,beta,f1,xns,dataset)

      Beta is the current coefficients.
      F1 is the file handle.
      XNS is a Kx1 list of variables used in the model.  Note that beta
         need not have same dimension as XNS.
      Dataset = name of dataset assigned to f1 .

      Preds should be a name of a routine, say figpred, and is called
       {elasts,yhat,idvals}=preds(beta,idindx);

             idindx= Jth var, in file dataset, to be extracted.

           elasts= TxK matrix of elasticites, corresponding to XNS.
           yhat =Tx1 matrix of E(Y)
           idvals = Tx1, of var indexed by idindx

NOTE FIGSTUFF2 VARIANT: Permits output of more then 1 predicted
value.  YHAT2 should be a character variable, it is passed for info
purposes only (if yhat output selected, then both yhat and
yhat2 will should be generated and will be output).

So, figstuff is a subset of figstuf2

@


 proc figstuff(&preds,beta,f1,xns,dataset);
  local preds:proc ;
  retp(figstuf2(&preds,beta,f1,xns,dataset,miss(1,1)));
endp;

  proc figstuf2(&preds,beta,f1,xns,dataset,yhat2);

   local preds:proc ;
   local foo,iout,oname,iyhat,curlist,idindx,idname,elids,ovars,outf ;
    local elasts,dtaout,yhat,idvals ;
   local ir,ic,allobs,idid,aa ;
 local yesno ; let yesno=NO YES ;

  iyhat=0;  idindx=0; idname="0"; curlist=0;

get1: aa=g_askopt(" Select output-file type ","GAUSS"|"OUTPUT"|"HELP","HELP",0,2,2);
     if scalmiss(aa);         retp(0);     endif;
     if aa$=="HELP";?;
        " Select output-file type for predictions:";
        "   OUTPUT) To current output file     ";
        "    GAUSS) To a GAUSS dataset ";
        goto get1 ;
    elseif  aa$=="OUTPUT";
        iout=1 ;
        " Output will be appended to current output file ";
    elseif aa$=="GAUSS";
        " Output will be written to a GAUSS dataset. ";
        "     Enter a name for this GAUSS dataset: ";;oname=upper(cons);?;
         if nodatset(oname)==0; format 1,8 ;
             foo=g_yesno(" Dataset exists. Overwrite ","N");
             if foo==0; goto get1 ; endif;
         endif;
         iout=0;
    endif;

agin: ;
?;
   " Select:";
   if missrv(yhat2,0) ne 0 ;
      "    Whether to computed expected value of Y == E[Y] and " ;; $yhat2 ;
   else ;
      "    Whether to computed expected value of Y == E[Y] ";
   endif;
   "    Name of an ID variable, used to identify the observation ";
   if not(scalmiss(xns));
     "    The X variables to compute elasticities of. "; ?;
   endif;

   format 1,8;
   if missrv(yhat2,0) ne 0 ;
      iyhat=g_yesno("     Compute E[Y]  & " $+ yhat2,iyhat);
   else ;
     iyhat=g_yesno("     Compute E[Y] ",iyhat);
   endif;
   {foo,idname}=  g_askv1(
                    "      Enter an 'ID' variable (0=Don't include):",IDNAME,
                    getname(dataset)|"0",
                    "   ....  Select from " $+ dataset $+ " ",0,0);
   if scalmiss(foo);
      idindx=0; idname="0";
    else;
      if idname$=="0" ;
         idindx=0; idname="0";
      else ;
         idindx=foo;
      endif;
   endif;

   cls ;
   if scalmiss(xns);
      elids=miss(1,1);
     " No independent variables specified. ";
   else;
                 "      Enter variables to compute elasticities for .. ";
      curlist=g_askvs(curlist,xns);?;?;?;?;?;?;?;?;?;
      if curlist$=="0";
         elids=miss(1,1);
         curlist=0;
      else;
         elids=indcv(curlist,xns);
      endif;
   endif;

   ovars=miss(1,1);
   if idindx > 0 ;
      ovars=idname ;
   endif;
   if iyhat > 0 ;
      ovars=ovars|"YHAT" ;
      if missrv(yhat2,0) ne 0 ;
         ovars=ovars|yhat2 ;
      endif;
   endif;
   if  not(scalmiss(elids));
       ovars=ovars|curlist ;
   endif;
   ovars=packr(ovars);

@ nothing selected, so duck out @
   if scalmiss(ovars);
      retp(0);
   endif;

   format 1,8;
   " Selections: ";
   "     Predictions: "  $yesno[iyhat+1] ;;
         if missrv(yhat2,0) ne 0 ;
               " (and ";;format 1,8; $yhat2 ; else ; ? ;
         endif;
   "     Id variable:  ";;
      if idindx == 0; " None " ; else; $idname ; endif ;
   "     X variables:  ";;
      if scalmiss(elids); " None " ; else ; $xns[elids]' ; endif;
     foo=g_yesno(" Is this okay ","Y");
     if foo==0; goto agin ; endif ;

   if iout==0;
      create outf=^oname with ^ovars,0,8 ;
   endif;

   foo=seekr(f1,1);

   if iout==1 ;
      output on ; screen off ;
      ?;?; "-------------------------------------";
      if idindx > 0;
        "   Id Var   ";;
      endif;
      if iyhat> 0;
        "    E(Y)      " ;;
         if missrv(yhat2,0) ne 0 ;
               format 8,8; $yhat2 ;;
         endif;

      endif;
      if not(scalmiss(elids));
        "  Elastciticies  ";
      else;
         ?;
      endif;
@      format 10,8 ; $ovars';
      ?;@
      output off; screen on ;
   endif;


   ? ;? ;
   ir=csrlin; ic=csrcol ; allobs=rowsf(f1); idid=0;

   do until eof(f1);
      {elasts,yhat,idvals}=preds(beta,idindx);

      if scalmiss(yhat);        @ no valid obs @
         goto endloop ;
      endif;

      dtaout=miss(1,1);
      if idindx > 0 ;
         dtaout=idvals ;
      endif;
      if iyhat > 0 ;
        if not(scalmiss(dtaout));
          dtaout=dtaout~yhat ;
        else;
          dtaout=yhat ;
        endif;
      endif;
      if not(scalmiss(elids));
         if scalmiss(dtaout);
           dtaout=elasts[.,elids] ;
         else;
           dtaout=dtaout~elasts[.,elids];
         endif;
     endif;

     if iout==1 ;
       screen off ;
       output on ;  format /rz /m1 10,5 ; dtaout ; output off ; screen on ;
     else;
        if writer(outf,dtaout) /= rows(dtaout) ;
           errorlog " Disk Full ";
           if iout==0 ; foo=close(outf); endif ;
           retp(0);
        endif;
     endif;

endloop: ;
     if not(eof(f1));
       locate ir,ic ;
       idid=idid+rows(dtaout);
       format /ld 10,0 ;
       " Completed " idid " ::  Enter A to abort..." ;
       foo=key;
       if foo==97 or foo==65 ;
          if iout==0 ; foo=close(outf); endif ;
         retp(2);
       endif;
     endif;

  endo;

  if iout==0 ; foo=close(outf); endif ;
  retp(1);

 endp;



/************************/
@ Compute CV. Note that we return CV even for parameters whose value equals
   a constraint.  This is probably not technically correct . @

   proc (2)=compcvc(b,vctype,freevar0,&model);
       local cv,tmp,hmtx,bhhh,hessnr,isbad ;
       local model:proc ;

       if ismiss(freevar0);            @if no 'constraints' @
          freevar0=seqa(1,1,rows(b));
       ELSE ;
           IF ( MINC(MISSRV(FREEVAR0,0)) LE 0 ) OR
              ( MAXC(MISSRV(FREEVAR0,0)) Gt ROWS(B)) ;
             if _error_out ne 0 ; output on ; endif;
              ERRORLOG " BAD FREEVAR LIST : ";
              " FREEVAR LIST= " ;; FORMAT /LZ 12,5 ; FREEVAR0';
              " B "; b' ;
             if _error_out ne 0 ; output on ; endif;
          " !! CV matrix unobtainable! Using Identity matrix ... ";
              cv=eye(rows(freevar0));
              hmtx=4 ;
              trap 0; output off ;
              call g_pause(5,0);
             retp(cv,hmtx);
          ENDIF;
       endif;


       trap 1; 
       if _error_out ne 0 ; output on ; endif;

eek: ;

    if vctype==1 ;           @ if hessian desired @
       isbad=0;
@       _CtCall[2]=_CtCall[2]+1 ;   _CtCall[1]=_CtCall[1]+1 ; @
       {tmp,tmp,cv}=model(b,2);
       if maxc(freevar0)>rows(cv) or minc(freevar0)<1 ;
          call grbl2_error("ERROR in COMPCVC: " $+ g_ftos2(rows(cv),1) $+ " rows, with rows(freevar0) = " $+ g_ftos2(rows(freevar0),1),1);
          stop ;
       endif;
       if not (badmtx(cv));
            cv=invpd2(-cv[freevar0,freevar0]);
       else ;
          isbad=1;
       endif;
       if scalerr(cv) or isbad==1;              @ can't, so use  BHHH matrix @
            if _error_out ne 0 ; output on ; endif;
             " Hessian uninvertible. Using BHHH matrix " ;
              vctype=2;
              goto eek ;
       endif;
       trap 0; output off ;
       retp(cv,1);

  elseif vctype==2;          @ bhhh desired, or neccessary @
       isbad=0;
@       _CtCall[2]=_CtCall[2]+1 ;       _CtCall[1]=_CtCall[1]+1 ; @
       {tmp,tmp,cv}=model(b,4);
       if not (badmtx(cv));
           if rows(cv) lt maxc(freevar0) ; format /ld 6,0;
		"ERROR: can not extract from VC matrix (of size ";;rows(cv);; " rows)";
		 isbad=1;
	    else;
                cv=invpd2(cv[freevar0,freevar0]);
            endif;
       else ;
          isbad=1;
       endif;
       hmtx=2;
       if badmtx(cv) or isbad==1 ;
            if _error_out ne 0 ; output on ; endif;
          " !! CV matrix unobtainable! Using Identity matrix ... ";
          cv=eye(rows(freevar0));
          hmtx=4 ;
       endif;

       trap 0; output off ;
       retp(cv,hmtx);

  elseif vctype==3;          @ robust 'white' estimator @
@        _CtCall[2]=_CtCall[2]+1 ;        _CtCall[1]=_CtCall[1]+1 ; @
         {tmp,tmp,hessnr}=model(b,2);
         cv=invpd2(-hessnr[freevar0,freevar0]);
         if scalerr(cv);                   @ can't, so use  BHHH matrix @
            if _error_out ne 0 ; output on ; endif;
             " Can't compute 'White' matrix -- using BHHH instead ";
              vctype=2;
           output off ;
              goto eek ;
         endif;

@ got hessian -- build 'white' matrix @
@         _CtCall[2]=_CtCall[2]+1 ;      _CtCall[1]=_CtCall[1]+1 ;  @
         {tmp,tmp,bhhh}=model(b,4);
          cv=cv*(bhhh[freevar0,freevar0])*cv ;
          if not(diag(cv)>0);
              " Hessian unacceptable (- elements on diagonal). Using BHHH";
              vctype=2;
           output off ;
              goto eek;
          endif;
          trap 0; output off ;
          retp(cv,3);               @ 1 signals 'used hessian' @

   endif;

 endp;

/********************************/
@ Compute stderr and tstat, given vc, beta, and freeindx.
  Note that tstat & stderr will be set to missing if freeindx doesn't point
  to a variable (freeindx= . or 0 means 'include all')  @

  proc (2)=do_tstat(beta,vc,freeindx);

     local stderr,tstat,tmp;
     freeindx=packr(freeindx);
          if missrv(freeindx,0) ==0 ;
                 stderr=sqrt(diag(vc));
                 tstat=beta./stderr ;
           else ;
                 tstat=miss(zeros(rows(beta),1),0); stderr=tstat ;
                 tmp=sqrt(diag(vc)); stderr[freeindx]=tmp ;
                 tstat[freeindx]=beta[freeindx]./stderr[freeindx] ;
           endif;
     retp(stderr,tstat);
  endp;



/***************/
@ proc to ask user for desired hill climbing technique @
@ the 6th argument is for user selected method, presumably NOT using
 MAXLIKC (note that 7 is random search ) @

   proc (3)=askmlalg(arg6,dtype,htype,steptype);

     local ans,ahh,hh,ans2,ans6,ii6,ans6a,ansvc,yesno,strlist;
     local tmp,ii,atmp,chnglst,foo,parmlist,messlist ;
      local atech,ith,tlist,ie1,ie2,ie3,ie4,ie5,a1,a2,a3,a4,a5,isfile;

@ note: mlmessxx initialized in declare statements above @

     let ansvc=HESSIAN BHHH WHITE ;
     let ans6a=N H S D G O R ;
     let ans2=DONE MAX_ALG GOLDEN MLE_LIST VC SET_OPT ;
     let yesno=NO YES ;

     let chnglst="__rows"  "_MAXITER" "_MLEPSG" "_MLEPSB" "_MLATMAX"
                 "_MLMNSTP" "_MLBHHH"  "_MLSTEEP" "_MLMARQ" "_USEALLB"
                 "_MLBKUP" "_MLDORND" "_MLDORN2" "_MLRNDVC" "_MLVRYLO" ;

    let strlist=DONE NONE SHOW SET HELP ;

    let string messlist=   " : # of obs per read, 0=Program determines "
   " :  Maximum # of iterations "
   " : Stopping criteria, gradient."
   " : Stopping criteria, beta ."
   " : # criteria to satisfy: 1,2,3 (dl is 3rd)."
   " : Minimum step size in non-Golden search."
   " : 1=Use BHHH, if required."
   " : Use steepest descent, if required. 0=No, 1=Yes, 2=Yes & 1-by-1 if needed."
   " : 1=Use Marquardt method (psd hessian). "
   " : Do constrained max: 0=No, 1=Yes, 2=Best. "
   " : 1=Try moving backwards if stuck. "
   " : # of random á, if necessary. 0=None."
   " : min # of successful random searches."
   " : VC for random method: <1.0=Fraction of á,  1=Hessian, 2= BHHH."
   " : Set to 500 to force convergence if dl very small, 1=standard";

@ do some less generic parameters @
   if not(scalmiss(typecv("_doelast")));
      messlist=messlist$|" : 1=Automatically compute elasticities.";
      chnglst=chnglst|"_doelast" ;
    endif;
    if not(scalmiss(typecv("_dotest")));
         chnglst=chnglst|"_dotest" ;
         messlist=messlist$|" : 1=Compute a specification test ";
    endif;

     foo=0; parmlist=zeros(rows(chnglst),1);

@ fill in globals @
     do until foo==rows(chnglst);
        foo=foo+1;
        parmlist[foo]=varget(chnglst[foo]);
     endo;

     tlist=varget("_mlelist");
     isfile=0;

     if scalmiss(arg6);
        ii6=6;
        let ans6=NR BHHH STEEPEST DFP BFGS RANDOM ;
     else ;
        ii6=7;
        let ans6=NR BHHH STEEPEST DFP BFGS OTHER RANDOM ;
        ans6[6]=arg6 ;
     endif;
     dtype=minc(ii6|dtype);

     ahh="HELP";


 do until 1==2 ;                        @ endless loop @
    ? ;
    ahh=g_askopt(" Change algorithim, etc...  (ESC to exit)",ans2|"HELP",
                 "DONE",0,3,3);
    if scalmiss(ahh); retp(dtype,htype,steptype); endif;

    if ahh$=="HELP";
cls;?;
call boldon ;"   Select one of the following options: ";
   call normalon ;
"        DONE) No changes ";
"     MAX_ALG) Maximization algorigthim            (Currently= " $ans6[dtype] ;
"      GOLDEN) GOLDEN step-search (NO=stepsize of 1.0)  (Currently= "
                                                      $yesno[steptype+1] ;
"    MLE_LIST) Custom control of maximization strategy (Currently:";;
       if cols(tlist)<5 ; " NO "; else ;" YES"; endif;
"          VC) VC matrix                           (Currently= " $ansvc[htype];
"     SET_OPT) Set optimization-control parameters  ";?;
       continue ;
    endif;

@ jump to proper section @
    hh=indcv(ahh,ans2)-1;
    if hh==0 ; ?; retp(dtype,htype,steptype);
      elseif ahh$=="MAX_ALG" ; goto a1 ;
      elseif ahh$=="GOLDEN" ; goto ss1;
      elseif ahh$=="MLE_LIST" ; goto setlist ;
      elseif ahh$=="VC" ; goto ggh;
      elseif ahh$=="SET_OPT" ; goto sett;
    endif;

/* SET ALG */
a1:      if dtype<1 or dtype>ii6; dtype=1 ; endif;
        foo=g_askopt(" Select maximization algorithim",ans6|"HELP",dtype,0,3,3);
        if scalmiss(foo); retp(dtype,htype,steptype); else ; dtype=foo; endif;
        if dtype$=="HELP";
         cls ;
          ?;  call boldon ;
         "    Select hill climbing alg; as in the D in B2=B1- step*inv(D)*G ";
              call normalon ;
           "           NR  ) Newton Raphson ";
           "         BHHH  ) BHHH " ;
           "      STEEPEST ) Steepest Descent ";
           "          DFP  ) DFP    ";
           "         BFGS  ) BFGS  " ;
           "       RANDOM  ) Random Search ";
           if not(scalmiss(arg6));
              "   or " $ans6[6]  "  ) " $arg6;
           endif;
           goto a1 ;
         endif;
         if dtype$=="RANDOM";
           dtype=7;
         else;
            dtype=indcv(dtype,ans6);
         endif;
         if dtype==7 ;
  "                       Note: You may want to set _MLRNDVC and _MLDORND ";
          endif;

          continue ;

/* SET STEP SIZE */
ss1:     ?;
         "        Select 'step size' algorithim  .. ";
      steptype=g_yesno("            Use GOLDEN Step Size ",steptype);
          continue  ;

/* SET VC MATRIX */
ggh:    foo=g_askopt(" Method to compute covariance matrix ",
                        ansvc|"HELP",htype,0,3,3);
        if scalmiss(foo); retp(dtype,htype,steptype); else ; htype=foo;endif;
        if htype$=="HELP";

       ?; call boldon ;
        "      Select method of computing covariance matrix: "; call normalon;
           "           HESSIAN) Use Hessian          (H) ";
           "              BHHH) Use the BHHH matrix  (G'G)";
           "             WHITE) Use robust 'White'   (Inv(H)*(G'G)*Inv(H))";
           goto ggh;
        endif;
        htype=indcv(htype,ansvc);
        continue ;

/* SET CONTROL PARAMTERS */
sett: ;
          cls; ?; call boldon;
          " Set 'optimization control' parameters...";; call normalon;?;
          parmlist  = g_usrset(" Set 'optimization control' parameters...",
				chnglst,parmlist,messlist,3,1,0);
          foo=0;
          do until foo==rows(parmlist);
            foo=foo+1;
            call varput(parmlist[foo],chnglst[foo]);
          endo;
          continue ;

/* SET _MLELIST */
setlist: ;
        foo=g_askopt(" Specify maximization strategy:",strlist,"HELP",0,3,3);
        if scalmiss(foo); retp(dtype,htype,steptype); endif;
@ Return to ALGSET menu @
        if foo$=="DONE" ;
                continue ;

        elseif foo$=="NONE" ;
           if cols(tlist) == 5 ; format /ld 5,0 ;
            " You have specified " rows(tlist)  " techniques.";
             foo=g_yesno("  Are you sure you want to cancel them","YES");
           else ;
              foo=1;
           endif;
           if foo==1 ;
                call blockon ; " There will be no user-defined techniques. ";;
                call normalon ; ?;
                call varput(0,"_mlelist");
                tlist=0;
           endif;
           goto setlist ;


        elseif foo$=="SHOW";
           call blockon ; " Current user-defined techniques are: ";
           call show_tch(tlist,dtype,steptype);
           goto setlist ;

        elseif foo$=="HELP" ;
           call helptech ;
           goto setlist ;

        endif;


@ Else, enter new ones @

   " Enter techniques, 1 per line. ENTER to stop. ESC to Cancel. ";
   ith=1 ;

   do until ith==0;             @ EXIT via BREAK @

     if isfile==0 ;
        format /ld 3,0 ;" Enter technique # " ith  ":";;
        atech=upper(g_cons2()) ; ?;
     endif;

   if scalmiss(atech) ; atech="" ; endif;
   if atech$=="" ;      @ DONE @
      if cols(tlist)<5;
           " Retaining prior list .... ";
      else ;
         if ith==1 ;
            " Retaining prior list....   ";
         else ;
            " Number of user defined techniques: " rows(tlist);
             format /ld 4,0;
             call varput(tlist,"_mlelist");
         endif;
      endif;
      goto setlist ;

   else ;

      if isfile==1 ;       @ FILE IO @
         {a1,ie2}=g_getwrd(atech,ie1+1);
         {a2,ie3}=g_getwrd(atech,ie2+1);
         {a3,ie4}=g_getwrd(atech,ie3+1);
         {a4,ie5}=g_getwrd(atech,ie4+1);
         {a5,ie1}=g_getwrd(atech,ie5+1);
         if (ie1*ie2*ie3*ie4*ie5==0);
            isfile=0;
            continue ;
          endif;

      else ;                @ CRT IO @

         {a1,ie1}=g_getwrd(atech,1);
         if a1$=="FILE";
            {a2,ie2}=g_getwrd(atech,ie1+1);
            if not(files(a2,0)$==0);
               atech=upper(getf(a2,0));
               atech= g_zapcmt(atech,"@");
               isfile=1; ie1=0;
            else;
               " No such file. ";
            endif;
            continue ;

         elseif a1$=="DOS";
             " Exit to DOS.  Enter EXIT to return to program ... ";
             dos ;
             continue ;
         endif;
         {a2,ie2}=g_getwrd(atech,ie1+1);
         {a3,ie3}=g_getwrd(atech,ie2+1);
         {a4,ie4}=g_getwrd(atech,ie3+1);
         {a5,ie5}=g_getwrd(atech,ie4+1);
         if (ie5*ie4*ie3*ie2)==0;
             " Not enough arguments "; continue ; endif;

       endif;

       tmp=a1 ;
       if a1$=="RANDOM";
          a1=7 ;
       else ;
         a1=indcv(a1,"NR"|"BHHH"|"STEEPEST"|"DFP"|"BFGS"|"DONE");
         if scalmiss(a1); errorlog " Bad technique: " $+ tmp ; continue; endif;
         if a1>5 ;a1=0; endif;       @ if DONE @
       endif;

       a2=indcv(a2,"BACKWARD"|"ONE"|"GOLDEN"|"ON");
       if scalmiss(a2); errorlog " Bad steptype: " $+ a2 ; continue; endif;
       a2=a2-2 ;

      if g_isnumber(a3); a3=stof(a3); else;
         errorlog " Bad frequency: " $+ a3 ; continue; endif;
      endif;

      if a1==0 ; a4=0;a5=0;
      else ;
         if g_isnumber(a4); a4=stof(a4); else;
           errorlog  " Bad lower iteration range: " $+ a4 ; continue;
         endif;

         if g_isnumber(a5); a5=stof(a5); else;
            errorlog  " Bad upper iteration range: " $+ a5 ; continue;
         endif;
      endif;

      if ith==1 ;
          tlist=a1~a2~a3~a4~a5;
      else ;
          tlist=appnd_it(tlist,a1~a2~a3~a4~a5);
      endif;

      format /ld 5,0 ;  if isfile==1; " Technique: " a1 a2 a3 a4 a5 ; endif;
      ith=ith+1;
   endo;
   goto setlist ;
endo;
endp;

/*********************************************/
@ help on how to enter technique list @

proc helptech ;
cls;
call boldon ;
"  Enter technique list."; call normalon ;
"  ";; call blockon ;
"Default techniques";; call normalon ;
call boldon ;
" (selected with MAX_ALG and SET_OPT) are done AFTER these.
"; call normalon ;
"  To suppress this (and use ";;
call normalon ; call blockon ; "default techniques";;call normalon;
" only), select NONE.

  You enter 5 components: Algorithim, Steptype, Frequency, Min_iter, Max_iter.
  Algorithim:  DONE NR BHHH STEEPEST DFP BFGS RANDOM
                (Note: DONE stop's estimation)
    Steptype:  BACKWARD ONE GOLDEN ON
                (Note: BACKWARD to move in -direction. ON used with DONE only.)
   Frequency:  n (a number):
                 n=1   do every iteration
                 n>1   do every nth iteration
                 Less frequent methods should precede more frequent methods.
   Min_iter & Max_iter == When to do it (inclusive)
 (hit any key to continue.. )";;
call g_pause(30,-27);
"
Alternatively, you can enter a filename containing these commands. To do this,
simply enter " ;; call boldon ;
" FILE filename.ext";; call normalon ;
".  The techniques in the file will be ADDED";
"to the list. Note: enter DOS for temporary exit to DOS.

Example: ";
call boldon ; "   BHHH ONE 1 1 1      ";; call normalon ;
"     (start with BHHH)";
call boldon ; "   FILE NREVERY5.TCH  ";; call normalon ;
"   (contains several NR techniques for periodic use)";
call boldon ; "   DFP GOLDEN 1 1 100000     ";; call normalon ;
" (DFP GOLDEN every iteration)";
call boldon ; "   STEEPEST ONE  1 1 100000   ";; call normalon ;
" (Steepest every iteration -- a backup technique)";
call boldon ;
"   STEEPEST BACKWARD 1 1 20 ";; call normalon ;
" (  another backup technique, but ";
"                                 only try it during first 20 iterations) ";
call boldon ; "   DONE ON 1 1 10000  " ;; call normalon ;
"           (If none of the above work, stop estimation) ";
" (hit any key to continue.. )";;
call g_pause(30,-27);

 retp(0);
endp;



/***** Show a list of current techniques ****/
proc show_tch(Tlist0,hessalg,dogolden);
 local ms,ss,tmp,im,techlist,i1,tlist;

   ms="NR"|"BHHH"|"STEEPEST"|"DFP"|"BFGS"|"DONE"|"RANDOM"|"STEEP_1";
   ss="BACKWRD2"|"BACKWARD"|"ONE"|"GOLDEN";

     tmp=_mlelist ; _mlelist=0;
     techlist=mk_techs(hessalg,dogolden,2) ;
     _mlelist=tmp ;

     i1=1 ;
     ? ; cls ;

foo1: ;
     if i1==3 ; retp(1); endif;       @done @

     if i1==2 ;  tlist=techlist ; call normalon ;
         ? ; " Displaying ";; call blockon ; "Default techniques.";
         call normalon ;  ?;
     endif ; @ display default techniques @

     if i1==1 ;  tlist = tlist0 ; call normalon ;
        ? ; " Displaying ";; call blockon ;
        "User defined (step-by-step) techniques.";
        call normalon ; ?;
    endif ;  @ display stepbystep techniques @


    if cols(tlist) < 5 ;
         "   .... None specified. " ;
         call normalon ;
         retp(0);
     else ;
        im=1 ;
        do until im>rows(tlist);
            tmp=tlist[im,1]; if tmp==0;tmp=6; endif;
           " Technique # ";; format /ld 4,0 ; im  " : " ;; format 1,8;
            call boldon ; $ms[tmp] ;;
            tmp=tlist[im,2]+3 ;
            if tmp<4 ; $ss[tmp] ;; else ; "    " ;; endif;
            call normalon ;
            tmp=tlist[im,3:5]; format 4,0 ;
            if tmp[1]==0;
              call boldon ; " first iter " ;; call normalon ;
            elseif tmp[1]<0 ;
               call boldon ;   -tmp[1] " iter only ";; call normalon;
            elseif tmp[1]==1 ;
               call boldon ;  " every iter ";; call normalon ;
            else ; " every " ;; call boldon; tmp[1] ;; call normalon ;
                  " iter ";;
            endif;
            format 4,0 ; " (range: " ;;
            call boldon; tmp[2] ;;
            call normalon ; " to" ;;
            call boldon ; tmp[3] ;;
            call normalon; " )";
            if (im%20)==0 ;
                call g_pause(20,-1);?;
            endif;
            im=im+1;
        endo;
        " (hit any key to continue..)";; call g_pause(20,-27);?;
 endif;
 i1=i1+1;
 goto foo1;
endp;


/*********************************/
@ process some grbl, mle, batch commands @
     proc  (6)=grblbtch(keywurd,acommand,ie,ditype,steptype,
                        vctype,bvarsb,cnststrg);

    local aword,nn,ok,aword2,aword3 ;
    ok=0;

    if keywurd$=="STRATEGY" or keywurd$=="METHOD" or keywurd$=="ALG" ;
              ok=1;
              {aword,ie}=g_getwrd(acommand,ie+1);
              if ie == 0 ; goto done ; endif;

              if aword$=="ALG" or aword$=="METHOD" or aword$=="STRATEGY" ;

                 {aword,ie}=g_getwrd(acommand,ie+1);
                 if ie == 0 ; goto done ; endif;
                 if aword$=="BHHH" or aword$=="G'G" or aword$=="2";
                        ditype=2;  " Using BHHH technique ";
                 elseif aword$=="NR" or aword$=="HESSIAN" or aword$=="1";
                        ditype=1;  " Using NR technique ";
                 elseif aword$=="STEEP" or aword$=="STEEPEST" or aword$=="3";
                      ditype=3 ; " Using Steepest Descent " ;
                 elseif aword$=="DFP"  or aword$=="4";
                      ditype=4 ; " Using DFP " ;
                 elseif aword$=="BFGS"  or aword$=="5";
                      ditype=5 ; " Using BFGS" ;
                 elseif aword$=="OTHER" or aword$=="6";
                      ditype=6 ; " Using non-mle technique " ;
                 elseif aword$=="RANDOM" or aword$=="7";
                     ditype=7 ; " Using RANDOM search " ;
                 endif;

             elseif aword$=="STEP" or aword$=="STEPTYPE" or aword$=="GOLDEN" ;
                   {aword,ie}=g_getwrd(acommand,ie+1);  @ grad a parameter@
                   if ie == 0 ; goto done ; endif;
                   if aword$=="1"  or aword$=="ONE" or aword$=="NO";
                            steptype=0;  " Not using GOLDEN ";
                   elseif aword$=="GOLDEN" or aword$=="YES" ;
                            steptype=1;   " Using GOLDEN ";
                   endif;

             elseif aword$=="VCTYPE" OR aWORD$=="CVTYPE" or
                         aword$=="VC";
                   {aword,ie}=g_getwrd(acommand,ie+1);
                   if ie == 0 ; goto done ; endif;
                   if aword$=="BHHH" or aword$=="G'G";
                         vctype=2 ; "  Using BHHH for VC matrix ";
                   elseif aword$=="NR" or aword$=="HESS" or
                                                     aword$=="HESSIAN";
                         vctype=1;   " Using Hessian for VC matrix ";
                   elseif aword$=="WHITE" or aword$=="HETERO";
                           vctype=3;
                   " Using 'White' heteroscedastic consistent matrix for VC ";
                   endif;

              endif;



    elseif keywurd$=="CONSTR" or keywurd$=="CONSTRAINT" or
                 keywurd$=="CONSTRAINTS" ;
                 ok=1;
                  {aword,ie}=g_getwrd(acommand,ie+1);
                  if ie == 0 ;
                       cnststrg=miss(1,1);
                  else ;
                     if files(aword,0)$/= 0 ;
                        cnststrg=getf(aword,0);
                        " Constrained estimation using " $aword ;
                     else ;
                        errorlog " No such CONSTRAINTS file: " $+ aword ;
                        cnststrg=miss(1,1);
                     endif;
                  endif;

     endif;
done:     retp(ditype,steptype,vctype,bvarsb,ok,cnststrg);
   endp;


/****************************************/
@ Process a string into a matrix of constraints-descriptors. This matrix can be given to MAXLIKC
  (as the 2nd arghment)
  
    Alist: string containing constraints.
	   each constraint is seperated by a ;
           constraints have the form:  n1 * var1 + n2 * var3  OPER n3
               n1.. are numbers
               * + and = are optional.  The important thing is the order.
               OPER can be one of: LE LT < GT GE > EQ
               Note: LE and GE are converted to LT and GT.
   
   Parlist: list of parameter names. MUST match order of PARAMETERS sent to maxlikc.
   EXACTXYZ: if 1, parameters names must exactly match parlist names. 0 means 'substring match is okay'

  Returns a #constraints x (K+2) matrix, where K = rows(parlist)

  The first K terms are multipliers on the ParList variables. By default (for variables not
  mentioned in a constraint), a multiplier of zero is used.

  The K+1 term is a "condition" flag: 1= lt, 2=eq.
  Note: gt constraints are changed to negative lt constraints
         ge and le  converted to gt and lt

  The K+2 term is the value (that the constraint should be LT, or EQ, to).

@

proc mle_constraints(alist,parlist,exactXYZ);


local nk,oldlist,cnmtx,aword1,ie,idx,aval,acond,condlist,idone,cnmtx2 ;
local gotval,isminus,adiv,foo,tmp,tmp2,ieq1,tmult,newlist ;


@ corrective hack: remove empty phrases @
newlist="";
tmp=alist;
do until 1==2 ;
  if tmp$=="" ; break ; endif;
  nk=strindx(tmp,";",1);
  if nk==0;  
     newlist=newlist $+ tmp $+ " ;  " ;
     break ;
  endif;
  aword1=g_strip(strsect(tmp,1,nk-1));
  if aword1$/="" ; newlist=newlist $+ aword1 $+ " ; " ; endif;
  tmp=g_strip(strsect(tmp,nk+1,1000000));
endo;
alist=newlist;  

alist=g_strip(alist);

let condlist=LT "<" GT ">" EQ "=" LE GE   ;

@ Use mle_ask_constraints to enter constraints from keyboard @
    if scalmiss(alisT);
       retp(miss(1,1));
   endif;

   nk=rows(parlist);
   alist=upper(g_zapcmt(alist,"@"));
   idx=0;idone=0;
   cnmtx2=miss(1,1);
   {aword1,ie}=g_getwrd(alist,1); aword1=g_strip(aword1);

   do until ie==0;

@ 4 possiblities: condition, +, number, or variable @
        IDX=0; aval=1; gotval=0; isminus=1 ;
        cnmtx=zeros(1,nk+2) ;
        DO UNTIL ie==0 ;        @ READ IN THIS CONSTRAINT EXIT VIA  break@
          ACOND=indcv(aword1,condlist);
          if not(scalmiss(acond));                 @ 1-- a condition, @
             if idx==0 and gotval==1 ;
                errorlog " No variable"; acond=miss(1,1);
             endif;
             BREAK ;
          ELSEif aword1$=="+" or aword1$=="-";   @ 2--GET NEXT VAL/VAR PAIR @
             if idx==0 and gotval==1 ;
                errorlog " No variable"; acond=miss(1,1); break ;
             endif;
             if aword1$=="-" ; isminus=-1; else; isminus=1; endif;
             IDX=0; aval=1; gotval=0;
@ 3--a value @
          ELSEIF g_ISNUMBER(AWORD1) and aword1$/="E" and aword1$/="PI";
             AVAL=STOF(AWORD1); gotval=1;
             if idx ne 0 ;
                cnmtx[idx]=aval*isminus;
             endif;
           else ;                       @ 4-- a variable @
                if exactxyz==1 ;
                  idx=indcv(aword1,parlist);
                else;
                  idx=g_indcv2(aword1,parlist);
                endif;
                if scalmiss(idx);
                    errorlog  " Bad variable: " $+aword1 ;
                    acond=miss(1,1); break ;
                elseif idx ne 0 ;
                    cnmtx[idx]=aval*isminus ;
                endif;
            endif;
            {aword1,ie}=g_getwrd(alist,ie+1);
         endo;

@ when here, either a condition, bad constraint, or ie=0 @
         if ie==0 ; break ; endif ;
         if scalmiss(acond);  gosub skipit ; continue ; endif ;

@ else, ok constraint @

         if aword1$=="EQ" or aword1$=="=" ;
               if sumc((cnmtx[1:nk]./=0)') gt 1 ;
                  errorlog " More than 1 variable in equality constraint!";
                  gosub skipit ;
                  continue ;
               else;
                  adiv=sumc(cnmtx[1:nk]');
                  cnmtx[nk+1]=2;
               endif;
         elseif aword1$=="GT" or aword1$=="GE" or aword1 $==">" or
                aword1$==">=" or aword1$=="=>" ;
                cnmtx[nk+1]=3;         @ a hack, used below @
         elseif aword1$=="LT" or aword1$=="LE" or aword1 $=="<" or
                aword1$=="<=" or aword1$=="=<" ;
                cnmtx[nk+1]=1;
         endif;

@ Get value @
         {aword1,ie}=g_getwrd(alist,ie+1); aword1=g_strip(aword1);
         if ie==0;  break ; endif;
	 tmult=1;
         if strsect(aword1,1,1)$=="-" ; tmult=-1; aword1=strsect(aword1,2,111) ; endif;
         if not(g_isnumber(aword1));
               errorlog " A Bad final value: " $+ aword1 ;
               continue ;
         else;
               cnmtx[nk+2]=tmult* stof(aword1) ;
               if cnmtx[nk+1]==2 ;        @ eq constraint @
                  cnmtx[nk+2]=cnmtx[nk+2]/adiv ;
                  cnmtx[1:nk]=cnmtx[1:nk]./=0 ;
               elseif cnmtx[nk+1]==3 ;     @ gt constraint .. - lt constraint@
                  cnmtx=-cnmtx;cnmtx[nk+1]=1 ;
               endif;
               if scalmiss(cnmtx2) ;cnmtx2=cnmtx ;
                 else; cnmtx2=cnmtx2|(cnmtx) ;
               endif;
         endif;
         gosub skipit ;              @ find ; or ie==0 @
     endo;                             @ get next constraint @

@ no more constraints @

     retp(cnmtx2);

@ search for ; @
skipit: ;
      do until ie==0 ;
         {aword1,ie}=g_getwrd(alist,ie+1);
         if aword1$==";";
            if ie ne 0 ;
              ieq1=ie;
              {aword1,ie}=g_getwrd(alist,ie+1);
              if g_strip(aword1)$==";"; ie=ieq1; continue; endif;
            endif;
            return;
         endif;
      endo;
      return ;
 endp;

/**********************************/
@ procedure to ask user to enter constraints @

proc mle_ask_constraints(parlist) ;
local foo,cnststrg ;

oy: ;
" Enter name of file containing constraint list,  
  or hit the ENTER key to enter from the keyboard .... ";
 call boldon ;"        ?  ";; call normalon ; foo=upper(cons);?;


 if foo$/="" ;
    if files(foo,0)$/= 0 ;
       cnststrg=getf(foo,0);
       retp(cnststrg);
    else ;
         " No such file: " $foo ;
       goto oy;
    endif;
endif;

@ otherwise, CRT input @


?;   " Enter a list of constraints. ESC when done. ? for further help ";

    cnststrg=" ";
    do until 2==1 ;             @ exit via RETP @
        call blueon ;call boldon ;" ? (ESC when done)> ";;call normalon ;
        foo=upper(g_cons2()) ; ?;
        if scalmiss(foo) ;      @ ESC @
             retp(cnststrg);

        elseif type(foo)==6 ;           @ ignore arrow keys, etc @
            continue ;

        elseif foo$=="?";
?;
"Enter a list of constraints. Each constraint consists of coefficients,
parameter names, arithemtic operations, relational operator, and a value.
For example:";
call boldon ;
"     2 * X1  +  Z1 * 3 < 10  ; X2 = 3 ;
      3  X1 - Z3 > 100 ;
";
call normalon ;
"Constraints are separated by semi-colons.  You can have several constraints
on a single line, or you can spread one constraint over several lines.";
?;
call boldon ; " Notes:";call normalon;
"  *  Before hitting ESC, hit ENTER key (any constraints on 'current line'
     will be ignored).
  *  To ease entry, you don't need to include the * (e.g.; 3 X1).
  *   Coefficient values of 1.0 are assumed when no coefficient is specified
      (e.g.; Z3 in 3rd constraint).
  *  You MUST leave a space between each item in the constraint.
     In other words: ";; call boldon ;" 2*X1+3*Z1<10";; call normalon;" is NOT permitted! ";

          call g_pause(20,-1);
          if not(scalmiss(parlist)) ; format 1,8;
              " Current parameters: "; $parlist' ;?;
          endif ;

          continue ;

        else ;
           cnststrg=cnststrg $+ " " $+ foo ;
        endif;

   endo ;
   endp ;


/************************/
@  Display constraints @

proc mle_show_constraints(cnstlist,parlist);

    local k,ithc,i,igot;

    if  scalmiss(cnstlist);
        " No constraints. ";
        retp(0);
    endif;

@ else, constraints... @
   k=rows(parlist);
   ithc=0;
   do until ithc == rows(cnstlist);            @ for each constraints @
     ithc=ithc+1;
     i=0;
     igot = 0;
     do until i==k;           @ scan each parameter @
        i=i+1;
        if cnstlist[ithc,i]==0 ;
           continue ;
        else ;
            if igot==0 ;
              format /lz 10,5 ;
             " ( " ;; cnstlist[ithc,i] ;;
            elseif cnstlist[ithc,i]<0  ; format /lz 10,5 ;
              ")  -  ( " ;; abs(cnstlist[ithc,i]) ;;
            else ;
              ")  +  ( " ;; cnstlist[ithc,i] ;;
            endif;

            " * " ;; format 1,8 ; $parlist[i] ;;
            igot=1 ;
        endif;

     endo ;

     if cnstlist[ithc,k+1]==1 ; ")  <  " ;; else ; ") = " ;; endif ;
     format /lz 10,5 ; cnstlist[ithc,k+2] ;

  endo;

retp(1);
endp;


/*****/
@ a shortcut @
proc badmtx(aval);
local aa,tf;
if rows(aval)==1 and cols(aval)==1 ;
  aa=scalinfnanmiss(aval);
else ;
  aa=isinfnanmiss(aval);
endif;
retp(aa);
endp;


/*****/
proc nodatset(afile);
local foop,aa ;
open foop=^afile ;
aa=close(foop);
if foop<0 ; retp(1); else; retp(0); endif;
endp;