/* create a new dataset, that contains various ID, NEW, and DUMMY variables (and observation selection) */
@ note: grbl2_create.dec should be included at top of calling program.
        (it initializes globals that calling programs can set)
        This file (grbl2_create.inc) should be included at the bottom
--------------------------------@

g_create: ;

clear wordlist,errlist,funcs,val_list,isneg;  @ globals used by equ_tree @
clear node_oper,node_neg,node_left,node_right ;
clear nr_orig,nr_sel,nr_use,used_obs;
clear got1,got2,got3,ismtx;
clear warning_mess ;
clear save_mess ;
clear firstprint ;

madevars={};
warning_mess={};
save_mess={};

foo=g_create_makex(1);    @ read data from input file @

iid={};

if nx==0;
  call grbl2_error("ERROR: you did not specify any X variables (reminder: X must follow INPUT).",1);
  stop;
else;
  dat=x;
  nam2=xname;
endif;

nobs=rows(dat);

if nxdummy>0;
  tmp=g_Vudum(dinfo,filename,1,0,xnewinfo.names);
  nam2=nam2|tmp;
  dat=dat~xdummy ;
endif;

if missrv(idindex,0) ne 0;
  dat=dat~id ;
  if got1>0 ;
     nam2=nam2|"ID";
     iid=iid|"ID";
  else;
     nam2=nam2|idname;
  endif;
endif;

if missrv(id2index,0) ne 0;
  dat=dat~id2 ;
  if got2>0 ;
     nam2=nam2|"ID2";
     iid=iid|"ID2";
  else;
     nam2=nam2|id2name;
  endif;
endif;

if missrv(id3index,0) ne 0;
  dat=dat~id3 ;
  if got3>0 ;
     nam2=nam2|"ID3";
     iid=iid|"ID3";
  else;
     nam2=nam2|id3name;
  endif;
endif;

if missrv(id1index,0) ne 0 ;        @ PANELS? @
  pids=1~rows(id1) ;
  oldid=id1[1];
  for j (2,rows(id1),1);
    if id1[j] ne oldid ;
   	   pids[rows(pids),2]= j-1 ;
       pids=pids|(j~rows(id1));
       oldid=id1[j];
    endif;
  endfor;
  " Number of panels " ;;format /ld 6,0; rows(pids) ;; " (using " ;;
  if id1name$=="#" ; id1index ; else ; format /ld 1,8; $id1name ; endif;
else ;
    pids=1~rows(dat) ;
endif;

madevars={};

if rows(equ_names)> 0;            @ for each panel (or just 1 panel if no BY specified @
  "Processing NEW and NEW IF keyphrases...";

  dat0=dat;
  keepdat=ones(rows(dat0),1);       @ default is to keep a row @
  equ_mtx=keepdat ;                 @ stub for matrix storage (max rows is rows in dat0 @
  equ_mtx_names=""~0~0~0~0~0 ;      @name, transpose-flag,row row2 col1 col2 @

  for j  (1,rows(equ_names),1);         @ for each equation ===================== @

  isthen=zeros(rows(dat0),1); iselse=isthen;

@ check basic structure of equations @
  addv1=0 ; addv2=0; addv1_use=0 ; addv2_use= 0;
  dokeep1=0; dokeep2=0;    @ flag to add vars, and keep obs @
  doelse=0;drop1="" ; drop2="";
  aequc=equ_conds[j];
  ismtx=0;
  dosay1=0;dosay2=0;
  if aequc$/="" ; " ... if ...  " ;;$aequc ; endif;

@ some preliminaries for this equation ...@
  if aequc$==""  ;       @ is this a non-conditional (specified using new if ), or a matrix definition ? @
       aname=0 $+ equ_names[j];
       aequ=equ_equs[j];

       ismtx=equ_ismtx[j];

       if ismtx==1 ;
           " ...  creating matrix ...  " ;;$aequ ;
       else;
           " ...  creating ...  " ;;$aequ ;
       endif;
       tmp=indcv(aname,nam2|equ_mtx_names[.,1]) ;
       if not(scalmiss(tmp));
           call grbl2_error("ERROR: A NEW variable, " $+  aname  $+
           ", has the same name as an existing, or other NEW, variable (or matrix) ",1);
           stop;
       endif;
       addv1=1;
       vv1=miss(zeros(rows(dat0),1),0) ;    @ uncond, or then, variable: default value is missing @
       isthen=isthen+1;         @ always true @

  else;                @ conditional @
      aname1=equ_names[j];
      aname1=0 $+ aname1 ;    @ THEN condition @
      aequ=equ_equs[j];

      if aequ$==".SAY" or aequ$==".SAY_MTX";
         dosay1=1;      dosay_it1=equ_names[j] ;
	 if aequ$==".SAY_MTX" ; dosay1=2 ; endif;
     elseif aequ$/=".DROP" and aequ$/=".KEEP" ;
         tmp=missrv(indcv(aname1,nam2),0) ;
         addv1=1;                           @ a THEN variable was specified (even if it is never TRUE @
         addv1_use=tmp ;
         vv1=miss(zeros(rows(dat0),1),0) ;    @ uncond, or then, variable: default value is missing @
      elseif aequ$==".KEEP";
         dokeep1=1;
      else;                 @ .DROP @
         drop1=aname1;
         if drop1$/="*" ;
            tmp=indcv(drop1,nam2) ;
            if scalmiss(tmp);
               call grbl2_error("ERROR: A IF .. THEN DROP variable, " $+ drop1 $+   ", does not exist",1);
               stop;
            endif;
         endif;
      endif;

      aname2=equ_names2[j];                 @ ELSE condition @
      aequ2=equ_equs2[j];
      if aname2$/="" or aequ2$/="";
         doelse=1;
         aname2=0 $+ aname2 ;

         if aequ2$==".SAY" or aequ2$==".SAY_MTX" ;
	    dosay2=1;
            dosay_it2=equ_names2[j] ;
            if aequ2$==".SAY_MTX" ; dosay2=2 ; endif;

         elseif aequ2$/=".DROP" and aequ2$/=".KEEP" ;
             tmp=missrv(indcv(aname2,nam2),0) ;
             vv2=miss(zeros(rows(dat0),1),0) ;    @ default value is missing @
             addv2=1;                           @ a else variable was specified (even if it is never TRUE @
             addv2_use=tmp;
         elseif aequ2$==".KEEP";
           dokeep2=1;
         else;                      @ DROP @
            drop2=aname2;
            if drop2$/="*";
              tmp=indcv(drop2,nam2) ;
               if scalmiss(tmp);
                 call grbl2_error("ERROR: A IF .. ELSE DROP variable, " $+ drop2 $+  ", does not exist",1);
                 stop;
               endif;
            endif;          @ drop2 @
         endif;             @ aequ =.DROP @
       endif;               @ aname2 /="" @
    endif;                  @ preliminaries (aequc="") @

@ done with preliminaries. Evalute the equations @
@ 4 vectors used: vv1,vv2, isthen, iselse (each rows(dat0),1)@

@ if its a matrix definition, do it (no IF condition allowed in matrix definitions,
  BY defintion of matrices not supported (though you can create matricies-on-the-fly
  using BY selected rows)@
 if ismtx==1 ;               @ NEW MATRIX  @
    istran=0;
    vv=doaneq(aequ,nam2,dat0,0,equ_mtx,equ_mtx_names);          @ unconditional @
    if cols(vv)>100 ;
       if rows(vv)>100 ;format /ld 6,0;
         "ERROR: matrix too large (max columns is 100): ";; cols(vv);;
       endif;
       istran=1 ;
       vv=vv';
    endif;
    if rows(vv) gt rows(equ_mtx) ;
         call grbl2_error("ERROR: matrix too large (max rows is " $+ g_ftos2(ROWS(equ_mtx),1) $+ " : "  $+ g_ftos2(rows(vv),1),1) ;
         stop;
    endif;
    tmp=aname~istran~(cols(equ_mtx)+1)~(cols(equ_mtx)+cols(vv))~1~rows(vv);
    equ_mtx_names=equ_mtx_names|tmp;
    arf=ones(rows(equ_mtx),cols(vv));
    arf[1:rows(vv),1:cols(vv)]=vv;

    equ_mtx=equ_mtx~arf ;
    continue ;          @ next equation @
endif;

@ if here, evaluating NEW or NEW IF (not NEW MATRIX) @

  for jp0 (1,rows(pids),1);           @ for each BY group:  @

    pbase=pids[jp0,1];             @ first row of this panel @
    pend=pids[jp0,2];
    inpanel=1+pend-pbase;
    dat=dat0[pbase:pend,.];
    aequ_mtx=equ_mtx[pbase:pend,.];

    if rows(pids)>1 ; jpanel=jp0 ; else ;jpanel=0; endif;

@ create new variables (possibly conditionally), or drop variables or rows ... using EQu_EQUS  @

    if addv1==1;       @ is this a non-conditional (specified using new)? @
      vv=doaneq(aequ,nam2,dat,jpanel,aequ_mtx,equ_mtx_names);                @ unconditional @
      vv=vv[.,1];                   @ discard other columns @
      if rows(vv)==1  ; vv=ones(inpanel,1).*vv ; endif;
      if rows(vv) gt  inpanel ; format /ld 5,0;
            vv=vv[1:inpanel];
            warning_mess=warning_mess$|"Truncating rows returned by THEN phrase " $+ g_ftos2(rows(vv),1) $+
                                       " " $+ g_ftos2(inpanel,1) $+ " (panel " $+ g_ftos2(jpanel,1);
      elseif rows(vv) lt inpanel;           @ pad if necessary  @
         tmp=zeros(inpanel-rows(vv),1);
         vv=vv|(miss(tmp,0)) ;
      endif;
      vv1[pbase:pend]=vv ;
    endif;

    if addv2==1 ;
          vv=doaneq(aequ2,nam2,dat,jpanel,aequ_mtx,equ_mtx_names);         @ else equation @
          vv=vv[.,1];                   @ discard other columns @
          if rows(vv)==1  ; vv=vv.*ones(inpanel,1); endif;
          if rows(vv) gt  inpanel ; format /ld 5,0;
              warning_mess=warning_mess$|"Truncating rows returned by ELSE phrase " $+ g_ftos2(rows(vv),1) $+
                                       " " $+ g_ftos2(inpanel,1) $+ " (panel " $+ g_ftos2(jpanel,1);
              vv=vv[1:inpanel];
          elseif rows(vv) lt inpanel;           @ pad if necessary @
             tmp=zeros(inpanel-rows(vv),1);
             vv=vv|(miss(tmp,0)) ;
          endif;
          vv2[pbase:pend]=vv ;
    endif;

    if aequc$/="" ;
       format /ld 1,8;
       vv=doaneq(aequc,nam2,dat,jpanel,aequ_mtx,equ_mtx_names);   @ condition equation @
       vv=vv[.,1];
       if rows(vv)==1 ; vv=vv.*ones(inpanel,1); endif;

       if rows(vv) gt  inpanel ; format /ld 5,0;
              vv=vv[1:inpanel];
              warning_mess=warning_mess$|"Truncating rows returned by IF phrase " $+ g_ftos2(rows(vv),1) $+
                                       " " $+ g_ftos2(inpanel,1) $+ " (panel " $+ g_ftos2(jpanel,1);
       elseif rows(vv) lt inpanel;           @ pad if necessary @
             tmp=zeros(inpanel-rows(vv),1);
             vv=vv|(miss(tmp,0)) ;
       endif;
       isthen[pbase:pend]= vv./=0; ;  @ the TRUEs (use THEN equation) @
       iselse[pbase:pend]= vv.==0 ; @ the FALSEs (use else equation) @
     endif;

  endfor;                   @ values, and true/false index, across all panels @

@ do the SAYs @

  if (dosay1 ne 0 or dosay2 ne 0)  and firstprint==0 ;
     mkout=" =============================================================== "  $+ chrs(10|13);
     screen off; output on ;    $mkout ;	output off; screen on ;
    firstprint=1;
  endif;

  if dosay1 ne 0 ;
      if sumc(isthen)>0 ;
	  mkout="" ;
          do while 1==1 ;
	      if dosay_it1$=="" ; break ; endif;
              dosay_it1=g_strip(dosay_it1);
              if strsect(dosay_it1,1,1)$/="\"" ;   @ no quote, must be var name @
		  {aword,dosay_it1}=token(dosay_it1); aword=upper(aword);
		  itmp=indcv(aword,nam2);
		  if not(scalmiss(itmp)) ;
			datdo=dat[.,itmp];
		  else;
                       ii=indcv(aword,equ_mtx_names[.,1]);
		       if scalmiss(ii);
			  call grbl2_error("ERROR: no such variable or matrix "  $+ aword,1) ; stop ;
		       endif;
                       ic1=equ_Mtx_names[ii,3];
                       ic2=equ_Mtx_names[ii,4];
                       ir1=equ_Mtx_names[ii,5];
                       ir2=equ_Mtx_names[ii,6];
                       datdo=equ_mtx[ir1:ir2,ic1:ic2];
		  endif;
                  if rows(datdo)>1 ;
		     mkout=mkout $+ chrs(10|13);
		  endif;

                  blen=0;
                  if dosay1==2 ;
	            for ir (1,rows(datdo),1);
	            for ic (1,cols(datdo),1);
   	         	  foo=g_ftos2(datdo[ir,ic],0);
			  blen=maxc(blen|strlen(foo));
                    endfor;
                    endfor;
		    foo=blen*cols(datdo);
                    if foo>179;
			outwidth 255 ;
		    elseif foo> 100 ;
		        outwidth 180 ;
	            elseif foo>78 ;
			outwidth 120 ;
		    endif;
                 endif;            


	          for ir (1,rows(datdo),1);
		  if blen>0 ; mkout=mkout $+ "       " ; endif;
	          for ic (1,cols(datdo),1);
   	         	  foo=g_ftos2(datdo[ir,ic],0);
			  if blen>0 ; foo=g_padwrd(foo,blen,1); endif;
                          mkout=mkout $+ " " $+ foo; 
                  endfor;
                  if rows(datdo)>1 ;
		     mkout=mkout $+ chrs(10|13);
		  endif;
                  endfor;
                  continue ; 
             endif;		
	     i2=strindx(dosay_it1,"\"",2);
	     if i2==0 ;  mkout=mkout $+  strsect(dosay_it1,2,100000) ; break ; endif; @ close an unclosed " string @
             tt=strsect(dosay_it1,2,i2-2);
             mkout=mkout $+ tt ;
             dosay_it1=strsect(dosay_it1,i2+1,100000);
         endo;
	 screen off; output on ;    $mkout ;	output off; screen on ;
     endif;

  endif;


  if dosay2 ne 0 ;
      if sumc(iselse)>0 ;
	  mkout="" ;
          do while 1==1 ;
	      if dosay_it2$=="" ; break ; endif;
              dosay_it2=g_strip(dosay_it2);
              if strsect(dosay_it2,1,1)$/="\"" ;   @ no quote, must be var name @
		  {aword,dosay_it2}=token(dosay_it2); aword=upper(aword);
		  itmp=indcv(aword,nam2);
		  if not(scalmiss(itmp)) ;
			datdo=dat[.,itmp];
		  else;
                       ii=indcv(aword,equ_mtx_names[.,1]);
		       if scalmiss(ii);
			  call grbl2_error("ERROR: no such variable or matrix "  $+ aword,1) ; stop ;
		       endif;
                       ic1=equ_Mtx_names[ii,3];
                       ic2=equ_Mtx_names[ii,4];
                       ir1=equ_Mtx_names[ii,5];
                       ir2=equ_Mtx_names[ii,6];
                       datdo=equ_mtx[ir1:ir2,ic1:ic2];
		  endif;
                  if rows(datdo)>1 ;
		     mkout=mkout $+ chrs(10|13);
		  endif;

                  blen=0;
                  if dosay2==2 ;
	            for ir (1,rows(datdo),1);
	            for ic (1,cols(datdo),1);
   	         	  foo=g_ftos2(datdo[ir,ic],0);
			  blen=maxc(blen|strlen(foo));
                    endfor;
                    endfor;
		    foo=blen*cols(datdo);
                    if foo>179;
			outwidth 255 ;
		    elseif foo> 100 ;
		        outwidth 180;
	            elseif foo>78 ;
			outwidth 120 ;
		    endif;

                 endif;

                  for ir (1,rows(datdo),1);
		  if blen>0 ; mkout=mkout $+ "       " ; endif;
	          for ic (1,cols(datdo),1);
   	         	  foo=g_ftos2(datdo[ir,ic],0);
			  if blen>0 ; foo=g_padwrd(foo,blen,0); endif;
                          mkout=mkout $+ " " $+ foo; 
                  endfor;
                  if rows(datdo)>1 ;
		     mkout=mkout $+ chrs(10|13);
		  endif;
                  endfor;


                  continue ; 
             endif;		
	     i2=strindx(dosay_it2,"\"",2);
	     if i2==0 ;  mkout=mkout $+  strsect(dosay_it2,2,100000) ; break ; endif; @ close an unclosed " string @
             tt=strsect(dosay_it2,2,i2-2);
             mkout=mkout $+ tt ;
             dosay_it2=strsect(dosay_it2,i2+1,100000);
         endo;
	 screen off; output on ;    $mkout ;	output off; screen on ;
     endif;
  endif;


@ now that equations (including condition) are evaluted, what to do ... @
  if dokeep1==1  ;
       keepdat=keepdat + (isthen./=0);       @ if then is true, be sure to keep @
  endif;
  if dokeep2==1 ;
      keepdat=keepdat+ (iselse./=0) ;        @ if else is true, be sure to keep @
  endif;
  if drop1$=="*";         @ drop obs @
      keepdat=keepdat.*(isthen./=1) ;    @ if true then =1, hence not =0, making keepdat=0 @
  else;
     if drop1$/="" ;
       if sumc(isthen)>0 ;          @ condition EVER true ? @
         if scalmiss(indcv(drop1,equ_drop_vars)); equ_drop_vars=equ_drop_vars|drop1 ;equ_do_drop=1 ; endif;
       endif;
     endif;
  endif;
  if drop2$=="*" ;
      keepdat=keepdat.*(iselse./=1) ;
  else;
     if drop2$/="" ;
       if sumc(iselse) gt 0 ;          @ condition EVER false ? @
         if scalmiss(indcv(drop2,equ_drop_vars)); equ_drop_vars=equ_drop_vars|drop2 ;equ_do_drop=1 ; endif;
       endif;
     endif;
  endif;
 
  if aequc$=="" ;               @ uncondtional never changes an existing variable @
     dat0=dat0~vv1 ;
     nam2=nam2|aname ;
     if scalmiss(indcv(aname,madevars)) ;  madevars=madevars|aname; endif;

  else;                 @ conditional: row values depend on whether THEN or ELSE is used @

    if ADDV1==1 ;
        if addv1_use==0 ;
           d1=miss(zeros(rows(vv1),1),0);
        else;
           d1=dat0[.,addv1_use];
        endif;
        for jz (1,rows(d1),1);
           if isthen[jz]./=0 ; d1[jz]=vv1[jz] ; endif;
        endfor;
     endif;

     if addv2==1 ;
       if aname2$==aname1 ;             @ THEN and ELSE point to same variable @
        for jz (1,rows(d1),1);
          if iselse[jz]./=0 ; d1[jz]=vv2[jz] ; endif;
        endfor;
       else;                            @ then and else point to different variables @
         if addv2_use==0 ;
             d2=miss(zeros(rows(vv2),1),0);
         else;
             d2=dat0[.,addv2_use];
         endif;
         for jz (1,rows(d1),1);
             if iselse[jz]./=0 ; d2[jz]=vv2[jz] ; endif;
         endfor;
       endif;
     endif;  

     if addv1==1 ;
         if addv1_use==0 ;
            dat0=dat0~d1;
            nam2=nam2|aname1 ;
         else;
            dat0[.,addv1_use]=d1 ;
         endif;
         if scalmiss(indcv(aname1,madevars)); madevars=madevars|aname1; endif;
         if addv2==1 and aname1$/=aname2 ;
            if addv2_use==0;
               dat0=dat0~d2;
               nam2=nam2|aname2 ;
            else;
               dat0[.,addv2_use]=d2 ;
            endif;
            if scalmiss(indcv(aname2,madevars));madevars=madevars|aname2; endif;
         endif;
      elseif addv2==1 ;
            if addv2_use==0;
               dat0=dat0~d2;
               nam2=nam2|aname2 ;
            else;
               dat0[.,addv2_use]=d2 ;
            endif;
            if scalmiss(indcv(aname2,madevars)) ; madevars=madevars|aname2; endif;
      endif;
             
  endif;        @ aequc=="" @


  endfor;           @ rows(equ_names) @

@ any observations to drop @
 keepdat=keepdat./=0 ;
 if sumc(keepdat) ne rows(keepdat) ;
     tmp=g_mkindx(keepdat);
     if scalmiss(tmp);
         call grbl2_error("Error: all observations are DROPped ",1);
         stop;
     endif;
     dat=dat0[tmp,.];
  else;
     dat=dat0;
  endif;


endif;              @ new vars @


if sort_do ne 0 ;		@ sort? @
   sid=indcv(sort_vars,nam2) ;
   if ismiss(sid);
     call grbl2_error("ERROR: sorting on a non-existent variable ( " $+ sort_vars,1) ;
      stop;
   endif;
   tmp=seqa(1,1,rows(dat))~dat[.,sid];
   tmp=sortmc(tmp,seqa(2,1,rows(sid))) ;
   tmp=tmp[.,1];
   if sort_do == -1 ; tmp=rev(tmp); endif;
   sortlst=tmp;
@   dat=dat[.,tmp];@
endif;

@ drop any vars @
keeplst=seqa(1,1,cols(dat));
if equ_do_drop ne 0 or equ_do_keep ne 0  ;
   if equ_do_keep==1 ;              @ ONLY keep variables that are in equ_keep_vars and NOT in equ_drop_vars @
       tmp={};
       for j (1,rows(equ_keep_vars),1);
          if scalmiss(indcv(equ_keep_vars[j],nam2)) ;format /ld 1,8;
              call grbl2_error("Error: variable in KEEP (" $+ equ_keep_vars[j] $+ ") not in saved variable list ",1);
              stop;
           endif;
       endfor;
       for j (1,rows(nam2),1);
          if scalmiss(indcv(nam2[j],equ_keep_vars)) ;
              tmp=tmp|nam2[j];
          endif;
       endfor;
       if equ_do_drop ne 0 ; tmp=tmp|equ_drop_vars ; endif;
       equ_drop_vars=unique(tmp,0);
   endif;
   droplist=indcv(equ_drop_vars,nam2);
   if ismiss(droplist);
        call grbl2_error("Error: variable in DROP not in saved variable list ",1);
        stop;
   endif;
   if rows(foo) eq rows(nam2);
        call grbl2_error("Error: all variables are dropped. ",1);
        stop;
   endif;
   foo=seqa(1,1,rows(nam2));
   foo[droplist]=foo[droplist]*0 ; foo=miss(foo,0);
   foo=packr(foo);
   nam2=nam2[foo];
@   dat=dat[.,foo];@
  keeplst=foo;
endif;

closeall ;

if c_outfil$/="" and c_outfil$/="0" ;
  create fhh=^c_outfil with ^nam2,0,4;

  for j(1,rows(dat),50000);
     j2=j+50000-1;
     j2=minc(j2|rows(dat));
     if sort_do ne 0;
        dat2=dat[sortlst[j:j2],keeplst];
     else;
       dat2=dat[j:j2,keeplst];
     endif;
     foo=writer(fhh,dat2);
     if foo ne rows(dat2);
       output on ;
             call grbl2_error("ERROR: unable to write all " $+ g_ftos2(rows(dat),1) $+ " rows of data to "  $+ c_outfil,1);
       output off;
       stop;
    endif;
  endfor;
  a=close(fhh);
endif;

save_mess=0~0;
if i_svlist ne 0 ;                @ save a matrix? @
  i2=1; ?;
  for i2 (2,rows(savevars),1);
    aword1=0 $+ savevars[i2];
    ii=indcv(aword1,equ_mtx_names[.,1]);
    if scalmiss(i1); format /ld 1,8 ;
                   call grbl2_error(" ERROR: can not SAVE undefined variable: " $+ aword1,1) ;
        stop ;
    endif;
    ic1=equ_Mtx_names[ii,3];
    ic2=equ_Mtx_names[ii,4];
    ir1=equ_Mtx_names[ii,5];
    ir2=equ_Mtx_names[ii,6];
    v1=equ_mtx[ir1:ir2,ic1:ic2];
    if equ_mtx_names[ii,2]==1 ; v1=v1'; endif;
    aword2=savefils[i2];
    save ^aword2=v1; format /ld 8,0;
output on ;
    " Saving " $aword1 ;; "(";; rows(v1);; " x ";; cols(v1);; ") to " ;;
output off;
     $aword2 ;
    save_mess=save_mess|(rows(v1)~cols(V1));

  endfor;
endif ;


output on ;
"-----------------------------------------";
if c_outfil$/="" and c_outfil$/="0"  ;
   " New GAUSS dataset created. ";?;
else ;
   " Details on an existing GAUSS dataset. ";?;
endif;

format /ld 8,0;
" :  Input file= " ;;$filename;
"    Number of rows= " ;; nr_orig ;
?;

if c_outfil$/="" and c_outfil$/="0"  ;
  " :  Newly created file= " ;; $c_outfil ;
  "     with ";; format /ld 8,0; rows(dat);;" rows retained.";
else ;
  " : An output dataset was NOT created. ";
endif;

?;" : Variables copied from input dataset: ";?;
if nx>0 ;
  call g_showvec(xname,6,10,"    ");
endif;


if not(ismiss(iid));
  ?;" : ID variables created: ";
  call g_showvec(iid,6,10,"    ");
endif;


if nxdummy >0 ;
  ?;" : Newly created DUMMY variables: ";
  call g_Vudum(dinfo,filename,0,2,xnewinfo.names) ;
endif;


if rows(madevars)>0;
  ?;" : Newly created  variables: ";
  call g_showvec(madevars,6,10,"    ");
endif;


if sort_do ne 0 ;
  ?;
  if sort_do < 0 ;
     " : Descending sort variables: " ;; format /ld 1,8 ; $sort_vars';
  else ;
     " : Sort variables: " ;; format /ld 1,8 ; $sort_vars';
  endif;
endif;



if missrv(id1index,0) ne 0 ;        @ PANELS? @
  ?;
  " : Number of panels (in original dataset) " ;;format /ld 6,0; rows(pids) ;; " (using " ;; 
  if id1name$=="#" ; id1index ; else ; format /ld 1,8; $id1name ; endif;
endif;


if equ_do_drop eq 1 ;
  ?;" : The following variables mentioned above were dropped: ";
  call g_showvec(equ_drop_vars,6,10,"    ");
endif;


if c_outfil$/=""  and c_outfil$/="0" ;
   ?;" : The created file has the following variables:  ";
else;
   ?;" : The following variables would of been written to an output file:  ";
endif;
call g_showvec(nam2,6,10,"    ");


if i_svlist ne 0 ;                @ save some globals? @
?;"The following matrices were saved to .FMT files.";
  for i2 (2,rows(savevars),1);
    aword1=savevars[i2];
    aword2=savefils[i2];
     format /ld 1,8 ; "    " $aword1 ;; "(";;
     format /ld 6,0 ; save_mess[i2,1];; " x ";; save_Mess[i2,2];; ") to " ;;
     format /ld 1,8 ; $aword2 ;
  endfor;
endif ;



if rows(warning_mess)>0;
? ;"The following warnings occured. These usually signal problems with your data.";
   tmp="   "$+ warning_mess ;
   $tmp;
  ?;
endif;



output off ;
closeall ;
return ;



/*******/
@ make global data:

 "X"|"DUMMY"|"ID1",  and ID2 and ID3 may be created

  NR_ORIG - # obs in dataset
  NR_SEL  - # obs after selecting observations
  NR_USE  - # obs after removing obs with missing values
@


proc g_create_makex(foo);


local n ;
local dname,f1;
local t1,t2,i;
local n2,n3,tmp,tmp3,tmp2,npanels;
local oldct,oldid;

closeall ;

open f1=^filename ;

nr_orig=rowsf(f1);  @ global @

nr_sel=nr_orig ;            @ slctlist selection NOT used - new if used instead @

t1=close(f1);

@ read data into global matrices @

used_obs=g_mkdata("X"|"DUMMY"|"ID1",miss(1,1),0,0);

nr_use=rows(x);

if missrv(idindex,0) ne 0 ;
  if idname $== "#" ;
       if idindex==-1 ; got1=rows(x); else ;  got1=idindex ; endif;
       tmp=rows(x) % got1 ;
       if tmp ne 0 ;
        call grbl2_error( "ERROR: rows of the data ( " $+ g_ftos2(rows(x),1) $+ " not even multiple of balanced ID "  $+  g_ftos2(got1,1),1);
        stop;
    endif;
   endif;
endif;


if missrv(id2index,0) ne 0 ;
  if id2name$=="#" ;
    if got1 eq 0;
      call grbl2_error("ERROR: specifying balanced ID2, but not a balanced ID1 ",1);
       stop;
    endif;
    if id2index==-1 ; got2=rows(x); else ;  got2=id2index ; endif;
    if (got1%got2) ne 0 ;
       call grbl2_error("ERROR: balanced ID1 (" $+ g_ftos2(got1,1) $+ ") not even multiple of balanced ID2 (" $+ g_ftos2(got2,0) $+ " ",1);
       stop;
    endif;
  endif;
endif;

if missrv(id3index,0) ne 0 ;
  if id3name $== "#" ;
    if got2 eq 0;
       call grbl2_error("ERROR: specifying balanced ID3, but not a balanced ID2 ",1);
       stop;
    endif;
    if id3index==-1 ; got3=rows(x); else ;  got3=id3index ; endif;
    if (got2%got3) ne 0 ;
       call grbl2_error("ERROR: balanced ID2 (" $+ g_ftos2(got2,1) $+ " not even multiple of balanced ID3 (" $+ g_ftos2(got3,1),1);
       stop;
    endif;
 endif;
endif;


if got1> 0 ;
   npanels=rows(x)/got1;
   id=seqa(1,1,npanels);
   id=id.*.ones(got1,1) ;
endif;

if got2>0 ;
  n2=got1/got2;         @ got2 is number of rows per nest, n2 is # of nest per panel @
  tmp=seqa(1,1,n2);
  tmp2=ones(got2,1);
  tmp3=ones(npanels,1);

  id2=tmp.*.tmp2;
  id2=tmp3.*.id2;
endif;

if got3>0 ;            @ got3 is number of rows per sub-nest @
  n3=got2/got3 ;         @ n3 is # of subnest per nest @
  tmp=seqa(1,1,n3);
  tmp2=ones(got3,1);
  id3=tmp.*.tmp2;
  tmp3=ones(n2,1);
  id3=tmp3.*.id3;

  tmp3=ones(npanels,1);
  id3=tmp3.*.id3;

endif;



retp(nr_use);

endp;


/************************************************************************/
@ Create a new variable, given a text equation @

proc doaneq(equ,varlist2,dat,nthpanel,aequ_mtx,aequ_mtx_names);
local j,a,vv;
local basenode,ith2,tlist;

clearg wordlist,funcs,isneg,val_list,errlist,node_n,node_oper,node_neg,node_left,node_right;

@ process the sentence into several arrays @

{wordlist,errlist,funcs,val_list,isneg}=equ_wordlist(equ);

@ check for valid variable names @
for j (1,rows(wordlist),1);
   tlist=varlist2|aequ_mtx_names[.,1];
   if val_list[j]==2 ;      @ check this variable name @
      if wordlist[j]$=="PI" or wordlist[j]$=="_N" or wordlist[j]$=="_PANEL" ; continue;  endif;
      if strlen(wordlist[j])>8 ;
        errlist=errlist$|("Variable name longer than 8 characters: " $+ wordlist[j]);
      else ;
        if scalmiss(indcv(wordlist[j],tlist));
          errlist=errlist$|("Unknown variable: " $+ wordlist[j] $+ " ... did you include this in the X keyphrase?" );
        endif;
      endif;
   endif;
endfor;

if rows(errlist)>0 ;
  a="ERROR: bad equation. " $+ g_ftos2(rows(errlist),1) $+ " errors found: ";
  a= "   "$|errlist ;
  call grbl2_error(a,1);
   stop;
endif;


@ display pretty version of equation @
@prettyq=equ_dispwordlist(wordlist,funcs,isneg);@


@ build the equation tree, using the globals set above
 Set up tree nodes using 4 "node_" arrays (rather then a structure
 Each node has 4 elements:
    operator  -- either a 0, a math opeartor, or function operator  (node_oper)
                   0 signals "terminal node"
    negative flag                                   (node_neg)
    left pointer   -- points to row in node_ arays  (node_left)
    right pointer -- points to row in node_ arrays, or 0  (node_right)

If var or value node, or if a function node, right pointer is not used
@


node_oper=zeros(5*rows(wordlist),1);
node_neg=node_oper; node_left=node_oper; node_right=node_oper;
node_n=1;        @pointer to last used @

node_oper[1]=".IDENTIT";
node_neg[1]=0; node_right[1]=0;

@ Get_node does the work -- it recursivly creates the equation tree, using
  the globals: wordlist,errlist,funcs,val_list, and isneg@

{basenode,ith2}=get_node(2,varlist2,1,aequ_mtx_names[.,1]);

node_left[1]=basenode;


/*
@ Regenerate equation from the parse tree (uses the node_ globals @
@foo=equ_regenerate(1,varlist2,equ_mtx_names[.,1]);
foo=g_strip(foo);
foo=strsect(foo,2,strlen(foo)-2);
@ send the finalized equation through pretty viewer @
{wordlist,errlist,funcs,val_list,isneg}=equ_wordlist(foo);
prettyq=equ_dispwordlist(wordlist,funcs,isneg);
*/

@ calculate value of the equation tree, given actual values (dat, equ_mtx) @
vv=equ_calc(1,node_oper,node_neg,node_left,node_right,varlist2,dat,nthpanel,aequ_mtx,aequ_mtx_names,0);
retp(vv);

endp;

/*******************************************************/
@ Step 1 for equation parsing: create a "word list" from an equation in
  a sentence @

proc (5)=equ_wordlist(equ);

local opers,opers1,opers2,numchars,funclist,equlen,equ0,wasspace;
local wordlist,errlist,funcs,val_list,isneg,anerr,is_val;
local aword,jj,j,usewords,oldword,neww,isn,activep,wasneg,newf,newn,newv,fflag;
local gotdec,gote,oldchar,equ00,achar;


@ Step 1: take the equation and identify parts
   Parts can be:
      1) numeric values. These contain 0..9,  ', period
         They may also start with - or + .
      2) ( or ) parenthesis.
      3) Names: either functions or variables.
      4) Operators: - + / * \ % ^   (65\6=10, 65%10=5)
    Spaces are ignored (hence, 10 06 is same as 1006 is same as 1,006)

Note: variables must NOT start with a -+.0...9
@

opers="-+*@/\\%^=><][&#:|~'";        @ math and logic opearators @
opers1="-+*@/\\%^=><][(&#:|~'";     @ used in +- assignements @
opers2="-+*@/\\%^=><][(),&#:|~'";   @ used in word seperation @
numchars="0123456789.";  @ note that leading -+ (for minus plus) are not considered part of number @

let funclist=LN EXP SIN COS TAN SIND COSD TAND LOG PDFN CDFN UNIQUE
             GAMMA TRIGAMMA LOGGAMMA DIGAMMA SQRT 
             CORRX VCX
             PDFBVN CDFBVN
             ROW COL INV ELEMENT ROWS
             MIN MAX SD VAR MEAN MEDIAN N NMISS SUM PROD MISS CLASS
             IN IN_RANGE IN_COL FACT LNFACT NOT AVG VAL QUANTILE RNDN RNDU RNDR TRUNC CEIL
             ARCSIN ARCCOS ARCTAN ARCCOST ARCSIND ARCTAND GCDIST ISMISS MISSRV ;

@ remove excess spaces @
equ0="";
wasspace=0;

equ=g_replacestring(equ," OR "," | ",1);   @ search now, before spaces are removed@
equ=g_replacestring(equ," AND "," & ",1);
equ=g_replacestring(equ," GT "," > ",1);
equ=g_replacestring(equ," LT "," < ",1);
equ=g_replacestring(equ," GE "," >= ",1);
equ=g_replacestring(equ," LE "," <= ",1);
equ=g_replacestring(equ," EQ "," == ",1);
equ=g_replacestring(equ," NE "," /= ",1);
equ=g_replacestring(equ,"[","#(",1);
equlen=strlen(equ);

for j (1,equlen,1);
  achar=upper(strsect(equ,j,1));
  if achar$=="[" ;  achar="("; endif;      @ note that { has special meaning @
  if achar$=="]" or achar$=="}" ; achar=")"; endif;
  if achar$=="`" ; continue ; endif;
  if achar$==" " or achar$=="" ;  continue; endif;     @ space  are spacers @
  if vals(achar)<32 ; continue ; endif;  @ crlfs, tabs, etc are skipped @
  equ0=equ0$+achar ;
endfor;

@ substitute ' with '0* if need be @
 equ00="";
 do while 1==1 ;
    i1=strindx(equ0,"'",1);
    if i1==0 ;
        equ00=equ00$+equ0;
        break;
    endif;
    if i1>1;
      equ00=equ00$+(strsect(equ0,1,i1-1));
      equ0=strsect(equ0,i1+1,1000000);
    endif;
    achar=strsect(equ0,1,1); if achar$=="" ; achar=" "  ; endif;
    if achar$=="'"; call grbl2_error("ERROR: '' not allowed ",1); stop ; endif;
    if strindx("-+*/\\%^=><][),&#:|~@ ",achar,1)>0;    @ if this follows a ', then ' DOES NOT imply multiplication@
      equ00=equ00 $+ "'1" ;
    else;                       @ else, does imply multiplication @
      equ00=equ00 $+ "'1*" ;
    endif;
endo;               @ ' fixup @

equ=equ00;

equ=g_replacestring(equ,">=","]",1);
equ=g_replacestring(equ,"/=","~",1);
equ=g_replacestring(equ,"<=","[",1);
equ=g_replacestring(equ,"==","=",1);
equ=g_replacestring(equ,"**","^",1);
equ=g_replacestring(equ,"{","ROW(1,",1);
equ=g_replacestring(equ,"!","!1",1);

equlen=strlen(equ);
wordlist="(";
errlist={};
aword="";  
j=0;
do while 1==1 ;
  j=j+1 ;
  if j>equlen ; break ;endif;
  achar=strsect(equ,j,1);

@ if aword is empty, @
  if aword$=="" ;   @ empty word @
    if strindx(opers2,achar,1)>0;   @ operator? add to wordlist @
       wordlist=wordlist$|achar;    @ note that +- at beginning of number is sorted out later @
       continue;
    endif;
    if strindx(numchars,achar,1)==0;   @ not beginning of a number, start a name @
       aword=aword $+ achar ;
       continue; 
    endif;

@else, beginning of a numbar @
    aword=aword $+ achar ;      @ this is first character of number @
    gotdec=0; gote=0; oldchar="";
    for jj (j+1,equlen,1);          @ get more characters @
      achar=strsect(equ,jj,1);
      j=jj;
      if strindx("0123456789",achar,1)>0;  @ a numeral, add it @
         aword=aword $+ achar ;
         oldchar=achar;
         continue;
      endif;
      if achar$==".";
        if gotdec==1 ;
           call grbl2_error("ERROR: more than one decimal point in a number ",1);
            stop;
        endif;
        gotdec=1;
        aword=aword $+ achar ;
        oldchar=achar;
        continue;
      endif;
      if achar$=="E" ;    @ sci notation ?@ 
         if gote==1 ;
            call grbl2_error("ERROR: more than one  E (scientific notation) in a number ",1);
            stop;
         endif;
         gote=1;
         aword=aword $+ achar ;
         oldchar=achar;
         continue;
      endif;                 
      if achar$=="+" or achar$=="-";    @ +- after E is part of sci notation ?@ 
         if oldchar$=="E" ;
             aword=aword $+ achar ;
             oldchar=achar;
             continue;
         else;              @ this is plus/minuse, and word over! @
             wordlist=wordlist$|aword;
             wordlist=wordlist$|achar ;
             aword="" ;
             break;
         endif;
      endif;

@ if here, number is over -- but must not contain characters @
      if strindx(opers2,achar,1)==0;   @ not an operator? ERROR @
             call grbl2_error("ERROR:a bad character in number (" $+ aword  $+ "," $+ achar $+ " "$+ strsect(equ,1,j),1);
             $equ;
             stop;
      endif;
      wordlist=wordlist$|aword;
      wordlist=wordlist$|achar ;
      aword="";
      break;
    endfor;          @ building a number ( since ) always ends, no need to check for overflow@
    
  else;         @ aword not blank @

   if strindx(opers2,achar,1)>0;   @ operator? add to wordlist @
       wordlist=wordlist$|aword ;
       wordlist=wordlist$|achar;    @ note that +- at beginning of number is sorted out later @
       aword="";
       continue;
    endif;

    aword=aword $+ achar ;      @ else, building a varname or function @

    continue ;
  endif;        @ aword blank @

endo;           @ j @ 

if aword$/="";  wordlist=wordlist$|aword ; endif;
wordlist=wordlist$|")";
wordlist=upper(wordlist);

@ 2) assign - and + to vars/values (if - or + is preceded by ( or an opers)@
usewords=ones(rows(wordlist),1);
isneg=zeros(rows(wordlist),1);
oldword=wordlist[1];
for j (2,rows(wordlist),1);
  aword=wordlist[j];
  if aword$=="-" or aword$=="+";   @ might be a pos/neg rather than plus/minus @
     if strindx(opers1,oldword,1)>0 ;  @ yes, prior word is not var or value @
        if aword$=="-"; isneg[j+1]=1 ; endif;  @ if -, next words is negative @
        usewords[j]=0;      @dont use this one @
     endif;
  endif;
  oldword=aword;
endfor;
neww={};isn={};
for j (1,rows(wordlist),1);        
   if usewords[j] == 1 ; neww=neww$|wordlist[j]; isn=isn|isneg[j]; endif;
endfor;
wordlist=neww;
isneg=isn;

@ 3) Identify functions -- non-opers preceding a (  @
funcs=zeros(rows(wordlist),1);
for j(2,rows(wordlist)-1,1);
  if wordlist[j+1]$=="(";   @ functions always have name then ( @
     if strindx(opers1$+",",wordlist[j],1)==0;  @ this word is NOT a math operator (and not a ( @
        funcs[j]=1;
     endif;
  endif;
endfor;
 
@ 3.1: check for valid functions @
for j (2,rows(wordlist)-1,1);
  if funcs[j]==0; continue; endif;
  aword=wordlist[j];
  if scalmiss(indcv(aword,funclist)) ;  
     anerr="NO SUCH FUNCTION: " $+ aword ; $aword; vals(aword);
     ERRLIST=ERRLIST$|anerr;
  endif;
endfor;
 

@ 4) Idendify values : 1=value, 2=variable, 3=matrix, 0 otherwise  @
val_list=zeros(rows(wordlist),1);

for j (1,rows(wordlist),1);

  achar=strsect(wordlist[j],1,1);
  if strindx(opers2,achar,1)>0;   @ an operator, hence not var or val, skip @
      continue; 
  endif;
  if funcs[j]==1 ;  
     continue;
  endif;            @ function name, cant be value, so skip @

  is_val=2;     @assume is var or matrix @
  if (strindx(numchars,achar,1)>0) ;   @if start with number or period @
      is_val=1;
  endif ; 

  val_list[j]=is_val;      @ if val_list and funcs not 1, must be a variable @
 
  if is_val==2 ;
   if strlen(aword)> 8 ;
      anerr="Variable name too long: " $+ aword ;
      errlist=errlist$|anerr;
   endif;
 endif;

endfor;

@ 5) Check for balanced parenthesis @
activep=0; wasneg=0;
for j (1,rows(wordlist),1);
  aword=wordlist[j];
  if aword$=="("; activep=activep+1; 
  elseif aword$==")" ; activep=activep-1;
  endif;
  if activep<0; wasneg=1 ; endif;
endfor;
if activep gt 0 ;
 errlist=errlist$|"Too many ( parenthesis ";
elseif activep lt 0 ;
 errlist=errlist$|"Too many ) parenthesis ";
elseif wasneg==1 ;
 errlist=errlist$|"Misplaced ) parenthesis ";
endif;

@ 6) make ( into "identity" function (if not preceded by an actual function) @
newf=1; neww=".IDENTIT"; newn=0; newv=0;
fflag="";


for j (2,rows(wordlist),1);
  if fflag$/="" ;      @fflag /="" only if function_name followed by ( @
    newf=newf|1 ;
    newn=newn|isneg[j];
    neww=neww$|fflag;
    newv=newv|0 ;
    fflag="";
    continue;
  endif;
         
  fflag="";  @ assume not a function @
  if funcs[j]==1 ;  fflag=wordlist[j]; continue; endif;  @ it is, skip @

  aword=wordlist[j];
  if aword$=="(";           @ convert simple ( to Idendity function @
    newf=newf|1 ;
    newn=newn|isneg[j];
    neww=neww$|".IDENTIT";
    newv=newv|0 ;
    continue;
  endif;

@ else, keep as is @
   newf=newf|0;
   newn=newn|isneg[j];
   neww=neww$|aword;
   newv=newv|val_list[j] ;
endfor;

isneg=newn;
wordlist=neww;
funcs=newf;
val_list=newv;

retp(wordlist,errlist,funcs,val_list,isneg);
endp;

/*******************************/
@ Display the equation, with parenthesis highlighted @
proc equ_dispwordlist(wordlist,funcs,isneg);
local plevel,parenlevels,ifunc,j,tt,ff,boo,ww,aword,fword,lastlevel;
parenlevels=zeros(rows(wordlist),1);
plevel=0;
for j (1,rows(wordlist),1);
 ifunc=funcs[j];
 aword=wordlist[j];
 if ifunc==1 ;  
    plevel=plevel+1;
    parenlevels[j]=plevel ;
 elseif aword$==")";
    parenlevels[j]=plevel;
    plevel=plevel-1;
 endif;
endfor;

tt="";
ff="";
boo="                                                   ";
@ww="1"|"2"|"3"|"4"|"5"|"6"|"7"|"8"|"9";@
ww="a"|"b"|"c"|"d"|"e"|"f"|"g"|"h"|"i"|"j"|"k"|"l"|"m"|"n"|"o";
lastlevel=1;
for j (1,rows(wordlist),1);

 ifunc=funcs[j];
 aword=wordlist[j];
 if isneg[j]; aword="-"$+aword; endif;
 if aword$==".IDENTIT"; aword="" ; endif;
 if aword$=="["; aword="<="; endif;
 if aword$=="]"; aword=">="; endif;
 if aword$=="~"; aword="/="; endif;
 if aword$=="|"; aword=" OR "; endif;
 if aword$=="&"; aword=" AND "; endif;
 if aword$=="#"; aword=" [ "; endif;

 fword=strsect(boo,1,strlen(aword));
 if funcs[j]==1 ;
    lastlevel=parenlevels[j];
    fword=fword $+ ww[lastlevel] ;
    aword=aword$+"(";
 endif;
 if aword$==")";
    lastlevel=parenlevels[j];
    fword=upper(ww[lastlevel]) ;
    lastlevel=lastlevel-1;
 endif;
 if aword$=="," ;
     fword="|"$+ww[lastlevel];
     aword=aword$+" ";
 endif;
 ff=ff$+ " " $+ fword;
 tt=tt$+" " $+ aword ; 
endfor;
retp(tt$|ff);
endp;


/*********************/
@GET_NODE is the recursive node builder (from parsed info in WORDLIST, etc)
 The IAT is where, in wordlist, to start building a node.
 VARLIST is list of variablenames
 ALEVEL equals 1 if this is a "begin this node" call.

 Nodes have fields: type, negative, left, right
      type:  VAR, VAL, a math operator, or a function name
      isneg:  1 if this is negative (i.e.; 7 * -6 )
      left:  pointer to left subnode, or value or varindex
      right: pointer to right subnode, or not used

 Process:
   i) check alevel -- alevel=1 means "start processing interior of function
                            =0 means "process the next element, which might
                              be a function (possibly IDENTITY)

     a)     if alevel= 0 (process element):
       a.1) if function, creation a function node (note use of the IDENTITY function):
          oper : function name
          neg: from the isneg array for the ith word
          left: node address from recursive call to get_node, which points
              to the "argument tree" for this function (this call sends alevel=1)
          right:  leave blank
       a.2) otherwise, create a val or var node:
          oper: VAR or VAL
          neg : from the isneg array
          left: the value, or pointer to the var (in varlist)
          right: leave blank

    b) alevel=1
        Generate a list of nodes, and the operators between nodes
         These nodes may be values, variables, or functions.
         Also generate a list of operators between the nodes

         If no operators (ie.; LN(100)
           Create a NULL node:
              oper: NULL
              neg: leave blank
              left: the root node of the function (as created in a.1 above)
              right: leave blank

          If one operator: i.e. 995/12
            Create an IDENTITY node, and an additonal operator node
            The identity node points to this operator node
            The operator node has:
                OPER: a math operation
                neg: from isneg array
                left: first element
                2nd element:

          If > one operator: i.e.; 88-2 / exp(2)   (exp(2) is treated as a single element)
             Go through list, creating IDENTITY-OPERATOR pairs of nodes using precedence rules
             When one makes the first such pair, remove the pair from the lists (of nodes and
             operators), and insert the IDENTITY node.
             Iterative processing of this list, as it shrinks, yields a single node
             pointing to the least precedent operatoin

           If multiple areas (comma delimited)
             As above, within each comma delmited area. The comma delimited nodes (which
             might have several operators) are connected into a LIST node.
              OPER: LIST
              NEG: # of areas
              LEFT: node of first area
              RIGHT: pointer to another LIST node
           The other LIST node is the same,
                OPER: LIST
                 NEG: # of areas-1
                 LEFT: node of second area
                 RIGHT: pointer to another LIST node
           This continues until the last area (NEG will be 1, and RIGHT will be 0)

Note: uses globals: wordlist,isneg,funcs,val_list
     (along with node_ globals

Returns pointer to row, in node_ arrays, for this node
 and index to last proceseed word (in wordlist)


----------------------------@

proc (2)= get_node(ith,varlist,alevel,mtxlist);

local anode,aword,ii,ii2,a1node,ii3,aval,ith0,use_node,ith2,ith1;
local localns,localos,j,ok,oldnode,basenode ;

aword=wordlist[ith];
if alevel ne 1 ;        @ not start of a identity or other function @
  use_node=node_n+1;
  node_n=use_node;

  if val_list[ith] eq 1 ;   @ value terminal node @
    node_oper[use_node]=".VAL";
    node_neg[use_node]=isneg[ith];
    aval=stof(aword);
    node_left[use_node]=aval;
    retp(use_node,ith);

  elseif val_list[ith] eq 2 ;     @ var terminal node @
     node_neg[use_node]=isneg[ith];
     node_oper[use_node]=".VAR";
     if aword$=="PI";
         aval=3.14;
     elseif aword$=="_N" ;
         aval=-1 ;          @ -1 means "obs counter" @
      elseif aword$=="_PANEL" ;
         aval=-2 ;          @ -2 means "panel counter" @
     else;                  @ var or matrix @
        aval=indcv(aword,varlist) ; @ pointer into varlist @
        if scalmiss(aval) ;
            aval=indcv(aword,mtxlist) ; @ pointer into matrix-name list @
            if scalmiss(aval);
                call grbl2_error("ERROR: no such variable or matrix " $+ aword,1); stop ;
            endif;
            node_oper[use_node]=".MTX";
        endif;
     endif;
     node_left[use_node]=aval;
     retp(use_node,ith);


 elseif funcs[ith] eq 1 ;     @ function node @
    node_oper[use_node]=aword;
    node_neg[use_node]=isneg[ith];
    {a1node,ii2}=get_Node(ith+1,varlist,1,mtxlist);
    node_left[use_node]=a1node;
    retp(use_node,ii2);  

 else; @ if here ERROR (operator should never be called with alevel ne 1. @
  call grbl2_error("ERROR: alevel ne 1 but operator "  $+aword,1) ;stop;
 endif;   

endif;
 
  
@ otherwise, start of interior of a ( ) (a "function" (often IDENTITY))
 Get list of node, operator, node,... )  -- store in localns and localos @

localns={};
localos=localns;


{anode,ith1}=get_node(ith,varlist,0,mtxlist);      @ note use of 0 @
ith1=ith1+1;
aword=wordlist[ith1];

@ simple case: no operators (say, exp(100) ) @
if aword$==")";  @ end of function, no operator -- use NULL operator @
    use_node=node_n+1; 
    node_n=use_node;
    node_oper[use_node]=".NULL";
    node_neg[use_Node]=0 ;        @ operators are NEVER negative @
    node_left[use_node]=anode;
    @use_node ;; " FCreated =";;  $node_oper[use_node]; @
    retp(use_node,ith1);
endif;

@ else, may be 2 or more elements within ( )  @
localns=localns|anode;
do while 1==1 ; 
  if strindx("-+*/\\%^=[]<>&#:|~,'@",aword,1)==0;    @ is next word NOT an operator. Then error @
    call grbl2_error("ERROR: expected a math operator but found " $+aword,1) ;stop;
  endif;
  localos=localos|aword;
  {anode,ith1}=get_node(ith1+1,varlist,0,mtxlist);      @ note use of 0 @
  localns=localns|anode;
  ith1=ith1+1;
  aword=wordlist[ith1];
  if wordlist[ith1]$==")"; 
     break ;        @ got all the elements within this ( ) @
  endif;
endo;    

@ multiple elements with multiple-1 operators. Assert precedence: ' ^ */%\  -+  ,  
 Do this by creating new IDENTITY nodes @

@ a00 : do  # (matrix element) @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  if localos[j]$/="#" ; continue; endif;
  use_node=node_n+1;
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0; 
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];

@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 




@ a0 : do  ' @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  if localos[j]$/="'" ; continue; endif;
  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0; 
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];

@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 


@ a : do  ^ @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  if localos[j]$/="^" ; continue; endif;
  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0; 
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];

@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 

@ b: do */\% @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  ok=0;
  if    localos[j]$=="*"   or localos[j]$=="/" 
      or localos[j]$=="\\"  or localos[j]$=="%"  ;
     ok=1;
  endif;
  if ok==0;
      continue;
  endif;

  use_node=node_n+1; 

  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0; 
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];
@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 

@ c: do -+ @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  ok=0;
  if localos[j]$=="+" or localos[j]$=="-"  ; ok=1; endif;
  if ok==0;
      continue;
  endif;
  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1;
  node_n=use_node;
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0;
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];
@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 

@ d: do : @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  ok=0;
  if localos[j]$==":"  ; ok=1; endif;
  if ok==0;
      continue;
  endif;
  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0; 
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];
@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo; 




@ e: do <>=[]&#:|~ ampserand ' @
j=0;
do while 1==1 ;
  j=j+1;
  if j>rows(localos); break ; endif;
  ok=0;
  if   localos[j]$=="=" or localos[j]$=="<"  or localos[j]$==">"
    or localos[j]$=="[" 
    or localos[j]$=="]" 
    or localos[j]$=="&" 
    or localos[j]$==":" 
    or localos[j]$=="#" 
    or localos[j]$=="|"  or localos[j]$=="~" or localos[j]$=="'"  or localos[j]$=="@"  ;
      ok=1;
   endif;

  if ok==0;
      continue;
  endif;
  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=".IDENTIT";
  node_neg[use_node]=0;
  
  node_right[use_node]=0;
  node_left[use_node]=use_node+1;

  use_node=node_n+1; 
  node_n=use_node;
  node_oper[use_node]=localos[j];
  node_neg[use_node]=0;
  node_left[use_node]=localns[j];
  node_right[use_node]=localns[j+1];
@ replace j,j+1 in localns with this new node @
  localns[j]=use_node-1;
  localns[j+1]=miss(1,1); localns=packr(localns);
  localos[j]=miss(1,1); localos=packr(localos);
  j=0;
endo;


if not(scalmiss(localos)) ;   @ must be commas @
  use_node=node_n+1; 
  node_n=use_node;
  basenode=use_node;
  node_oper[use_node]=".LIST";
  node_neg[use_node]=rows(localns);
  node_left[use_node]=localns[1];
  for j (2,rows(localns),1);
    oldnode=use_node;    use_node=node_n+1;     node_n=use_node;
    node_right[oldnode]=use_node ;
    node_oper[use_node]=".LIST";
    node_neg[use_node]=1+rows(localns)-j;
    node_left[use_node]=localns[j];
  endfor;
  node_right[use_node]=0;    @ terminal node@

else;               @ not a list @
   basenode=use_node;

endif;  

retp(basenode,ith1);


endp;


@-----------------------------@
@ Regenerate textual equation from the parse tree (the NODE_ arrays @
proc equ_regenerate(anode,varlist,mtxlist);
local lnode,ltt,rnode,rtt,tt,aoper;
aoper="" $+  node_oper[anode];

if aoper$==".VAL";
  tt=g_ftos2(node_left[anode],1);
  if node_neg[anode]==1 ; tt="-" $+ tt; endif;
@"val is: " $tt;@
  retp(tt);
elseif aoper$==".VAR";
  if node_left[anode]==3.14;
    tt="PI";
  elseif node_left[anode]==-1 ;
    tt="_N" ;
  elseif node_left[anode]==-2 ;
    tt="_PANEL" ;
  else;        
     tt="" $+ varlist[node_left[anode]];
  endif;
  if node_neg[anode]==1 ; tt="-" $+ tt; endif;
@"var is: " $tt;@
  retp(tt);

elseif aoper$==".MTX";
  tt="" $+ mtxlist[node_left[anode]];
  if node_neg[anode]==1 ; tt="-" $+ tt; endif;
@"var is: " $tt;@
  retp(tt);

elseif aoper$==".LIST";
  ltt=equ_regenerate(node_left[anode],varlist,mtxlist);
  rnode=node_right[anode];
  do while rnode ne 0 ;
      rtt=equ_regenerate(node_left[rnode],varlist,mtxlist);
      ltt=ltt $+ " , " $+ rtt ;
      rnode=node_right[rnode];
  endo;
  retp(ltt);


elseif aoper$==".NULL";
  rnode=node_left[anode];
@  "get left null equ " rnode " of " anode ;@
  ltt=equ_regenerate(rnode,varlist,mtxlist);
@"null is: " $ltt;@
  retp(ltt);

endif;

if strindx("-+*/\\%^=[]<>&#:|~'@",aoper,1)>0;  

  lnode=node_left[anode];
@"get left equ " lnode " of " anode ;@
  ltt=equ_regenerate(lnode,varlist,mtxlist);

  rnode=node_right[anode];
@"get right equ " rnode  " of " anode;@
  rtt=equ_regenerate(rnode,varlist,mtxlist);
  if aoper$=="[" ; aoper="<=" ;
  elseif aoper$=="]" ; aoper=">="; 
  elseif aoper$=="~" ; aoper="/="; 
  elseif aoper$=="&" ; aoper=" AND ";
  elseif aoper$=="#" ; aoper="[ ";
  elseif aoper$=="|" ; aoper=" OR ";
  endif;
  if aoper$=="'";
     tt=ltt $+ " " $+ aoper $+ " "  ;    @ right node was added by equ processor @
  else ;
     tt=ltt $+ " " $+ aoper $+ " " $+ rtt ;
  endif;
  retp(Tt);
endif;

@ else, function @
if aoper$==".IDENTIT" ; aoper="" ; endif;
tt=aoper $+"( " ;
if node_neg[anode]==1 ; tt="-" $+ tt; endif;
rnode=node_left[anode];
@"get func node " rnode ;; " of " ;;anode;@
rtt=equ_regenerate(rnode,varlist,mtxlist);
tt=tt$+rtt $+ " ) ";
@"func is " $tt;@

retp(tt);
endp;

/*************************************/
@ calulate equation (in the NODE_ arrays), given varlist and dat arrays
  This is recursive, depending on values of oper
   VAL - return value of node_left
   VAR - return node_left column of dat 
   MTX - return matrix of data as
   .NULL - return equ_calc call on node_left
   .IDENTITY - same as NULL
   .LIST - return NxK matrix of values (as drawn from a linked list of nodes)
   operator  - equ_calc(left) + equ_calc(right)
   function_name  - equ_calc to get list of argument, then equ_func on arguments
  listdo : used when processing .LIST type (range of argument to return )
@


proc equ_calc(anode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,listdo);

local v1,v2,vval,opers;
local boper,aoper,lnode,rnode,ineg,donode,specfuncs,plist,plist_all;
local argvals,nxt,ndo,tmp,vv2;
local ithmtx,ir1,ir2,ic1,ic2,ind1,ind2;

opers="-+*/\\%^=><][|&#:~@'";        @ math and logic opearators @
specfuncs="ROW"|"COL"|"ELEMENT"|"VAL"|"ROWS";

aoper=node_oper[anode];
@"aoper " anode;; $aoper ;waitc;@
ineg=1;
if node_neg[anode]==1 ; ineg=-1 ; endif;

@ a simple value @
if aoper$==".VAL" ;        @ value node, just return left_node @
   retp(node_left[anode]*ineg );
endif;

@ the value of a variable (possibly a constant) @
if aoper$==".VAR" ;  @ var node, return node_left column of dat (or PI, or obsnum, or panelnum) @
    if node_left[anode]==3.14;
       v1=ones(rows(dat),1).*pi;
    elseif node_left[anode]==-1 ;    @ obs number @
        v1=seqa(1,1,rows(dat));
    elseif node_left[anode]==-2 ;    @ panel number @
       tmp=maxc(nthpanel|1);
        v1=ones(rows(dat),1).*tmp;
    else;
       v1=dat[.,node_left[anode]];
    endif;
    retp(v1*ineg);
endif;

@ the value of a matrix (possibly a constant) @
if aoper$==".MTX" ;  @ matrix, return node_left matrix @
   ithmtx=node_left[anode];
   ic1=aequ_Mtx_names[ithmtx,3];
   ic2=aequ_Mtx_names[ithmtx,4];

   if nthpanel==0 ;             @ non panel mode, so allow for math-matrices @
      ir1=aequ_Mtx_names[ithmtx,5];
      ir2=aequ_Mtx_names[ithmtx,6];
   else;
      ir1=1 ; ir2=rows(aequ_mtx);       @ all rows in this panel @
   endif;

   v1=aequ_mtx[ir1:ir2,ic1:ic2];
   if aequ_mtx_names[ithmtx,2]==1 ; v1=v1'; endif;
   retp(v1*ineg);
endif;

@ IDENTITY or NULL -- get value of left node @
if aoper$==".IDENTIT" or aoper$==".NULL" ;
    lnode=node_left[anode];
    v1=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,0);
    retp(v1*ineg) ;
endif;

@ LIST of arguments -- return list of values (in columns@
if aoper$==".LIST"; 
   if rows(listdo)==1 ; listdo=listdo|10000000 ; endif;
   ndo=node_neg[anode];  @ number of arguments @
   nxt=anode; argvals={};
   for j (1,minc(ndo|listdo[2]),1);
     if (j ge listdo[1]) ;
       donode=node_left[nxt];
       vval=equ_calc(donode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,0);
       if rows(vval)==1 and rows(dat)>1 ;
         vval=ones(rows(dat),1).*vval ;
       endif;
       argvals=argvals~vval;            @note: error if different sized matrices are the arguments @
     endif;
     nxt=node_right[nxt];    @ nxt=0 at end, but thats okay@
   endfor;
   retp(argvals);
endif;

@ math operator (or logic) operator? @

if strindx(opers,aoper,1)>0;   @ operator. Get left and right side, and then do the operation @

    lnode=node_left[anode];
    v1=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,0);
    rnode=node_right[anode];
    v2=equ_calc(rnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,0);
    if aoper$=="-";
       vval=v1-v2 ;
    elseif aoper$=="+";
       vval=v1+v2 ;
    elseif aoper$=="/" ;
         vval=v1./v2 ;
    elseif aoper$=="\\" ;   @ remove remainder @
         vval=trunc(v1./v2) ;
    elseif aoper$=="%" ;    @ remainder @
         vval=(v1%v2) ;

    elseif aoper$=="@" ;        @ element by element multiplication @
         vval=v1.*v2 ;          @ no dimension checking -- maybe add later @

    elseif aoper$=="*" ;        @ might be a matrix multiplicaiton? @
      if rows(V1)==1 and cols(v1)==1 ;      @ v1 is scalar, v2 can be anything @
            retp(v1.*v2);
      endif;
      if rows(V2)==1 and cols(v2)==1 ;      @ v2 is scalar, v1 can be anytying @
            retp(v1.*v2);
      endif;
      if cols(v1)==1 and cols(v2)==1 ;      @ both  of them are vectors @
         if rows(V1) ne rows(v2);
             call grbl2_error("Error: vector multiplication of different lengths: " $+ g_ftos2(rows(v1),1) $+ " " $+ g_ftos2(rows(v2),1),1);
             stop;
          endif;
          retp(v1.*v2);
      endif;

      if cols(v1) ne rows(V2) ;         @ matrix multiplication @
            call grbl2_error("ERROR: matrix multiplication: " $+ g_ftos2(rows(v1),1) $+  " x " $+ g_ftos2(cols(v1),1) $+ " * "  $+ g_ftos2(rows(v2),1) $+  "x" $+ g_ftos2(cols(v2),1),1) ;
            stop;
      endif;
      vval=v1*v2 ;              @ matrix multiply @

    elseif aoper$=="^";
       vval=v1.^v2 ; 

    elseif aoper$=="!";
       vval=v1! ; 
     
    elseif aoper$=="=" ;   @ logical (0/1) operations @
         vval=(v1.==v2) ;
    elseif aoper$==">" ;
         vval=(v1.>v2) ;
    elseif aoper$=="<" ;
         vval=(v1.<v2) ;
    elseif aoper$=="]" ;
         vval=(v1.>=v2) ;
    elseif aoper$=="~" ;
         vval=(v1./=v2) ;
    elseif aoper$=="[" ;
         vval=(v1.<=v2) ;
    elseif aoper$=="&" ;
       retp((v1.*v2)./=0);

    elseif aoper$==":" ;
       i1=v1[1]; i2=v2[1];
       vval=miss(zeros(rows(dat),1),0);
       tmp=seqa(i1,1,1+i2-i1);
       if rows(tmp)>rows(dat); format /ld 10,0;
	 call grbl2_error("ERROR: too many rows generated by " $+ g_ftos2(i1,1) $+ " : " $+ g_ftos2(i2,1),1) ;
	 stop;
       endif;
       vval[1:rows(tmp)]=tmp;
       retp(vval);

    elseif aoper$=="#" ;		@ extract elements @
       if cols(v2) lt 2 ;
         call grbl2_error("ERROR: did not specify  row(s) and column(s) for matrix element ",1);
         stop;
       endif;

       ind1=packr(v2[.,1]); ind2=packr(v2[.,2]);
       if scalmiss(ind1) or scalmiss(ind2);
         call grbl2_error("ERROR: bad indices (all missing) ",1) ;
	 stop;
       endif;
       if ind1[1]==ind1 ; ind1=ind1[1]; endif;   @ just one element @
       if ind2[1]==ind2 ; ind2=ind2[1]; endif;   @ just one element @

       if rows(v1) lt maxc(ind1) ; format /ld 6,0;
          call grbl2_error("ERROR: row index (" $+ g_ftos2(maxc(ind1),1) $+ " too large (#rows= " $+ g_ftos2(rows(V1)),1);
          stop;
       endif;
       if cols(v1) lt maxc(ind2); format /ld 6,0;
          call grbl2_error("ERROR: row index (" $+ g_ftos2(maxc(ind2),1) $+ " too large (#cols= " $+ g_ftos2(cols(v1),1),1);
          stop;
       endif;
       tmp=v1[ind1,ind2];
       retp(tmp);

    elseif aoper$=="'" ;      @ ' is ALWAYS teamed with 1. no need to check dimensions @ 
       retp(v1');

    elseif aoper$=="|" ;
       retp((abs(v1)+abs(v2))./=0);
    endif;
    retp(vval);
endif;


@ else, a function. First get list of arguments, then apply function  @
lnode=node_left[anode];  

plist_all=0;

@ for ROW, COL, VAL, ROWS : first argument is specification
      ELEMENT : first two are @
if not(scalmiss(indcv(aoper,specfuncs)));   @ need to get param specs? @
   if aoper$=="ELEMENT";  
     plist=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,1|2);
     plist=plist[1,.]';
     argvals=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,3);
   else;
      plist=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,1|1);
@      clearg hoy ; hoy=plist ; call waitc ;  @
      plist_all=plist ;           @ used by rows @
      plist=plist[1,.]';
      argvals=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,2);
   endif;
else ;          @ all arguments are data @
   plist=0;
   argvals=equ_calc(lnode,node_oper,node_neg,node_left,node_right,varlist,dat,nthpanel,aequ_mtx,aequ_mtx_names,1);
endif;
vval=equ_dofunc(aoper,argvals,plist,plist_all);
retp(vval);

endp;

/**********************/
@ compute functions
 LN EXP SIN COS TAN SIND COSD TAND LOG PDFN CDFN  CLASS UNIQUE
 GAMMA LOGGAMMA DIGAMMA TRIGAMMA
 CORRX VCX
 PDFBVN CDFBVN
 MIN MAX SD VAR MEAN MEDIAN N NMISS SUM PROD MISS IN FACT LNFACT AVG VAL QUANTILE
 RNDR RNDU RNDN TRUNC CEIL ROW COL INV ELEMENT ;
@
proc equ_dofunc(afunc,vvals,plist,plist_all);

local tmp,ido,vval,j,atmp,i1,itmp,tmpb,tmpl,tmp1;
local mults,nv,cval,k,rsum,oldtmp,oldval;
local sd1,sd2,rho ;
local rand_values,rr ;


i1=ones(rows(vvals),1);

@ the trivial case -- return first argument @
if afunc$==".IDENTIT" ; retp(vvals[.,1]); endif;

if afunc$=="MISSRV"  ;
  vval=missrv(vvals[.,2],vvals[.,1]);
  retp(vval);
endif;

if afunc$=="MISS" or afunc$=="ISMISS" ;
   tmp1=ones(rows(vvals),1);  @ assume all are missing @
   tmp=seqa(1,1,rows(vvals));
   tmp=packr(vvals~tmp);
   if scalmiss(tmp); retp(tmp1); endif ;      @ ALL rows have missing values @
   if rows(tmp)== rows(vvals) ; retp(tmp1.*0) ;endif ;  @ NO rows have missing values @
   tmp=tmp[.,cols(tmp)];
   tmp1[tmp]=tmp1[tmp].*0;     @ these are the non-mssing rows @
   retp(tmp1);
endif;

if afunc$=="IN" ;
   if cols(vvals)<2 ; retp(zeros(rows(vvals),1)); endif;
   tmp=zeros(rows(vvals),1);
   for j (1,rows(tmp),1);
     atmp=indnv(vvals[j,1],vvals[j,2:cols(vvals)]');
     tmp[j]=missrv(atmp,0);
   endfor;
   retp(tmp);
endif;

if afunc$=="IN_RANGE" ;
   if cols(vvals)<2 ; retp(zeros(rows(vvals),1)); endif;
   rr=seqa(1,1,cols(vvals));
   tmp=zeros(rows(vvals),1);
   for j (1,rows(tmp),1);
     vv2=(vvals[j,2:cols(vvals)]~1e100)' ;
     tmp[j]=subscat(vvals[j,1],vv2,rr);
   endfor;
   retp(tmp);
endif;

if afunc$=="IN_COL" ;
  vval=zeros(rows(vvals),1) ;
  tmp=vvals[1,1];
  vval[1]=missrv(indnv(tmp,vvals[.,2]),0);
  tmp1=vvals[.,1].==tmp;
  if sumc(tmp1)==rows(vvals);   @ if first arg is a scalar @
      vval=vval[1];
     retp(vval);
  endif;

  oldtmp=tmp; oldval=vval[1];  @ the "Cache" @
  for j (2,rows(vvals),1);
     tmp=vvals[j,1];
     if tmp==oldtmp ; vval[j]=oldval ;
     else;
       vval[j]=missrv(indnv(tmp,vvals[.,2]),0);
       oldtmp=tmp; oldval=vval[j];
     endif;
  endfor;
  retp(vval);
endif;


if afunc$=="GCDIST" ;
   if cols(vvals)<4 ;
       call grbl2_error("ERROR: GCDIST requires 4 columns (Lat1, Long1, Lat2, Long2)",1);
       stop;
   endif;
   tmp=gcdistance(vvals[.,1],vvals[.,2],vvals[.,3],vvals[.,4]);
   retp(Tmp);
endif;




if afunc$=="INV" ;
   if cols(vvals) ne rows(vvals);
      call grbl2_error("ERROR: inv of non-square matrix: " $+ g_ftos2(rows(vvals),1) $+  "x"  $+ g_ftos2(cols(vvals),1),1);
       stop;
   endif;
   trap 1 ;
   tmp=inv(vvals);
   trap 0;
   if scalerr(Tmp) ne 0 ;
     warning_mess=warning_mess$|"Singular matrix (will use missings) ";
      tmp=miss(zeros(rows(vvals),cols(vvals)),0);
   endif;
   retp(tmp);
endif;

if afunc$=="ROW" ;
   j=plist[1,1];
   if j== 0;
     retp(rows(vvals));
   endif;
   if j>rows(vvals);
      call grbl2_error("ERROR: no such row " $+ g_ftos2(j,1) $+ " (# rows= "  $+ g_ftos2(rows(vvals),1),1);
       stop;
   endif;
   retp(vvals[j,.]');    @ return as Kx1 vectr @
endif;

if afunc$=="ROWS" ;

   j=packr(plist_all)  ;   @ remove missing @
   if sumc(j.==0)==rows(j); retp(rows(vvals)) ; endif;  @ special case: 0s @

   if scalmiss(j); retp(miss(1,1)); endif;     @ no good rows @
   if maxc(j)>rows(vvals) ;
       call grbl2_error("ERROR: a ROWS index is too large:  " $+ g_ftos2(maxc(j),1) $+ " (# rows= "  $+ g_ftos2(rows(vvals),1),1);
       stop;
   endif;
   if minc(j)<1 ;
       call grbl2_error("ERROR: a ROWS index  is < 1 :  " $+ g_ftos2(minc(j),1) );
       stop;
   endif;
   retp(vvals[j,.]);    @ return as Rows(j) x cols(vvals) (often Nx1) matrix @
endif;

if afunc$=="COL" ;

   j=plist[1,1];
   if j== 0;
     retp(cols(vvals));           @ since first col is argument,drop it @
   endif;
   if j>(cols(vvals));
      call grbl2_error("ERROR: no such col " $+ g_ftos2(j,1) $+  " (# cols= "  $+ g_ftos2(cols(vvals),1),1);
       stop;
   endif;
   retp(vvals[.,j]);
endif;


if afunc$=="ELEMENT" ;
   j=plist[1,1];
   k=plist[2,1];
   if j>rows(vvals) or k>cols(vvals) or j<1 or k<1;
      call grbl2_error("ERROR: no such col or row " $+ g_ftos2(j,1) $+ " " $+ g_ftos2(k,1),1);
       stop;
   endif;
   retp(vvals[j,k]);
endif;



if afunc$=="CLASS" ;
   cval=zeros(rows(vvals),cols(vvals));
   nv=zeros(cols(vvals),1);
   for k (1,cols(vvals),1);
     tmp=unique(vvals[.,k],1);
     nv[k]=rows(tmp);
     cval[1:nv[k],k]=tmp;
   endfor;
@trivial case: one column has no duplicates @
  if maxc(nv) eq rows(vvals) ;
      retp(seqa(1,1,rows(vvals)));
  endif;
  mults=ones(rows(nv),1);
  for j (rows(nv)-1,1,-1);
     tmp=mults[j+1]*nv[j+1]; ;
     mults[j]=tmp;
  endfor;

  vval=zeros(rows(vvals),1);
  for j  (1,rows(vvals),1);
     rsum=0;
     for k (cols(vvals),1,-1);
       tmp=indnv(vvals[j,k],cval[1:nv[k],k]);
       rsum=rsum+(tmp*mults[k]);
     endfor;
     vval[j]=rsum;
  endfor;
  retp(vval);
endif;

if afunc$=="EXP" ;
  retp(exp(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="FACT" ;
  retp(vvals!);  @ ignore other arguments @
endif;

if afunc$=="LNFACT" ;
  retp(lnfact(vvals));  @ ignore other arguments @
endif;

if afunc$=="NOT";
   retp(vvals.==0);
endif;

if afunc$=="LN";
  if cols(vvals)>1 ;
     vvals=substute(vvals,vvals .le 0,vvals[.,2]);
  endif;
  retp(ln(vvals[.,1]));
endif;

if afunc$=="GAMMA";
  tmp=gamma(vvals[.,1]);
  retp(tmp);
endif;

if afunc$=="DIGAMMA";
  retp(digamma(vvals[.,1]));
endif;

if afunc$=="TRIGAMMA";
  retp(trigamma(vvals[.,1]));
endif;

if afunc$=="LOGGAMMA";
  tmp=vvals.>100 ;
  tmpb=g_mkindx(tmp);
  tmpl=g_mkindx(tmp.==0);
  tmp=zeros(rows(vvals),1);
  if not(scalmiss(tmpb)) ;
     tmp[tmpb]=lnfact(vvals[tmpb,.]-1);
  endif;
  if not(scalmiss(tmpl)) ;
     tmp[tmpl]=ln(gamma(vvals[tmpl,.]));
  endif;
  retp(tmp);
endif;

@ lnfact(3)=(ln(gamma(4)));@

 if afunc$=="LOG";
  if cols(vvals)>1 ;
     vvals=substute(vvals,vvals .le 0,vvals[.,2]);
  endif;
  retp(log(vvals[.,1]));
endif;

if afunc$=="SQRT";
  tmp=sqrt(vvals[.,1]);
  itmp=vvals[.,1].<0;
  tmp=substute(tmp,itmp,-1);
  tmp=miss(tmp,-1);
  retp(tmp);
endif;


if afunc$=="RNDU";
     if cols(vvals) ge 2 ;
        rand_values=vvals[.,1:2];
     else;
        rand_values=zeros(rows(vvals),1)~ones(rows(vvals),1) ;   @ min , max @
     endif;
     tmp=rndu(rows(vvals),1);
     tmp=tmp.*(rand_values[.,2]-rand_values[.,1]) + rand_values[.,1] ;
     retp(tmp);
endif;

if afunc$=="RNDN";
     if cols(vvals) ge 2 ;
        rand_values=vvals[.,1:2];
     else;
        rand_values=zeros(rows(vvals),1)~ones(rows(vvals),1) ;    @ mean sd @
     endif;
     tmp=rndn(rows(dat),1);
     tmp=(tmp.*rand_values[.,2]) + rand_values[.,1] ;
     retp(tmp);
endif;

if afunc$=="RNDR" ;
     tmp=rndu(rows(vvals),1).*cols(vvals);
     itmp=ceil(tmp);
     vval=zeros(rows(vvals),1);
     for j (1,rows(vval),1);
        vval[j]=vvals[j,itmp[j]];
     endfor;
     retp(vval);
endif;

if afunc$=="TRUNC";
  retp(trunc(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="CEIL";
  retp(ceil(vvals[.,1]));  @ ignore other arguments @
endif;

@ trig @

if afunc$=="TAN";
  retp(tan(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="TAND";
  TMP=2*PI/360;
  retp(tan(TMP.*vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="SIN";
  retp(sin(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="SIND";
  TMP=2*PI/360;
  retp(TMP.*sin(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="COS";
  retp(cos(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="COSD";
  TMP=2*PI/360;
  retp(TMP.*cos(vvals[.,1]));  @ ignore other arguments @
endif;


@ -- arc trig @

if afunc$=="ARCTAN";
  retp(atan(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="ARCTAND";
  TMP=180/PI;
  retp(tmp.*atan(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="ARCSIN";
  retp(tmp.*arcsin(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="ARCSIND";
  TMP=180/PI;
  retp(TMP.*arcsin(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="ARCCOS";
  retp(arccos(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="ARCCOSD";
  TMP=180/PI;
  retp(TMP.*arccos(vvals[.,1]));  @ ignore other arguments @
endif;



if afunc$=="PDFN";
  retp(pdfn(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="CDFN";
  retp(cdfn(vvals[.,1]));  @ ignore other arguments @
endif;

if afunc$=="CDFBVN";
  if cols(vvals)<2 ;            @ requires at least two rows @
     call grbl2_error("ERROR: must specify at least two columns with CDFBVN ",1) ;
     stop;
  endif;

  if cols(vvals)<2 ;            @ requires at least two rows @
     call grbl2_error("ERROR: must specify at least two columns with PDFBVN ",1) ;
     stop;
  endif;
  sd1=1; sd2=1; rho=0;
  if cols(vvals) ge 3 ;
     sd1=vvals[.,3];
  endif;
  if cols(vvals) ge 4;
    sd2=vvals[.,4];
  endif;
  if cols(vvals) ge 5 ;
    rho=vvals[.,5];
  endif;

  tmp=cdfbvn(vvals[.,1]./sd1,vvals[.,2]./sd2,rho);
  retp(tmp);
endif;


if afunc$=="PDFBVN";
  if cols(vvals)<2 ;            @ requires at least two rows @
     call grbl2_error("ERROR: must specify at least two columns with PDFBVN ",1) ;
     stop;
  endif;
  sd1=1; sd2=1; rho=0;
  if cols(vvals) ge 3 ;
     sd1=vvals[.,3];
  endif;
  if cols(vvals) ge 4;
    sd2=vvals[.,4];
  endif;
  if cols(vvals) ge 5 ;
    rho=vvals[.,5];
  endif;
  tmp=a_pdfbvn(vvals[.,1],vvals[.,2],sd1,sd2,rho);
  retp(tmp);
endif;

@ Stat FUNCTIONS 
  if multiple columns, compute statistics across columns (each row is an observation, each 
  column a seperate variable)
  if one column, compute across rows (across the observations).

  ..@ 

if afunc$=="MIN";
  if cols(vvals)==1 ;
      retp(i1.*minc(vvals));
  else;
     retp(minc(vvals'));
  endif;
endif;

if afunc$=="MAX";
  if cols(vvals)==1 ;
      retp(i1.*maxc(vvals));
  else;
     retp(maxc(vvals'));
  endif;
endif;

if afunc$=="MEAN" or afunc$=="AVG" ;
  if cols(vvals)==1 ;
      retp(i1.*meanc(vvals));
  else;
     retp(meanc(vvals'));
  endif;
endif;

if afunc$=="SUM";
  if cols(vvals)==1 ;
      retp(i1.*sumc(vvals));
  else;
     retp(sumc(vvals'));
  endif;
endif;

if afunc$=="PROD";
  if cols(vvals)==1 ;
      retp(i1.*prodc(vvals));
  else;
     retp(prodc(vvals'));
  endif;
endif;

if afunc$=="N";
  if cols(vvals)==1 ;
      retp(i1.*rows(vvals));
  else;
     retp(cols(vvals));
  endif;
endif;


if afunc$=="UNIQUE";
  if cols(vvals)==1 ;
      tmp=rows(unique(vvals,1));
      retp(i1.*tmp);
  else;
     vval=ones(rows(vvals),1);
     for j (1,rows(vvals),1);
       tmp=unique(vvals[j,.]',1);
       vval[j]=rows(tmp);
     endfor;
     retp(vval);
  endif;
endif;


if afunc$=="NMISS";
  if cols(vvals)==1 ;
      tmp=rows(vvals)-rows(packr(vvals));
      retp(i1.*tmp);
  else;
     vvalS=abs(vvalS)+2;  
     vvalS=missrv(vvalS,0);
     vval=sumc((vvalS.==0)');
     retp(vval);
  endif;
endif;

if afunc$=="QUANTILE";
  if cols(vvals)==1 ;
    call grbl2_error("ERROR: first argument of QUANTILE must be the level (0.0 to 1.0) ",1);
    stop;
  endif;

  if cols(vvals==2);
     ido=ceil(rows(VVALS).*vvals[1,1]);
     ido=minc(maxc(ido|1)|rows(vvals));
     vval=vvals[.,2];
     tmp=sortc(vval,1);
     vval=vval|tmp[ido];
     retp(vval.*i1);
  else;
     ido=ceil(rows(VVALS).*vvals[1,1]);
     ido=minc(maxc(ido|1)|cols(vvals));
     vval={};
     for j (2,rows(VVALS),1);
       tmp=VVALS[j,.]';
       tmp=sortc(tmp,1);
       vval=vval|tmp[ido];
     endfor;
     retp(vval);
  endif;
endif;


if afunc$=="VAL";
  if cols(vvals)==1 ;
     call grbl2_error("ERROR: first argument of VAL must be an observation number  ",1);
    stop;
  endif;

  if cols(vvals==2);
     ido=vvals[1,1];
     if scalmiss(ido)==1; ido=rows(vvals) ; endif;
     ido=minc(maxc(ido|1)|rows(vvals));
     vval=vvals[.,2];
     tmp=sortc(vval,1);
     vval=vval|tmp[ido];
     retp(vval.*i1);
  else;
     ido=vvals[1,1];
     if scalmiss(ido)==1; ido=cols(vvals) ; endif;
     ido=minc(maxc(ido|1)|cols(vvals));
     vval={};
     for j (2,rows(VVALS),1);
       tmp=VVALS[j,.]';
       tmp=sortc(tmp,1);
       vval=vval|tmp[ido];
     endfor;
     retp(vval);
  endif;
endif;

if afunc$=="MEDIAN";

  if cols(vvals)==1 ;
     ido=ceil(rows(VVALS)/2);
     tmp=sortc(vvals,1);
     vval=vval|tmp[ido];
     retp(vval.*i1);
  else;
     ido=ceil(cols(VVALS)/2);
     vval={};
     for j (1,rows(VVALS),1);
       tmp=VVALS[j,.]';
       tmp=sortc(tmp,1);
       vval=vval|tmp[ido];
     endfor;
     retp(vval);
  endif;
endif;

if afunc$=="SD";
  if cols(vvals)==1 ;
      retp(i1.*stdc(vvals));
  else;
    retp(stdc(vvals'));
  endif;
endif;

if afunc$=="VAR";
  if cols(vvals)==1 ;
      retp(i1.*(stdc(vvals).^2) );
  else;
     retp(stdc(vvals').^2);
  endif;
endif;

if afunc$=="VCX";
  retp(vcx(vvals));
endif;


if afunc$=="CORRX";
  retp(corrx(vvals));
endif;


 call grbl2_error("ERROR: unknown function " $+afunc,1) ;
stop;

 
endp;  




/***************************************/
@ pdf for bivariate normal  with 0 means, but not necessarily sd1=sd2=1
  x1 and x2 are upper bounds of integration
  all terms may be Nx1 
  Typically, x1 and x2 are Nx1
  sd1,sd2, and rho are 1x1, or Nx1
@

proc a_pdfbvn(x1,x2,sd1,sd2,rho);
local t1,t2,z1,z2,z3,z,t2d,p;
local p2;

t1=2.*pi.*sd1.*sd2.*sqrt(1-(rho.*rho));

z1=(x1.^2)./(sd1.*sd1);
z2=2.*rho.*(x1.*x2);
z2=z2./(sd1.*sd2);
z3=(x2.^2)/(sd2.*sd2);
z=(z1+z3-z2) ;
t2d=2.*(1-(rho.*rho));
t2=exp(-z./t2d);

p=t2./t1 ;
retp(p);
endp;



/***************************************/
@ Procedure to process a CREATE ... RUN code block. Sets globals @

proc g_equation_parse(acommand);

 local comlen,aword,foo,slist,ietmp,foo1,a;
 local idv,aaa,mmm,mm;
 local modelx,arf,iee ;
 local exactxyz;
 local tmpc,aword1,aword2,aword3,aword4,t1,t2,t3,t4,nn,ndumvar;
 local bvals;
 local alog,ALOGS,ALOG2,ALOG1, rr,cc,tmp,rr1,rr2;
 local probchi,adiff,df,aww;
 local vtype0,vtypei1,bb,qtext,eek ;
 local skipto,i1 ;
 local fgoo ;
 local defword,defword2,addrho,jjmp;
 local term1,term2,condn,isthen,ii1,iik;
 local ien,iff;
 local nthval,topers,acomm2,ac ;
 local i2,i3,i1a,tmpr,tt;
 local acommandt,oldvar,newvar;
 local tmp1,tmp2,tmp3,tmp4,ifoo1,cond1,tmp5,ie0;
  local aname,aequ,acond,atodo,atodo2,atodoelse,ii2,ismtx;
 local com1a,com1b,commands,keywurd,ie;

 exactxyz=1; 


@ these globals are used by the G_create subroutine. They are DECLARED in grbl2_create.dec @
equ_names={}; equ_equs={}; equ_conds={};
equ_names2=equ_names ; equ_equs2=equ_equs; equ_ismtx=equ_conds;
equ_drop_vars={}; equ_keep_vars={}; 
equ_do_drop=0; equ_do_keep=0 ;

@ get rid of embedded crlfs  @
tmp=vals(acommand);
tmp=substute(tmp,(tmp.==10) + (tmp.==13),32);
commands=chrs(tmp);
com1a=1;

do while 1==1 	;		@ do forever @

 if commands$=="" ; retp(1); endif;

 com1b=strindx(commands,";",1);
 if com1b==0 ;           retp(1);         endif;

 acommand=upper(strsect(commands,com1a,com1b-1)) ; @drop end ; @
 commands=strsect(commands,com1b+1,1000000);


 {KEYWURD,ie}=g_getwrd(acommand,1) ;

 keywurd=g_strip(upper(keywurd)) ;

@ Depending on KEYWURD, do different things @

if keywurd$=="CREATE";
   continue ;

elseif keywurd$=="RUN" ;		@ done @
  retp(1);

elseif keywurd$=="X" or keywurd$=="DUMMY" or keywurd$=="INPUT" or keywurd$=="FILE" or keywurd$=="OUTPUT" ;

   call g_batch1(acommand,0);

elseif keywurd$=="DROP" ;
   {foo,equ_drop_vars,tmp,ie }=g_Nopts("" $+ keywurd $+ " variables",
                                acommand,ie,exactxyz,"*","0"|"NO");
   equ_do_drop=0;
   if tmp gt 0 ;
      equ_do_drop=1 ;
    endif;

elseif keywurd$=="KEEP" ;
    equ_do_keep=0 ;

    {foo,equ_keep_vars,tmp,ie}=g_Nopts("" $+ keywurd $+ " variables",
                                acommand,ie,exactxyz,"*","0"|"NO"|"*");

   if equ_keep_vars$=="*" ;
       equ_do_drop=0; equ_do_keep=0; equ_drop_vars={}; equ_keep_vars={};
   elseif tmp gt 0 ;
      equ_do_keep=1;
      equ_do_drop=0; equ_drop_vars={};  @ keep overrides drop @
   endif;

elseif keywurd$=="SORT" ;
    sort_do=1;
    tmp=strsect(acommand,ie+1,1111111);

    if tmp$=="" ;
        "ERROR: No Sort variable  specified ";
        sort_do=0;
        call g_pause(10,0);
        retp(0);              @ return to user input @
    endif;
    foo=vals(tmp); foo2=substute(foo,foo.==9,32);
    tmp2=upper(chrs(foo2));
    tmp={};
    do while 1==1 ;
           if tmp2$=="" ; break ;endif;
           {aword,tmp2}=token(tmp2) ;
           aword=g_strip(upper(aword));
           if aword$=="-D" ;
               sort_do=-1;
           else;
               tmp=tmp|aword;
           endif;
     endo;
     if rows(tmp) > 3 ;
           " ERROR: Too many Sort variables specified ";
           sort_do=0;
           call g_pause(10,0);
           retp(0);              @ return to user input @
     endif;
     sort_vars=tmp ;

elseif (keywurd$=="BY") ;

          {id1index,id1name,foo,ie}=g_1opt("BY variable ",
                                  acommand,ie,exactxyz,varlist,1,"NO") ;
            if foo==3 ;
               id1index=0 ; id1name="" ;
          elseif foo==2 ;
             id1index=id1name ; id1name="#" ;
          elseif scalmiss(id1index); 
		     "ERROR: No such BY variable: " ;; $acommand ;
              retp(0);
           endif;
 
elseif (keywurd$=="OUTDATA")  ;
         {make_OUT,ie}=g_getwrd(acommand,ie+1);
         c_outfil=make_out ;
         if c_outfil$=="0";
              "Note: results will NOT be written to a new datafile ";
         else;
            if strindx(upper(make_out),".",1) == 0 ;
                make_out=make_out $+ ".DAT" ;
             endif;
             make_out=g_fspec2("A",make_out,work_dir);
             if not(filesa(make_out) $== "") ;
                       "Warning: overwriting dataset: " $make_out ;
             else;
                        " Creating gauss dataset: "  make_out;
             endif;
         endif;

elseif keywurd$=="NEW" ;
  {aa,foo}=g_getwrd(acommand,ie+1);
   aa=upper(g_strip(aa));
   if aa$/="IF";             @ NOT CONDITIONAL @
         ismtx=0;
         if aa$=="MATRIX";             @ create a matrix NOT CONDITIONAL @
            ismtx=1;
            ie=foo;         @ skip over the MATRIX word @
         endif;
         ii=strindx(acommand,"=",ie+1);
         if ii==0 ;
	            "ERROR: no = in a NEW command:  ";;$acommand ;   retp(0);

         endif;
         aname=g_strip(strsect(acommand,ie+1,(ii-(ie+1))));
         if strlen(aname)>8 ;
 	            "ERROR: new variable names must be 8 or fewer characters : " ;;$acommand ;
         endif;
         aequ=g_strip(strsect(acommand,ii+1,1000000));
         equ_conds=equ_conds$|"" ;
         equ_names=equ_names$|aname;
         equ_equs=equ_equs$|aequ ;
         equ_names2=equ_names2$|"" ;
         equ_equs2=equ_equs2$|"" ;
         equ_ismtx=equ_ismtx|ismtx ;

    else;                          @ CONDITIONAL @

        ie=foo;
        ii=strindx(acommand," THEN ",ie);
        if ii==0 ;
	       "ERROR: no THEN in an IF command:  ";;$acommand ; retp(0);
         endif;

         acond=strsect(acommand,ie+1,(ii-(ie+1)));
         equ_conds=equ_conds $| acond ;
         equ_ismtx=equ_ismtx|0;

         ii2=strindx(acommand," ELSE ",ie+1);
         if ii2==0 ;
           atodo=strsect(acommand,ii+6,10000);
           atodoelse="" ;
         else;
           atodo=strsect(acommand,ii+6,ii2-(II+6));
           atodoelse=strsect(acommand,ii2+6,10000);
         endif;

@ the THEN equation @
         {aword,foo}=upper(token(atodo));
         if aword$=="DROP"  ;   @ next word is the variable @
             equ_equs=equ_equs $|".DROP"  ;  @ add . just to be safe @
             if foo$/="";
                {aname,foo}=upper(token(foo));
                if strlen(aname)>8 ;
       	            "ERROR: DROP variable name (in IF THEN) must be 8 or fewer characters : " ;;$acommand ;
                    retp(0);
                endif;
                equ_names=equ_names $| aname;
             else;
                equ_names=equ_names $|"*";       @ hing means drop obs @
             endif;

         elseif aword$=="KEEP" ;    @ keep only used for observations @
             equ_equs=equ_equs $|".KEEP"  ;  @ add . just to be safe @
             equ_names=equ_names $|"*";       @ hing means drop obs @

         else;                         @ must be a new variable @
            ii=strindx(atodo,"=",1);
            if ii==0 ;
	            "ERROR: no = in an IF ... then x=..  command:  ";;$acommand ; retp(0);
            endif;
            aname=g_strip(strsect(atodo,1,ii-1));
            if strlen(aname)>8 ;
 	            "ERROR: new variable name (in IF THEN) must be 8 or fewer characters : " ;;$acommand ;
                 retp(0);
             endif;
             equ_names=equ_names $| aname;
             aequ=g_strip(strsect(atodo,II+1,1000000));
             equ_equs=equ_equs $|aequ  ;
         endif;

@ the ELSE equation @
         if atodoelse$=="" ;       @No else ? @
              equ_names2=equ_names2$|"";
              equ_equs2=equ_equs2 $| "" ;
         else;
           atodo=atodoelse;
           {aword,foo}=token(atodo);aword=upper(g_strip(aword));
           if aword$=="DROP" ;   @ next word is the variable @
             equ_equs2=equ_equs2 $|".DROP"  ;  @ add . just to be safe @
             if foo$/="";
                {aname,foo}=upper(token(foo));
                if strlen(aname)>8 ;
       	            "ERROR: DROP variable name (in IF .. ELSE ) must be 8 or fewer characters : " ;;$acommand ;
                    retp(0);
                endif;
                equ_names2=equ_names2 $| aname;
             else;
                equ_names2=equ_names2 $|"*";       @ nothing means drop obs @
             endif;

         elseif aword$=="KEEP" ;    @ keep only used for observations @
             equ_equs2=equ_equs2 $|".KEEP"  ;  @ add . just to be safe @
             equ_names2=equ_names2 $|"*";       @ hing means drop obs @

          else;                         @ must be a new variable @
            ii=strindx(atodo,"=",1);
            if ii==0 ;
   	            "ERROR: no = in an IF ... then x=..  command:  ";;$acommand ; retp(0);
            endif;
            aname=upper(g_strip(strsect(atodo,1,ii-1)));
            if strlen(aname)>8 ;
 	           "ERROR: new variable name (in IF THEN) must be 8 or fewer characters : " ;;$acommand ;
                retp(0);
            endif;
            equ_names2=equ_names2 $| aname;
            aequ=g_strip(strsect(atodo,II+1,1000000));
            equ_equs2=equ_equs2 $|aequ  ;
          endif;
        endif;

   endif;               @ NEW @


else;                 @ unknown keywurd @
               " WARNING: Unknown keyword: " $keywurd;
               retp(0);

endif;

endo;		@ processing keyphrases @

endp;



/*******************************************/
@ Compute great circle (as the crow flys) distance between two lat/long pairs @
proc gcdistance(lat1,long1,lat2,long2);
local orlatr,dslatr,lndifr,tmp,gcdist,arcdst;
local alat1,along1,alat2,along2;
local res,j ;

res=zeros(rows(lat1),1);
 long1=abs(long1);
 long2=abs(long2);
  lat1=abs(lat1);
  lat2=abs(lat2);

for j (1,rows(lat1),1);
   along1=long1[j];
   along2=long2[j];
   alat1=lat1[j];
   alat2=lat2[j];
   if ismiss(along1|along2|alat1|alat2) ; res[j]=miss(1,1) ; continue; endif;
  if (ABS(alat1-alat2) Lt 0.001  and ABS(along1-along2) LT 0.001) ;  continue;   endif;

   ORLATR = alat1 * (pi./180);
   DSLATR = alat2 * (pi./180) ;
   LNDIFR = (ABS(along1-along2))* (pi./180) ;

@     COMPUTE GREAT CIRCLE DISTANCE. @
    tmp = ((SIN(ORLATR) .* SIN(DSLATR)) +
            (COS(ORLATR) .* COS(DSLATR) .* COS(LNDIFR))) ;

      if (abs(tmp) gt 1.0)  ;
               gcdist=0 ;
          else ;
             arcdst=arccos(tmp) ;
             GCDIST = 3953.0 .* ARCDST ;
       endif ;
       res[j]=gcdist;
endfor;

retp(res);
endp;
