/*****  Contains 'math' routines of various sorts ****
**
**  proc dolog(x,expln0) ;
**   Return ln of all values in matrix x. For all values <= 0,  use
     ln(expln0). If expln0=miss(1,1), return missing for values <=0.
**
** proc inttobin(num)
**   convert a decimal integer to a binary number, return in an array
**      eg:  vv=inttobin(7) 
**      vv would be (4x1): 0 1 1 1 
**
** proc appnd_it(mtx,row)
**       Append row to matrix. Matrix is NxK1, row is 1xK2 or K2x1.
**        Differs from | because!!
**            If mtx missing, K1 ne K2 , then set mtx=row
**     Use aa=appnd_it(mtx',row)' to append a column !
**

Proc pdf_Mvn(x,mu,vc):

  vals=pdf_mvn(x,mu,vc)

       x: NxK vector of 1xK vectors to compute pdf's of
      mu: Kx1 vector of "means"
      vc: KxK covariance matrix
   
    vals: N x 1 vector of probabilities

 if Mu=scalar, use a Kx1 vector (each element = to mu)
 if VC=scalar, use a KxK diagonal matrix (each diagonal element=vc)

  If you'll be repeatedly calling pdf_mvn with the same values of mu and vc, you
  can use the following faster version:

      vals=pdf_mvn_q(x,mu,invvc,tdenom);

  where invvc and tdenom are computed using;

      {invvc,tdenom}=pdf_mvn_factors(vc);


proc (2)=rndn2(u,cv,n);
     Generates a matrix of multivariate vectors.
     U is the mean, CV is the covariance, and N is the # of vectors.
     Each row of the resulting NxK (K=rows(u)) matrix is a N(U,CV) vector.

   eg;   {Xn,cvtmp}=rndn2(U3,CV3,N)  ;
             with rows(u3)=rows(cv3)=3
       will return a 10x3 matrix.  Each row is a random vector
       distributed as N(U3,CV3).

       CVTMP is a matrix with property: A*A'=CV.  It can be used to
       speed up evaluation of rndn2.

       Note: If cv3[1,1] were < 0, then rndn2 assumes that
             -CV3 is the "A" matrix, and NOT the "-CV" matrix.  This
             can speed things up in the case of repetitive calls to rndn2.

proc rndn2_q(u,ccv,n);
    Like rndn2, but instead of providing the  variance matrix, provide the cholesky decomposition of the 
    variance matrix. 
       Ie.
	  vc=something;
	  ccv=chol(vc);
	  ps=rndn2_q(u,ccv,100);
    This is useful if you are calling rndn2 repetitively (rndn2_q is a bit quicker then rndn2)



**   proc pdfpois(lambda,k);
**        Compute PDF of Poisson random variates with parameter lambda,
**         at value k.
**
**  proc into_cat(avals,minct,maxct);
**       Divide avals into categories ranging from minct to maxct.
**       Each category has equal length.
** proc cdfpois(lambda,k);
**        Compute CDF of Poisson random variates with parameter lambda,
**          at value k.
**
** proc rndntrun(n,min,sd)
**   Returns n*1 vector of truncated normals; with min value min and sd=sd
**
**
**
**  proc inv_cdf(prob,density)
**    A general, empirical inverse CDF function.
**    You supply the densitys (defined at 1...N)
**    and you supply a probability
**    inv_cdf will return the smallest I s.t. the cumulative distribution
**    is greater then prob. I.e.: min I s.t. sumc(density[1:I))>prob.
**
**   proc rvpois(lambda)
**    Generate row(lambda) Poisson RVs, given parameter(s) lambda.
**
**    proc cdft1ev(x)   CDF of type 1 extreme value distribution
**    proc pdft1ev(x)   PDF of type 1 extreme value distribution
**    proc rndt1ev(n)   Return Nx1 vector of rvs from CDF of type 1
**                       extreme value distribution
**
**    proc lncdfn(z)   Log of CDFN or CDFNC.  Currently, this just traps
**     and lncdfnc(z)  the argument, and sets it equal to 8.3 if
**                     > 8.3 (CDFN) or 37 if > 37 (cdfnc).
**

** Proc varyindx(x)
Usage:
   okindx=varyindx(x)
where:
        x: a matrix
and
   okindx: an index into the columns of x

okindx will point to the columns of x that are "not constant" --
that have (non-missing) values that are NOT all the same.

If okindx=miss(1,1), then there are NO "non-constant" columns


**  proc hazardf(vals)
**
**             Returns the inverse mill's ratio = pdfn/cdfn, evaluated
** at vals.  If -16<vals<36, returns pdfn/cdfn.  For vals>36, returns
** 1e-300, for values <-16, returns abs(vals)+ 1/abs(vals) (an approximation
** with accuracy of < 0.0001 for values between -16 and -36,  so it seems
** to be very good)

Note: for pdfn(x)/(1-cdfn(x)); call hazardf(-x) 
**

Proc g_meanc(x,weights)
  Returns a weighted mean of x

Proc g_meanc2(x,weights)
  Returns a weighted mean of x -- with "rowwise" deletion of missing values.
  That is, all available information is used (just the missing values are removed)


proc g_vcx(x,weights)
  Returns a weighted covariance matrix of x

proc g_Vcx2(x,weights) 
  Returns a weighted covariance matrix of x -- with "rowwise" deletion of missing values.
  That is, all available information is used (just the missing values are removed)


proc (3)=g_chol_to_vc(chols,chols_vc)
   procedure to compute variance matrix, expected value of variance matrix,
   and the stderr^2 of each cell in this matrix,
   given estimates of  upper triangular cholesky decomposition of this matrix AND the
   VC matrix of these estimates.
   Each return is KxX, where K is the number of rows of chol (the number of "Variables")



proc g_rnd_halton(nobs,ndraws,nvars,options);
   Return an array, or matrix, of random numbers generated using a "halton sequence".
 
  Syntax:  eps_array=g_rnd_halton(nobs,ndraws,nvs,options);

  If nvs>1, this will be a 3d array.
   i.e.   foo=g_rnd_halton[100,40,3,""] will create an array with largest element foo[100,40,3]
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.

  The options are:

     TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.

    TYPES=typelist	       : typelist is a string of type identifiers.  Each id is a single character.
				 Valid identifiers are:
					N = normal
					E = extreme
					U = uniform
					T = triangle (-1 to 1 range)
					1 = uniform (-1 to 1 range)

			          If specified, the length of TYPES must equal NVS.
				  If specified, TYPE is ignored.
				
				  Example: TYPES=NEU1NN 



    RANDOM or RANDOM=seed     : Do a "randomized" halton numbers. A uniform random variable is added to
                                all "halton numbers"; the same value is used for all rows & draws for
                                a variable.
				Random is also used to create a scramble sequence for the SCRAMBLE option.

    SCRAMBLE                  : "Scramble" the halton sequence. 

    PRIMES=a,b,c,..           : Specify exactly which primes to use to create halton sequence.
                                a,b,c, etc must be prime numbers; with NO SPACES between the
                                = sign, the numbers, and the commas.
                                Example:  PRIMES=31,47,49
                                You MUST specify  NVS number of primes (otherwise, an error 
                                occurs). Thus, in the above example, NVS must be 3.

    PRIME_START=ith          : Specify where, in the list of 100 predefined primes, to start.
                               for example, PRIME_START=4 means "use 7 with the first variable, 11
                               for the second variable, etc."
                               By default PRIME_START=1 is used.
                               Note that if PRIMES is used, then PRIME_START will use the ith
                               prime specified in PRIMES as the first prime. In other words,
                               it is silly to use PRIME_START and PRIMES together.
     SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
     MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!

    ADD=n		    : Add, then remove, this multiple (n) draws to each observation
			      Thus, number of generated "draws" is NDRAWS*(n+1), but NDRAWS*n "columns" 
			      are tossed out. 
		              In particular, every nth column is used.
 			      If n=0, or is not specified, then n=3

    DISCARD=n		    : Discard n draws from beginning of the halton sequence (for each variable).



  Notes:
   *  A seperate vector of halton numbers is generated for each variable.
      Thus, for 100 rows, 10 draws, and 3 variables, there will be three
      1000 element halton-number vectors created.
   *  These vectors are chopped into contiguous subvectors, one subvector per observation
      (the subvector is used for the draws of a single observation & single variable)

   *  Each variable uses a different "prime number" to create the halton-number vector
   *  A max of 100 variables can be specified  (add numbers to the PRIM constant if you
       need more variables).

   *  To pull out a matrix of epsilons for all draws of the oth observation
         foo=arraytomat(eps_array[oth,.,.]);
   *  To pull out a vector of epsilons for the dth draw of the oth observation
         foo=arraytomat(eps_array[oth,dth,.]);
   *  To pull out a Nobs X Ndraws matrix of epsilons for the Vth variable:
         foo=arraytomat(atranspose(eps_array[.,.,Vth],3|1|2)); 

   * For larger Halton sequences, the SCRAMBLE option is ignored for low values of primes.
     In particular, if the # of numbers to be drawn  (for a given variable) is  > prime^8, then 
     scrambling will not be done.
     To force scrambling, you can use large prime numbers (use the PRIMES or PRIME_START options)
     Note that scrambling is most useful for shorter sequences (small number of draws) and larger
     primes (high values of NVS). 

   * For more than 25 variables, the behavior of the SCRAMBLE option is not quite kosher --
     the sorting sequence is not strictly increasing. This doesn't seem to diminish the SCRAMBLE
     features ability to fix the clumping problem when NVS is large (hence use of larger prime numbers)
     and a small number of draws, but the solution is something of a hack.
   

proc g_rnd_anti(nobs,ndraws,nvars,options);

   Return an array, or matrix, of random numbers generated using a "antithetics".

  Syntax: eps_array=g_rnd_anti(nrows,ndraws,nvs,options);

  if nvs>1, this will be a 3d matrix.
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.
  The options are:

  TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.
  ANTI= n           : the type of antithetic "replications"
     You can use 
        0 = None         --  no antithetic (one value per generated rv). This is included to 
			     facilitate comparisons, and to simplify creation of 3d arrays of rvs
        1 = All Negative -- generate 2 values per random  vector (the negative
  		            of the draw)
        2 = Negative     -- generate 2^NVS multiple values (per generated rv) by setting singles,
			    pairs, etc to their negative.
	3 = Shift        -- similar to 2, but also shift positions of last two variables
        4 = Scramble     -- Similar to 3, but scramble all positions

        The default type is 1 (all negative).

   SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
   MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!



proc g_rnd_mlhs(nrows,ndraws,nvs,options);

  Syntax: eps_array=g_rnd_mlhs(nrows,ndraws,nvs,options);
  If nvs>1, this will be a 3d array.
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.
  The options are:

     TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.


    RANDOM=seed            : Set the seed value for the randomized sequence.


   SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
   MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!


  This is a implementation of the method described in see Hess, Train and Polak
      "On the use of Modified Latin Hypercube Sampling (MLHS) Method in the Estimation of a Mixed 
       Logit model for vehicle choice".
  
  Basically, for the NDRAWS values (for one variable for one observation
    a) Generate a sequence of equally spaced values (0,1/NDRAWS,...,(NDRAWS-1)/NDRAWS)
    b) Add x, a uniform random value, to each number in the sequence (same value added to each);
       where 0 < x  < 1/NDRAWS
    c) Randomly reorder this sequence

   A collection of NROWS*NVS of these sequences is used to form the NROWS x NDRAWS x NVS array.




proc g_choose_noreplace(plist,ndo);


    syntax:  is_chosen=g_choose_noreplace(plist,ndo)

    where:
	  plist: Mx1 vector: probability of chosing an item (from a list)

	  ndo  : # of items to choose
          is_chosen: Mx1 0/1 dummy -- 1 if chosen (sumc(is_chosen) will equal NDO)

    This is used to "chose NDO items, from a Mx1 vector of items, without replacement.
    Plist is the probability of chosing an item -- higher values mean greater chance of 
    the item being chosen.

    
    Note that PLIST will be normalized.

    To choose J items out of M, with equal probability, 

         clist=g_choose_nreplace(ones(M,1),J)  ;





proc gc_distance(alat1,along1,alat2,along2);
    syntax:  gcdist=gc_distance(alat1,along1,alat2,along2);
  
    Computes the "great circle" (as the crow flies) distance between two latitude/longitude
    pairs (in the Western Hemisphere).  





**********************************/
proc g_dolog(x,expln0);
 local tmp;
 tmp=dolog(x,expln0);
 retp(tmp);
endp;


@ proc for taking log, with substitution for elements < 0  @
@ nb: all elements <=0 set to ln(expln0)  @
proc dolog(x,expln0);

local tmp,ind,k;
tmp=miss(zeros(rows(x),cols(x)),0);
if x>0;  retp(ln(x)); endif ;		@ trivial case @

for k (1,cols(x),1);
  ind=g_mkindex(x[.,k].>0);
  if not(scalmiss(ind)) ;
    tmp[ind,k]=ln(x[ind,k]);
  endif;
endfor;
if not(scalmiss(expln0)); tmp=missrv(tmp,ln(expln0)); endif;
retp(tmp);
 


         
@         retp(ln( ( (x .<= 0) .* expln0 ) + (x.*(x.>0)) )  );@
endp ;


/********************/
@ append a row to a matrix. IF matrix is empty or doesn't have right # of
  columns, then use the row as first line of the matrix @
  proc appnd_it(oldmtx,arow);
        if cols(arow)==1 ; arow=arow' ; endif ;  @ convert column to a row @
        if scalmiss(oldmtx);
           oldmtx=arow ;
        elseif cols(oldmtx) ne cols(arow);
           oldmtx=arow ;
        else;
            oldmtx=oldmtx|arow;
        endif;
     retp(oldmtx);
 endp;

/**********************************/
@ Generates multivariate normal, with mean U and variance matrix of CV. @
@ Returns an NxK matrix, K=rows(U).  Each row of N is distributed N(U,CV) @


proc (2)=rndn2(u,cv,n);

    local a,k;

@ Note: u=0 (scalar) or u=missing signals that U=0 @

    if scalmiss(u); u=0; endif;

    k=rows(u);
    if k==1 and u[1]==0;
        k=rows(cv); u=zeros(k,1);
    elseif k ne rows(cv);
        " % Error: mean vector and CV matrix of different size ";
        retp(miss(1,1),miss(1,1));
    endif;

@ to save time, if cv has negative on diagonal, then it is actually the
  A matrix, with A*A'=CV.  Otherwise, compute A @

    if cv[1,1]<0 ;
       a=-cv ;
    else ;
      trap 1 ;
       a=chol(cv) ;
       if scalerr(a);
          call grbl2_error(" Chol not work on vc  ",1); cv ; stop ;
        endif ;
        trap 0 ;
    endif;

    retp( ((rndn(n,k)*a)+u'),a);
endp;



/**********************************/
@ Generates multivariate normal, with mean U and cholesky decomposition of 
  variance matrix of CCV.
NOte: CCV*CCV'=VC
Note: ccv=chol(vc);

 @
@ Returns an NxK matrix, K=rows(U).  Each row of N is distributed N(U,CV) @


proc rndn2_q(u,ccv,n);

    local k;

@ Note: u=0 (scalar) or u=missing signals that U=0 @

    if missrv(u,0)==0; u=0; endif;

    k=rows(u);
    if k==1 and u[1]==0;
        k=rows(ccv); u=zeros(k,1);
    elseif k ne rows(ccv);
        " % Error: mean vector and cholesky matrix of different size ";
        retp(miss(1,1),miss(1,1));
    endif;

    retp( ((rndn(n,k)*ccv)+u'));

endp;


/*************/
@ pdf of mvn 
  vals=pdf_mvn(x,mu,vc)
      x: NxK vector of 1xK vectors to compute pdf's of
      mu: Kx1 vector of "means"
      vc: KxK covariance matrix
   
    vals: N x 1 vector of probabilities

 if Mu=scalar, use a Kx1 vector (each element = to mu)
 if VC=scalar, use a KxK diagonal matrix (each diagonal element=vc)
@

proc pdf_mvn(x,mu,vc);
local nx,kx,detvc,t1,t2,p1,tnum,tdenom,invvc,p1a,tmp,tmp1;

kx=cols(x);
nx=rows(x);

if kx==1 and nx>1 ; x=x'; endif;

if rows(mu)==1 and cols(mu)==1 ;
  mu=ones(kx,1)*mu ;
endif;
if rows(mu)==1 and cols(mu)==kx ;
  mu=mu';
endif;

if rows(vc)==1 and cols(vc)==1 ;
  vc=eye(kx)*vc ;
endif;

if kx ne rows(mu);
   call grbl2_error("Error in pdf_mvn: cols(X) ne size of mu " $+ g_ftos2(kx,1) $+ " " $+ g_ftos2(rows(mu),1),1);
  stop;
endif;

if kx ne cols(vc) or kx ne rows(vc);
  call grbl2_error("Error in pdf_mvn: cols(X) ne rows of cols of vc  " $+ g_ftos2(kx,1) $+ " " $+ g_ftos2(rows(vc),1) $+ 
      " " $+  g_ftos2(cols(vc),1),1);
  stop;
endif;

trap 1;
  invvc=invpd(vc);
  if scalerr(invvc) ;
    call grbl2_error("Error in pdf_mvn: unable to take inverse of VC matrix ",1);
    stop;
  endif;
trap 0;
detvc=detl ;  
 
t1=(2*pi)^(kx/2);
t2=sqrt(detvc);
tdenom=t1.*t2 ;
p1=x-mu';
 
tmp=zeros(rows(p1),1);
for j (1,rows(p1),1);
   p1a=p1[j,.]';
   tmp[j]=(-p1a'invvc*p1a)/2 ;
endfor;
p1a=exp(tmp)./tdenom;
retp(p1a);

endp;


/*************/
@ quicker pdf of mvn 
  vals=pdf_mvn_q(x,mu,invvc,tdenom)
      x: NxK vector of 1xK vectors to compute pdf's of
      mu: Kx1 vector of "means"
      invvc: inverse of  KxK covariance matrix
      tdenom: denominator of mvn = (2 Pi )^k/2 * det(vc) ^-0.5

   vals: N x 1 vector of probabilities

 if Mu=scalar, use a Kx1 vector (each element = to mu)
 if VC=scalar, use a KxK diagonal matrix (each diagonal element=vc)
@

proc pdf_mvn_q(x,mu,invvc,tdenom);
local nx,kx,detvc,t1,t2,p1,p1a,tmp,tmp1;

kx=cols(x);
nx=rows(x);

if kx==1 and nx>1 ; x=x'; endif;

if rows(mu)==1 and cols(mu)==1 ;
  mu=ones(kx,1)*mu ;
endif;
if rows(mu)==1 and cols(mu)==kx ;
  mu=mu';
endif;

if rows(invvc)==1 and cols(invvc)==1 ;
  invvc=eye(kx)*invvc ;
endif;

if kx ne rows(mu);
  call grbl2_error("Error in pdf_mvn_q: cols(X) ne size of mu " $+ g_ftos2(kx,1) $+ " " $+ g_ftos2(rows(mu),1),1);
  stop;
endif;

if kx ne cols(invvc) or kx ne rows(invvc);
  call grbl2_error("Error in pdf_mvn_q: cols(X) ne rows or cols of invvc " $+ g_ftos2(kx,1) $+ " " $+ g_ftos2(rows(invvc),1) 
    $+ " " $+ g_ftos2(cols(invvc),1),1) ;
  stop;
endif;

if missrv(tdenom,0)==0 ;
  call grbl2_error("Error in pdf_mvn_q: tdenom =0 " ,1);
  stop;
endif;

p1=x-mu';
tmp=zeros(rows(p1),1);
for j (1,rows(p1),1);
   p1a=p1[j,.]';
   tmp[j]=(-p1a'invvc*p1a)/2 ;
endfor;
p1a=exp(tmp)./tdenom;
retp(p1a);

endp;


proc (2)=pdf_mvn_factors(vc);
local t1,t2,tdenom,invvc,detvc,kx ;

kx=rows(vc);

trap 1;
invvc=invpd(vc);
if scalerr(invvc) ;     @problem, return missing @
  retp(miss(1,1),miss(1,1));
endif;
trap 0;
detvc=detl ;  
t1=(2*pi)^(kx/2);
t2=sqrt(detvc);
tdenom=t1.*t2 ;

retp(invvc,tdenom);
endp;



/*******************************************/
@ assign values to a category, from mincat to maxcat
 eg; if avals=0 5 4 2 3 1  and mincat=1 maxcat=3
      the result=1 3 3 2 3 1  (equal lengthed intervals)
@

proc into_cat(avals,minct,maxct);

local ncts,minv,maxv,diff,dvis,foo;

ncts=(maxct-minct)+1;
minv=minc(avals);
maxv=maxc(avals)+abs(minv/100000);
diff=maxv-minv ;
dvis=diff/ncts ;
foo= trunc((avals-minv)./dvis)+1 ;
retp(foo);
endp;


/******************/
@ Compute PDF of Poisson variates, given lambda and given K;

  pdf=exp(lambda)*lambda^k / k!

  If k is a scalar, then same k used for all lambda
  If lambda is a scalar, then same lambda used with all k
  Else, rows(lambda) MUST match rows(k)
  Both lambda and K MUST be Nx1 (or 1xN).
  Result is Nx1
@

proc pdfpois(lambda,k);
  local denom,t1,t2,pdfs,a,cdfs,lam,avals,i;

  if cols(k)>1 ; k=k'; endif;
  if cols(k)>1 ;
     call grbl2_error(" % CDFPOIS -- K must be a vector, not a matrix " ,1);stop;
  endif;
  if maxc(k)>100;
    call grbl2_error(" % CDFPOIS -- K too large (overflow above 100! ) ",1) ;
    stop;
 endif;

  if cols(lambda)>1 ; lambda=lambda'; endif;
  if cols(lambda)>1 ;
    call grbl2_error( " % CDFPOIS -- Lambda must be a vector, not a matrix ",1) ;stop;
  endif;


  if rows(k)==1  ;          @ uniform k value to evaluate at @
     t1=exp(-lambda).*(lambda.^ k);
     retp(t1./(k!)) ;

   else ;
     avals=zeros(rows(k),1);
     if (rows(k) ne rows(lambda)) and (rows(lambda) ne 1);
         call grbl2_error(" % POISCDF: rows(lambda) ne rows(k) ",1); stop ;
     endif;
     i=1;
     do until i>rows(k);
 @ many k, one lambda??  @
        if rows(lambda)==1 ; lam=1 ; else; lam=lambda[i]; endif ;
        avals[i]=exp(-lam)*(lam^k[i])/(k[i]!) ;
        i=i+1;
     endo;
     retp(avals);
    endif;
 endp;

/******************************/
@ generate CDF of poisson distribution, with parameter(s) lambda, at value=k

  cdf = Sum from j=0 to K : exp(-lambda)*lambda^j / j!

  If k is a scalar, then same k used for all lambda
  If lambda is a scalar, then same lambda used with all k
  Else, rows(lambda) MUST match rows(k)
  Both lambda and K MUST be Nx1 (or 1xN).
  Result is Nx1
@

proc cdfpois(lambda,k);

local denom,t1,t2,pdfs,a,cdfs,lam,avals,i;

@ just sum from 0..k, for each k, given lambda @

if cols(k)>1 ; k=k'; endif;
if cols(k)>1 ;
     call grbl2_error(" % CDFPOIS -- K must be a vector, not a matrix ",1) ;stop;
endif;

if rows(k)==1  ;          @ uniform value to evaluate at @
    if k==0;
      a=0;
    else ;
      a=seqa(0,1,k+1);
    endif;
    denom=a! ;
    t1=exp(-lambda);
    t2=lambda.^(a') ;
    pdfs=(t1.*t2)./(denom') ;
    cdfs=sumc(pdfs');
    retp(cdfs);

 else ;                 @ do 1 at a time @
     avals=zeros(rows(k),1);
     if (rows(k) ne rows(lambda)) and (rows(lambda) ne 1);
         call grbl2_error(" % POISCDF: rows(lambda) ne rows(k) ",1); stop ;
     endif;
     i=1;
     do until i>rows(k);
         a=seqa(0,1,k[i]+1);
         if rows(lambda)==1 ;
             lam=lambda ;
         else ;
             lam=lambda[i];
         endif;
         pdfs=exp(-lam).*(lam.^a)./(a!) ;
         avals[i]=sumc(pdfs);
         i=i+1;
      endo;
      retp(avals);
 endif;
endp;

/******************************/
@ Generate RVs from a Poisson distribution @
@ Works best with small lambda.  Uses the Poisson process assumption --
  N low PROBability uniform RVs are generated, s.t. N*prob=lambda.
  The number of sucesses (values of uniform RV < PROB) is used as the
  realization of the poisson RV @

  proc rvpois(lambda0,iout);

  local pois,foo,n,psucc,kk1,kk2,kk,idiv;
  local dotransp,dolam,lambda,amin,amax,ii,ichunk,igot,ar;

   igot=coreleft/48 ;       @ available workspace @

   ar=csrlin ;

   if rows(lambda0)>1 and cols(lambda0)>1 ;
       call grbl2_error( " % RV Poisson can only do vectors  ",1);
       stop ;
   elseif cols(lambda0)>1 ;
       lambda0=lambda0';
       dotransp=1 ;
    else;
       dotransp=0;
   endif ;

   idiv=8 ;       @ increase Idiv to approach ideal poisson process @
                  @ (but increases execution time @

   amax=ceil(maxc(lambda0));        @ max per read @
   ichunk=ceil((igot/idiv)/amax) ;

@ Algorithim: Divide lambda by a large N, s.t. lambda/N is small.
  Then, generate PS, a Nx1  vector of uniform RVs.  Realization is the number
  of times PS > lambda/N @

   kk=rows(lambda0);
   dolam=zeros(kk,1);

   kk1=1 ;
   do until kk1 ge kk ;
      if iout==1 ;
         locate ar,1 ;
        " Poisson RV .. " kk1;;
      endif;
      kk2=minc((kk1+ichunk)|kk);
      lambda=lambda0[kk1:kk2];

      amin=minc(lambda);
      amax=maxc(lambda);

      n=ceil(amax)*idiv;

      psucc=lambda/n ;
      foo=rndu(n,rows(lambda));
      pois=sumc(foo.<psucc');
      dolam[kk1:kk2]=pois ;

      kk1=kk1+ichunk ;
    endo;
    if iout==1; ? ; endif ;
    if dotransp==1 ;
       retp(dolam');
    else;
       retp(dolam);
    endif;
endp;


/******************************************/
@ General inverse cdf program.
  Given a list of "densitys" and a probability, find the n
  at which sum(densities from 1 to n) > prob.  Retp(n).
  Prob is between 0 and rows(list)
@


   proc inv_cdf(prob,wts);
    local inlist,lowerb,upperb,wtsum,wtcand,ncand,lowsum,totsum,tmp,do1,hisum;
    inlist=rows(wts);
    if not (wts ge 0 ) ;
        errorlog " % Negative weight in inv_cdf ";
        retp(miss(1,1));
    elseif prob le 0 ;
        retp(1);
    elseif prob ge 1 ;
        retp(inlist);
    endif;

    totsum=sumc(wts);


@ check this ... @
   if (wts[1]/totsum ge prob)  ;   @ if inlist==1 this also holds @
       retp(1);
   endif ;

   if inlist==2 ;
     retp(2);           @cause retp(1) done above @
   endif;

@ simple cases accounted for. Now, do a binary search for the proper n @
@ Note that ncand==1 will NEVER happen (see above @

    upperb=inlist ; lowerb=1; lowsum=wts[1] ; do1=0;

foo1: ;                         @ jump here for next candidate @
    if upperb-lowerb le 1 ;
       retp(upperb);
     endif ;

    ncand=trunc((upperb-lowerb)/2)+lowerb ;

    if do1==0;
       wtsum=lowsum+ sumc(wts[(lowerb+1):ncand]);
    else ;
      wtsum=hisum-sumc(wts[ncand+1:upperb]);
    endif ;
    wtcand=wtsum/totsum ;

    if wtcand > prob ;            @ gone too far @
        upperb=ncand ;
        hisum=wtsum;
        do1=1;
    else ;
        lowerb=ncand;
        lowsum=wtsum ;
        do1=0 ;
    endif ;
    goto foo1 ;

    endp ;

@ CDF of type 1 extreme value distribution @
proc cdft1ev(v);
   retp(exp(-exp(-v)));
endp;

@ pdf of type 1 extreme value distribution @
proc pdft1ev(v);
  retp(exp(-v-exp(-v)));
endp;

@ generate nx1 RVs from type 1 extreme value distribution @
proc rndt1ev(n);
  local foo,foo1;
  foo=rndu(n,1);
  retp(-ln(-ln(foo)));
endp;

@ return the log of the cdfn.  If z>8.3, return ln(cdfn(8.3)) @
proc lncdfn(z);

   if (z<8.2) ;
       retp(ln(cdfn(z))) ;

   else          ;   @  possible problems @

      retp( ln(cdfn( substute(z,z.>8.2,8.2))) ) ;

   endif;
   endp;

proc lncdfnc(z);

   if (z<37) ;
       retp(ln(cdfnc(z))) ;

   else          ;   @  possible problems @

      retp( ln(cdfnc( substute(z,z.>37,37))) ) ;

   endif;
   endp;



/*********************************
 add to reflist all elements of alist that are not in reflist.
  Also, augment counter of elements.
  If alist and counter are missing, initialize them.
  If new elements are added, expand counter.
  Note that counter and reflist ALWAYS maintain a row-wise match, even
  after reflist has been added to.
*/


      proc (2)=aug_list(alist,reflist,counter);
       local ind1,ind2,foo1,foo2,foo3,foo4,mm;

      if scalmiss(reflist) ; reflist=alist[1]; counter=0;  endif;
      if scalmiss(counter); counter=ones(rows(reflist),1,); endif;
      ind1=indnv(alist,reflist);
      if ismiss(ind1) ;   @ new elements.  Augment reflist and counter @
          foo1=ind1.==miss(1,1);
          foo2=g_mkindx(foo1);
          foo3=unique(alist[foo2],1);
          reflist=reflist|foo3 ;
          foo4=zeros(rows(foo3),1);
          counter=counter|foo4 ;
      endif;

@ augment counters @
      foo1=coreleft/16 ;
      if rows(reflist) * rows(alist) < foo1 ;
          foo2=sumc(reflist'.==alist) ;
          counter=counter+foo2;
     else;
         mm=1 ;
         do until mm> rows(reflist);
            counter[mm]=counter[mm]+sumc(alist.==reflist[mm]);
            mm=mm+1;
         endo;
     endif;

       retp(reflist,counter);

   endp;





/**********************************/
@ ---- check invmill for overflow condition --- @
proc hazardf(av);
    local tmp,tmpid,tmpsmall,tmpbig,amint,id,foo ;
@clearg hoy ; hoy=av;@
  if abs(av) lt 16  ;
      retp(pdfn(av)./cdfn(av)) ;
  endif;
  foo=miss(1,1);
 if ismiss(av) ;
    foo=av*0 ; foo=g_mkindx(missrv(foo,1));
    av=missrv(av,0) ;
  endif;
@" foo " foo';@
  tmp=av*0;
  amint=-16;

  tmpsmall= av .lt amint ;
  tmpid=g_mkindx(tmpsmall);
  if not(scalmiss(tmpid)) ;
      tmp[tmpid]=-av[tmpid]-(1./av[tmpid]) ;
  endif ;

  tmpbig= av .> 36 ;
  tmpid=g_mkindx(tmpbig);
  if not(scalmiss(tmpid)) ;
     tmp[tmpid]=1e-306+(tmp[tmpid]*0) ;
  endif ;

  tmpsmall=1-(tmpbig+tmpsmall);
  tmpid=g_mkindx(tmpsmall);
  if not(scalmiss(tmpid)) ;
      tmp[tmpid]=pdfn(av[tmpid])./cdfn(av[tmpid]) ;
  endif ;

  if not(scalmiss(foo)) ;           @reset missings @
     tmp[foo]=miss(zeros(rows(foo),1),0);
  endif;

 retp(tmp);

endp;


/**********************************/
@ Point to columns of x that are non-constant @
proc varyindx(x);

local ix,ind1,ac;
ix=0;ind1=miss(1,1);
do until ix==cols(x);
  ix=ix+1;
  ac=packr(x[.,ix]);
  if scalmiss(ac)==1 ; continue ; endif;
  if stdc(ac)==0 ; continue ; endif;
  ind1=ind1|ix;
endo;
retp(packr(ind1));
endp;


/***********************************/
proc rndntrun(n,amin,sd);

local co1,in1,ind2,mm,rng1,refs1,goon,goon2,nvs;

local denom,compto,rng,a,musd,epsx,refs,corrfact,stp,r1,r2,indx,arf;
musd=amin/sd;
if musd>7 ;     @ larg min -- all probability is on  musd @
   retp(ones(n,1)*musd);
endif;
if musd<-7 ;    @ small min -- ignore @
   retp(rndn(n,1));
endif;

corrfact=cdfn(musd);
denom=1-corrfact;

compto=rndu(n,1)     ;  @ inv cdf on this to generate rvs @

@ range to search between (Normallized by sd) @
r1=musd;
r2=maxc(5|(r1+2));
stp=(r2-r1)/10;

rng=seqa(r1,stp,11);

refs=0|((cdfn(rng)-corrfact)./denom)|1.1 ;
rng=-1e9|rng|1e9;

indx=indnvlt(compto,refs);

@format /ld 6,3;
arf=compto~indx~rng[indx]~rng[indx+1]~refs[indx]~refs[indx+1];
arf[1:minc(10|rows(arf)),.];@

mm=0;
goon=miss(1,1);goon2=goon;
do until mm==n ;
  mm=mm+1;
  co1=compto[mm];
  in1=indx[mm];
  r1=rng[in1]; r2=rng[in1+1];
  stp=(r2-r1)/33;
  rng1=seqa(r1,stp,34);
  refs1=((cdfn(rng1)-corrfact)./denom) ;
  ind2=indnvlt(co1,refs1);
  goon=goon|rng1[ind2];
  goon2=goon2|rng1[ind2+1];

@?;  co1~ind2~rng1[ind2]~rng1[ind2+1]~refs1[ind2]~refs1[ind2+1]~goon;@

endo;
goon=packr(goon); goon2=packr(goon2);
nvs=goon+(rndu(n,1).*(goon2-goon)) ;
retp(nvs);

endp;

@ ------------------------------------------- @
@ convert a decimal number to a binary number, return in an array
eg:  vv=inttobin(7) 
     vv would be (4x1): 0 1 1 1 
@
PROC INTTOBIN(N);
  if n==0 ; retp(0); endif;
  local t2,eek,t1;
  t2=n;
  eek=miss(1,1);
  do while 1==1;
     if t2==0 ; break ; endif;
     t1=t2%2 ;
     if t1 ==1 ;
        eek=eek|1;
     else;
        eek=eek|0;
     endif;
     t2=trunc(t2/2);
   endo;
retp(rev(packr(eek)));
endp;

/**************************/
/* Compare a vector of numbers to a sorted reference vector, return
a vector of values with each element containing a
pointer to the element of the reference vector that is just smaller
(or equal).

  indx=indnvlt(vals,refvals)


Example:
when
      vals=0  5  2  61 10  2311 622  133 9
   refvals=1 10 100 1000
then
   indx=indnvlt(vals,refvals)
yields
     indx= 0  1  1  2  2   4    3    3   1

Thus,
     refvals2[indx+1].<=vals, and
     refvals2[indx+2].>vals,
where
     refvals2=-infinity|refvals|infinity

Notes:
   * if refvals is unsorted, missing is returned.
   * if refvals and vals are not vectors, missing is returned
   * Implicitly, "refvals[0] is set equal to - infinity",
     hence all values of vals < refvals[1] yield indx values of 0

*/
proc indnvlt(avals,refvals);
local tmp,i1,goof2;

if cols(refvals)>1 ; refvals=refvals'; endif;
if cols(Refvals)>1 ; errorlog "Error: INDNVLT refvals argument not a vector " ;
   retp(miss(1,1));
endif;

tmp=sortc(refvals,1);
if not(tmp eq refvals);
  errorlog "Error: INDNVLT refvals not sorted ";
  retp(miss(1,1));
endif;

if cols(avals)>1 ; avals=avals' ; endif;
if cols(avals)>1 ; errorlog "Error: INDNVLT vals argument not a vector " ;
   retp(miss(1,1));
endif;


i1=ones(rows(avals),1);
tmp=i1.*.refvals';
goof2=dotfge(avals,tmp) ;
retp(sumc(goof2'));

endp;





/****** THE FOLLOWING RNDMVN functions are deprecated (use pdf_mvn instead)  ***/
/****** THE FOLLOWING RNDMVN functions are deprecated (use pdf_mvn instead)  ***/
/****** THE FOLLOWING RNDMVN functions are deprecated (use pdf_mvn instead)  ***/

/**********************************/
@ Generates multivariate normal, with mean U and variance matrix of CV. @
@ Returns an NxK matrix, K=rows(U).  Each row of N is distributed N(U,CV) 
   Note: checks for 0 variance items, returns the constant value (from U)@


proc rndMVN(u,cv,n);

    local a,k,fixit,iid,dd,tmp,thediag,cv0,arebad,t1,t2,imin;

    fixit=0;

    if rows(u)==1 and cols(u) gt 1 ; u=u' ; endif;

    arebad=miss(1,1);

@ Note: u=0 (scalar) or u=missing signals that U=0 @

    if scalmiss(u); u=0; endif;

    k=rows(u);
    if k==1 and u[1]==0;
        k=rows(cv); u=zeros(k,1);
    elseif k ne rows(cv);
        " % Error: mean vector and CV matrix of different size ";
        retp(miss(1,1),miss(1,1));
    endif;

    cv0=cv;
    thediag=diag(cv);
    dd=abs(thediag) .gt 0.00001 ;
    if dd eq 0 ; retp(u); endif;   @ degenerate case (no variance) @
    iid=seqa(1,1,k);

goo1: ;  @ jump here if chol fails @

    if not(dd eq 1) or not(scalmiss(arebad)) ; 
      fixit=diag(cv0) .== 0; 
      if not(scalmiss(arebad)); 
        fixit[arebad]=fixit[arebad]*0; 
        dd[arebad]=dd[arebad]*0;
      endif;
      fixit=g_mkindx(fixit);
      if rows(fixit)==k ;      @ all rows of cv are somehow bad @
	 call grbl2_error("Can not decompose variance matrix: " ,1); stop; 
      endif;
      iid=g_mkindx(dd) ;
      tmp=cv0[iid,.];
      cv=tmp[.,iid];
    endif;


@ Compute A, where A*A'=CV.@
     trap 1 ;  a=chol(cv) ;  trap 0 ;
     if scalerr(a);
        t1=diag(cv);
        imin=minindc(t1);
        t2=iid[imin];
        arebad=packr(arebad|t2);
	goto goo1 ;
     endif ;
 
    if fixit==0 ;
       retp( ((rndn(n,k)*a)+u') );
    endif;
    tmp=zeros(n,k);
    tmp[.,iid]=u[iid]'+(rndn(n,rows(a))*a) ;
    tmp[.,fixit]=tmp[.,fixit]+(u[fixit]');

    retp(tmp);


endp;

/**********************************/
@ Quickly generates multivariate normal, with mean U and variance matrix of A'A. 
  You should use rndMVN_a to compute A (from the VC matrix) @


proc rndMVN_q(u,a,n);

    local k;

    if rows(u)==1 and cols(u) gt 1 ; u=u' ; endif;

@ Note: u=0 (scalar) or u=missing signals that U=0 @

    if scalmiss(u); u=0; endif;

    k=rows(u);
    if k==1 and u[1]==0;
        k=rows(a); u=zeros(k,1);
    elseif k ne rows(a);
        " % Error: mean vector and CV matrix of different size ";
        retp(miss(1,1),miss(1,1));
    endif;

    retp( ((rndn(n,k)*a)+u'));
endp;


/**********************************/
@ Compute A, where A'A=CV.
  A can then be used in rndMVNq (quick version of rndMVN)@


proc rndMVN_a(cv);

  local a,k;
   trap 1 ;
   a=chol(cv) ;
   if scalerr(a);
          call grbl2_error(" Chol not work matrix ",1);  cv ; stop ;
     endif ;
   trap 0 ;
    retp(a);
endp;


/****** END OF RNDMVN functions  ***/


/*************************/
@ weighted vc matrix @
proc g_Vcx(x,w);
local tx,tm,dev,dev2;

if rows(w)==1 ; w=ones(rows(x),1).*w ; endif;

tx=x.*w ;
tm=sumc(tx)./(sumc(w));

dev=x-tm';
dev2=dev'(dev.*w) ;
retp(dev2/(sumc(w)));
endp;



/*****************/
@ weighted mean @
proc g_meanc(x,w);
local tx,tm;
if rows(w)==1 ; w=ones(rows(x),1).*w ; endif;

tx=x.*w ;
tm=sumc(tx)./(sumc(w));
retp(Tm);
endp;



/*************************/
@ weighted vc matrix, with listwise deletion (use all possible data, just excluding rows with missings @
proc g_Vcx2(x,w);
local tx,tm,dev1,dev2;
local k,j,xx,aa;
if rows(w)==1 ; w=ones(rows(x),1).*w ; endif;

tx=g_meanc2(x,w) ;

aa=miss(zeros(cols(x),cols(x)),0);
for j (1,cols(x),1);
  for k (1,cols(x),1);
    xx=packr(x[.,j|k]~w) ;
    if scalmiss(xx) or scalmiss(tx[j]) or scalmiss(tx[k]) ; 
	continue ; 
    endif ;
    dev1=xx[.,1]-tx[j]; 
    dev2=xx[.,2]-tx[k];
    aa[j,k]=( (dev1'(dev2.*xx[.,3]))  )/ (sumc(xx[.,3]));
  endfor;
endfor;

retp(aa);

endp;



/*****************/
@ weighted mean -- listwise deletion of missings @
proc g_meanc2(x,w);
local tx,tm;
local aa,j,x0,xw ;
if rows(w)==1 ; w=ones(rows(x),1).*w ; endif;

aa=miss(zeros(cols(x),1),0);

for j (1,cols(x),1);
 x0=x[.,j];
 xw=packr(x0~w);
 if scalmiss(xw) ; continue; endif;
 tx=xw[.,1].*xw[.,2] ;
 tm=sumc(tx)./(sumc(xw[.,2]));
 aa[j]=tm;
endfor;
retp(aa);
endp;


/**************/
@ Ln, with 0 check and "don't do dummies" check@

proc g_ln(x,expln0);
local i,a1,ia1,lnx,xnew,isdummy;
if upper(expln0)$=="E" ; expln0=exp(1); endif;
for i (1,cols(x),1);
  a1=x[.,i];
  ia1=unique(a1,1);
  isdummy=1;
  if rows(ia1) gt 2 ;
     isdummy=0;
  else;
     if ia1[1] ne 0 and ia1[1] ne 1; isdummy=0 ; endif;
     if isdummy==1 and rows(ia1)==2 ;
        if ia1[2] ne 0 and ia1[2] ne 1; isdummy=0 ; endif;
     endif;
  endif;

 if isdummy==0 ;
      if a1>0;
          lnx=ln(a1);
       else;
         lnx=(ln( ( (a1 .<= 0) .* expln0 ) + (a1.*(a1.>0)) )  );
       endif;
  else ;
     lnx=a1;
  endif;
  if i==1 ; xnew=lnx ; else ; xnew=xnew~lnx ; endif;

endfor;
retp(xnew);

  
endp ;


/*************************/
@ Gamma(X+1)= X! 
  LogGamma=Ln(gamma(x+1))  = ln(X!) = lnfact(x) 
Hence logGamma(x) = lnfact(x-1)

HOWEVER: 
   lnfact(z<0) is undefined. Hence, for x<1.0, directly compute ln(gamma(x)) 
AND
   gamma(z>169) overflosw. Hence. for X>100, use lnfact(x-1)

@

proc g_loggamma(x);
local res,ish,isl,i1;

if maxc(x)<100 ;
  retp(ln(gamma(x))) ;
elseif minc(x) ge 1.0 ;
     retp(lnfact(x-1)) ;

else ;			@ some lt 1 or some gt 120 @
  i1=x.<2  ;
  isL=g_mkindx(i1);
  isH=g_mkindx(i1.==0);
  res=zeros(rows(x),1);
  if not(scalmiss(isL));
     res[isL]=ln(gamma(x[isL]));
  endif;
  if not(scalmiss(isH));
     res[isH]=lnfact(x[isH]);
  endif;
  retp(res);
endif;

endp;


/******************************/
/* procedure to compute variance matrix, and the stderr of each cell in this matrix,
   given estimates of  upper triangular cholesky decomposition of this matrix ANd the
   VC matrix of these estimates.

  Syntax:
      {vc,vc_exp,vc_se2}=g_chol_to_vc(chols,chols_vc)
   where
      vc = K x K variance
      vc_exp = expected value of vc
     vc_se2 = K x K matrix of the squared standard errors of the corresponding element of vc

     chols =  Either a K x K upper triangular matrix (the cholesky decomposition to use to create VC)
	      or a (K+1)*K/2 vector containing the row-wise VEC of this upper triangle.
	      For example, for a 3 x 3 cholesky upper triangular matrix
			  c11 c12 c13
			  0   c22 c23
		          0   0   c33
	      the 6x1 vector-version would be: {c11 c12 c13 c22 c23 c33}

    chols_vc  = KK x KK covariance matrix of elements in chols; where KK = (K+1)*K/2
	        If each element is independent, you can use a KK x 1 matrix (the variance terms)
		For example, for the 3 x 3 cholesky matrix, which contains 6 terms, the 6 x 6 matrix would be
		
		   s11,11    s11,12  s11,13  s11,22  s11,23  s11,33	        
		   s12,11    s12,12  s12,13  s12,22  s12,23  s12,33	        
		   s13,11    s13,12  s13,13  s13,22  s13,23  s13,33	        
		   s22,11    s22,12  s22,13  s22,22  s22,23  s22,33	        
		   s23,11    s23,12  s23,13  s23,22  s23,23  s23,33	        
		   s33,11    s33,12  s33,13  s33,22  s33,23  s33,33	        

                where sij,kl is the "estimated covariance term between the cij and ckl coefficients (of chols)"


	      Note: the cij terms are "estimates", with chols_vc being the VC of these estimates
                
Example of usage: 

i5=4;
let l[4,4]=11 12 13 14 0 22 23 24 0 0 33 34 0 0 0 44 ;
ndo=(i5+1)*i5/2;
arf=rndn(20,ndo);
goo2=vcx(arf);        
{a,ae,b}=g_chol_to_vc(L,goo2);
end;
============================================= */


@----------------@
Proc (3)= g_chol_to_vc(chols,chols_vc);

local t,t2,neps,ii;
local nchols,epsvc,epsvc_exp;
local ltable;
LOCAL EPSVC_SE,vir,vic,pairlist,ir,ic,ip,t1r,t1c,t2r,t2c,c1,c2;
local l1,l2,v_c1c1,v_C2c2,v_c1c2,tt;
 
if cols(chols)==1 ;   	@ in vector form @
  t=2*rows(chols);
  t2=t*2 ;
  neps=(sqrt(1+4*t)-1)/2 ;
  ii=0;
  t=zeros(neps,neps);
  for ir (1,neps,1) ;
    for ic (ir,neps,1) ;
      ii=ii+1;
      t[ir,ic]=chols[ii];
    endfor;
  endfor;
  chols=t; 
 
else;
  neps=rows(chols);  	 @ number of "epsilons" @
  t=lowmat1(chols);
  if not (t eq eye(neps));
    call grbl2_error("ERROR: cholesky matrix must be upper triangular ",1); stop ;
  endif;
endif;


nchols=(neps+1)*neps/2 ;   @ # of terms in chol matrix @

if cols(chols_vc)==1 ;       @ iid terms, create matrix with this on diagonal@
  if rows(chols_vc) ne nchols ; format /ld 3,0;
    call grbl2_error("ERROR: number of elements in diagonal chols_vc (" $+ g_ftos2(rows(chols_vc),1) $+ 
         ") ne to number of terms in chols (" $+ g_ftos2(nchols,1),1);
    stop ;
  endif;
  t=eye(nchols);
  chols_vc=diagrv(t,chols_vc);
else;
  if rows(chols_vc) ne nchols ;
    call grbl2_error("ERROR: rows of chols_vc (" $+ g_ftos2(rows(chols_vc),1) $+
        ") ne to number of terms in chols (" $+ g_ftos2(nchols,1),1);
  endif;
endif;
 
@ done with checking @

@ Easy part: compute the variance matrix@ 

epsvc=chols'chols      ;


/*
 Computing EPSVC_SE: the standard errors (actually, square of standard errors) of each term in epsvc
 Several assumptions are used:
 1) Each term in epsvc is a dot product--- a sum or multiplicative pairs
       Example:  epsvc_23 =  chol_12*chol_13 +  chol_22* chol_23
 2) Simplfication: the variance of each term in epsvc is the sum of the variances of its pairs.
    That is, we ignore covariances between the pairs.
 3) The variance of each pair is estimated using a first order linear approximation (Mood et al pg 181)
    	    Var(g(X,Y)=  [var(X)  (dg/dx)^2 ]   +   [var(Y) (dg/dy)^2 ]   +   [2 cov(X,Y) dg/dx dg/dy ]
 4) Given that G(X,Y)= X*Y, 
		dg/dx = Y  and dg/dy = X 
 
 5)  Example, for the first pair in epsvc_SE_23: 
       Var(chol_12*chol_13) =
	     CHOLVC_12,12 * CHOL _13^2  +  CHOLVC_13,13 * CHOL_12^2  +  2* CHOLVC_12,13 * CHOL_13 * CHOL_12
     where cholvc_xx,yy is the xx,yy term from ovariance matrix of the estimated cholesky factors


    Reminder:
      If there are J  "varying coefficients" there will be j "epsilons"
      The variance (and covariances) of these j epsilons can be computed from an estimated a matrix of
      "cholesky" factors
      There are j rows in this matrix of cholesky factors
      As an upper diagonal matrix, there are (J+1)*J/2 unique elements in this matrix
      Each of these unique elements corresponds to a cholesky factor.
      Each element is estimated as if it were a seperate coefficient.
      Hence, the beta vector produced by the MLE contains the "beta" terms (for both varying and nonvarying coefficents)
      and these "cholesky factors"
      Furthermore, the VC matrix produced by the MLE contains variances for the beta terms and the cholesky factors
      It also contains covariance terms between cholesky factors.
      Example: if there are 3 varying parameters and 4 non varying parameters
	   There will 3 epsilons
	   There will be 6 cholesky factors
           The cholesky matrix used to generate estimates of the variance of the epsilons would
           be an upper diagonal 3x3 matrix, with the 6 cholesky factors filling the upper triangle.
           MLE will generate a 13x1 vector of coefficients
           MLE will generate a 13 x 13 variance/covariance matrix of these coefficients
           One can extract the 6x6 submatrix corresponding to the variance/covariance of the cholesky factors.
           CHOLVC is this submatrix


  Computing EPSVC_EXP, the expected value of the VC matrix (EPSVC).
  EPSVC is computed by a matrix product of the estimated cholesky matrix.
  However, since these factors are estimated, the expected value will depend on covariances.
  In particualr we use E(XY)= ux + uy + cov(x,y),
  and E(A + B) = E(A) + E (B)

   
    
*/

@ create a lookup table -- from ir,ic element of chols to corresponding row of chols_vc @
ltable=eye(neps);
t=0;
for vir (1,neps,1);
  for vic (vir,neps,1);
    t=t+1;
    ltable[vir,vic]=t;
 endfor;
endfor;


@ For each element in epsvc_se (actually, in upper triangle) @
epsvc_se=zeros(neps,neps); epsvc_exp=epsvc_se;
 
for vir (1,neps,1);
  for vic (vir,neps,1);
 
@ figure out the list of pairs contributing to this element
  Example: for epsvc_23:  chol_12*chol_13 +  chol_22* chol_23 
          vir=2, vic=3; hence pairlist will be:
             1 2 1 3
  	     2 2 2 3 
@
  pairlist=miss(zeros(1,4),0);
    for ir (1,vir,1) ;
	pairlist=pairlist|(ir~vir~ir~vic)  ;  
    endfor;
    pairlist=packr(pairlist);
 
@ for each row in pairlist, compute variance of its two terms @
   for ip (1,rows(pairlist),1);
     t1r=pairlist[ip,1];   t1c=pairlist[ip,2];
     t2r=pairlist[ip,3];   t2c=pairlist[ip,4];

     c1=chols[t1r,t1c] ; c2=chols[t2r,t2c];

     l1=ltable[t1r,t1c] ; l2=ltable[t2r,t2c];
     v_c1c1=chols_vc[l1,l1] ; v_c2c2=chols_vc[l2,l2];
     v_c1c2=chols_vc[l1,l2];

@ reminder: terms are of the form
     CHOLVC_12,12 * CHOL_13^2  +  CHOLVC_13,13 * CHOL_12^2  + 2 * CHOLVC_12,13 * CHOL_13 * CHOL_12
@
     tt=(v_c1c1*l2*L2) + (v_c2c2*l1*l1) + (2* v_c1c2*l1 *l2);
     epsvc_se[vir,vic]=epsvc_se[vir,vic]+ tt ;   

     epsvc_exp[vir,vic]=epsvc_exp[vir,vic]+((c1*c2) + v_c1c2 );

@l1 l2 v_c1c2;@

   endfor;		@ ip @
   epsvc_se[vic,vir]=epsvc_se[vir,vic] ;
   epsvc_exp[vic,vir]=epsvc_exp[vir,vic] ;

 endfor;
endfor;

retp(epsvc,epsvc_exp,epsvc_se);


endp;

/********************/
@ returns a matrix of random numbers generated using antithetic draws 

  eps_array=g_rnd_anti(nrows,ndraws,nvs,options);

  If nvs>1, this will be a 3d array.
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.
  The options are:

  TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.
  ANTI= n           : the type of antithetic "replications"
     n specifies the type. You can use 
        0 = None      --  no antithetic (one value per generated rv). This is included to 
			   facilitate comparisons, and to simplify creation of 3d arrays of rvs)
        1=All Negative -- generate 2 values per random  vector (the negative
  		           of the draw
        2=Negative    --  generate 2**NVS multiple values per vector by setting singles,
			  pairs, etc to their negative.
	3=Shift       -- similar to 2, but also shift positions of last two variables
        4=Scramble   -- Similar to 3, but scramble all positions

        The default type is 1 (all negative).

   SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
   MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!


@
proc g_rnd_anti(nrows,ndraws,nvs,opts);

local tlist,tsd,tmean,dtype,i1,i2,i3,atype,antitype;
local tmp,eps,eps1,eps2,iscramble,todo;

local eek,t1,t2,iat;

let tlist=NORMAL EXTREME UNIFORM  ;
tsd=1; tmean=0;
opts=" " $+  upper(opts)$+ " " ;
for i3 (1,rows(tlist),1);
  tmp=" " $+ tlist[i3] $+" ";
  i1=strindx(opts,tmp,1);
  if i1>0 ;
       opts=opts $+ " TYPE=" $+ tlist[i3] $+ " " ;
   endif;
endfor;   

dtype=1 ;  @default of normal rvs @
i1=strindx(opts," TYPE=",1);
if i1>0 ;
   i2=strindx(opts,"=",i1);
   i3=strindx(opts," ",i2);
   atype=strsect(opts,i2+1,i3-(i2+1)) ;
   dtype=missrv(indcv(atype,tlist),1);   @ default is normal @
endif;
 
antitype=0;
i1=strindx(opts," ANTI=",1);
if i1>0 ;
  i2=strindx(opts,"=",i1);
  i3=strindx(opts," ",i2);
  antitype=strsect(opts,i2+1,i3-(i2+1)) ;
  antitype=missrv(indcv(antitype,"1"|"2"|"3"|"4"|"0"),1);
endif;

i1=strindx(opts," SD=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tsd=stof(tmp) ; endif;
endif;

i1=strindx(opts," MEAN=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tmean=stof(tmp) ; endif;
endif;

@ generate one uniform draw, of nvs rvs, at a atime -- 
  then augment it using the antithetic rule
  then transform 
  then store in proper obs/draws @

eps={} ;

todo=nrows*ndraws ;
eps=zeros(todo+2,nvs);
iat=0;

do until iat gt todo  ;

eps1=rndu(nvs,1);

if antitype gt 1 and antitype lt 5 ; @negative combos, or negative with swaps @
  if antitype==4 ; iscramble=sortind(rndu(nvs,1)); endif;

  for i1 (0,(2^nvs)-1,1) ;
    iat=iat+1;
    if iat gt todo ; break ; endif;
@ create binary representation of a number --- @
    eek={};
    t2=i1;
    for i2 (1,10000,1);
      if t2==0 ; break ; endif;
      t1=t2%2 ;
      if t1 ==1 ;
         eek=1|eek;
      else;
         eek=0|eek;
      endif;
      t2=trunc(t2/2);
   endfor;
 
   if rows(eek) < nvs ;
     eek=zeros(nvs-rows(eek),1)|eek;
   endif;
   eek=eek-(eek.==0);

   if antitype==2 or ( (antitype>2 and antitype<5) and (nvs eq 1)  ) ;
      eps[iat,.]=(eek.*eps1)';

   elseif antitype==3 ;
      eps2=eps1;
 
      if sumc(eek[nvs-1:nvs])==0 ;
         tmp=eps2[nvs] ; 
         eps2[nvs]=eps2[nvs-1]; 
         eps2[nvs-1]=tmp;
      endif;
      eps[iat,.]=(eek.*eps2)';

   elseif antitype==4 ;
      eps2=eps1[iscramble];
      eps[iat,.]=(eek.*eps2)';

   endif;
 endfor;

elseif antitype==5;  @ no antithetic @
    iat=iat+1;
    if iat gt todo ; break ; endif;
    eps[iat,.]=eps1' ;

else 	;	@ full negative @
    iat=iat+1;
    if iat gt todo ; break ; endif;
    eps[iat,.]=(eps1');
    iat=iat+1;
    if iat gt todo ; break ; endif;
    eps[iat,.]=(-eps1');
endif;
endo;			@ rows(eps)> todo @
eps=eps[1:todo,.];

if antitype ne 5 ;
  eps=(eps/2)+0.5;	  @ put into 0-1 range @ 
endif;

@ sort into a matrix, or array @
if dtype==1 ;  @ normal random variable @
     eps=cdfni(eps-(1e-10));		@ hack to avoid resolusion glitches @
     eps=eps.*(eps .<= 10)+10.*(eps .> 10);          /*truncates the normal dist. at 10 and -10*/
     eps=eps.*(eps .>= -10)-10.*(eps .< -10);
elseif dtype==2 ;   @ extreme valued @
     eps=-ln(-ln(eps));
     eps=eps.*(eps .<= 10)+10.*(eps .> 10);          /*truncates the extreme value dist. at 10 and -10*/
     eps=eps.*(eps .>= -10)-10.*(eps .< -10);
endif;
if tsd ne 1 ; eps=eps*tsd ; endif;
if tmean ne 0 ; eps=eps+tmean ; endif;

if nvs==1 ;
  eps=reshape(eps,nrows,ndraws);
else ;
  eps=areshape(eps,nrows|ndraws|nvs);
endif;


retp(eps);
endp;


/********************/
@ returns a matrix of random numbers generated using halton numbers.

Syntax:
   eps_array=g_rnd_halton(nrows,ndraws,nvs,options);

  If nvs>1, this will be a 3d array.
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.
  The options are:
     TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.
    TYPES=typelist	       : typelist is a string of type identifiers.  Each id is a single character.
				 Valid identifiers are:
					N = normal
					E = extreme
					U = uniform
					T = triangle (-1 to 1 range)
					1 = uniform (-1 to 1 range)

			          If specified, the length of TYPES must equal NVS.
				  If specified, TYPE is ignored.				
				  Example: TYPES=NEIU1NN 


    RANDOM or RANDOM=seed     : Do a "randomized" halton numbers. A uniform random variable is added to
                                all "halton numbers"; the same value is used for all rows & draws for
                                a variable.
				Random is also used to create a scramble sequence for the SCRAMBLE option.

    SCRAMBLE                  : "Scramble" the halton sequence. 

    PRIMES=a,b,c,..           : Specify exactly which primes to use to create halton sequence.
                                a,b,c, etc must be prime numbers; with NO SPACES between the
                                = sign, the numbers, and the commas.
                                Example:  PRIMES=31,47,49
                                You MUST specify  NVS number of primes (otherwise, an error 
                                occurs). Thus, in the above example, NVS must be 3.

    PRIME_START=ith          : Specify where, in the list of 100 predefined primes, to start.
                               for example, PRIME_START=4 means "use 7 with the first variable, 11
                               for the second variable, etc."
                               By default PRIME_START=1 is used.
                               Note that if PRIMES is used, then PRIME_START will use the ith
                               prime specified in PRIMES as the first prime. In other words,
                               it is silly to use PRIME_START and PRIMES together.
     SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
     MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!

    ADD=n		    : Add, then remove, this multiple (n) draws to each observation
			      Thus, number of generated "draws" is NDRAWS*(n+1), but NDRAWS*n "columns" 
			      are tossed out.
		              In particular, every nth column is used.
 			      If n=0, or is not specified, then n=3

    DISCARD=n		    : Discard n draws from beginning of the halton sequence (for each variable).

Notes:
   *  A seperate vector of halton numbers is generated for each variable.
      Thus, for 100 rows, 10 draws, and 3 variables, there will be three
      1000 element halton-number vectors created.
   *  These vectors are chopped into contiguous subvectors, one subvector per observation
      (the subvector is used for the draws of a single observation & single variable)

   *  Each variable uses a different "prime number" to create the halton-number vector
   *  A max of 100 variables can be specified  (add numbers to the PRIM constant if you
       need more variables).

   *  To pull out a matrix of epsilons for all draws of the oth observation
         foo=arraytomat(eps_array[oth,.,.]);
   *  To pull out a vector of epsilons for the dth draw of the oth observation
         foo=arraytomat(eps_array[oth,dth,.]);
   *  To pull out a Nrows X Ndraws matrix of epsilons for the Vth variable:
         foo=arraytomat(atranspose(eps_array[.,.,Vth],3|1|2));

   * The SCRAMBLE and ADD options are manipulations that are supposed to
     cure problems that arise when NVS is large and NDRAWs is small.

   * For larger Halton sequences, the SCRAMBLE option is ignored for low values of primes.
     In particular, if the # of numbers to be drawn  (for a given variable) is  > prime^8, then
     scrambling will not be done.
     To force scrambling, you can use large prime numbers (use the PRIMES or PRIME_START options)
     Note that scrambling is most useful for shorter sequences (small number of draws) and larger
     primes (high values of NVS).
  * For more than 25 variables, the behavior of the SCRAMBLE option is not quite kosher --
    the sorting sequence is not strictly increasing. This doesn't seem to diminish the SCRAMBLE
    features ability to fix the clumping problem when NVS is large (hence use of larger prime numbers)
    and a small number of draws, but the solution is something of a hack.

@
proc g_rnd_halton(nrows,ndraws,nvs,opts);

local res,hh,hm1,phi,i,y,x,j,s;
local prim,todo;
local dtype,i1,i2,i3;
local dorandom,gotseed,arand;
local aprimes,ap1,ap2,plist,a1;
local iprime,tlist,atype,tmp,tsd,tmean,doscramble;
local iat,iat2;
local scramble,xscramble,yscramble,ascrs,foo ;
local eps,maxscram;
local ndiscard,moredraws,diddraws ;
local goo,tlist2,typelist;
if nvs >100 ;
 call grbl2_error("Error: maximum of 100 variables can be created by proc rndn_halton (specified= " $+ g_ftos2(nvs,1),1);
 stop;
endif;

if nvs==0 ;   
 call grbl2_error(" ERROR: g_rnd_halton called with nvs=0 ",1); stop; 
endif;

let tlist=NORMAL EXTREME UNIFORM TRIANGLE ONE  ;   @ last 2 are -1 to 1 range @
let tlist2="N" "E" "U" "T" "1" ;

tsd=1; tmean=0;

@ Some additional prime numbers ---
569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 
683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 
827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 
971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049 1051 1061 1063 1069 
1087 1091 1093 1097 1103 1109 1117 1123 1129 1151 1153 1163 1171 1181 1187 1193 
1201 1213 1217 1223 1229 1231 1237 1249 1259 1277 1279 1283 1289 1291 1297 1301 
1303 1307 1319 1321 1327 1361 1367 1373 1381 1399 1409 1423 1427 1429 1433 1439 
1447 1451 1453 1459 1471 1481 1483 1487 1489 1493 1499 1511 1523 1531 1543 1549 
1553 1559 1567 1571 1579 1583 1597 1601 1607 1609 1613 1619 1621 1627 1637 1657 
1663 1667 1669 1693 1697 1699 1709 1721 1723 1733 1741 1747 1753 1759 1777 1783 
1787 1789 1801 1811 1823 1831 1847 1861 1867 1871 1873 1877 1879 1889 1901 1907
1913 1931 1933 1949 1951 1973 1979 1987 1993 1997 1999 2003 2011 2017 2027 2029
 2039 2053 2063 2069 2081 2083 2087 2089 2099 2111 2113 2129 2131 2137 2141 2143
 2153 2161 @



let prim =2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 
          59 61 67 71 73 79 83 89 97 101 103 107 109 
          113 127 131 137 139 149 151 157 163 167 173 
          179 181 191 193 197 199 211 223 227 229 233 239 
          241 251 257 263 269 271 277 281 283 293 307 311 
          313 317 331 337 347 349 353 359 367 373 379 383 
          389 397 401 409 419 421 431 433 439 443 449 457 
          461 463 467 479 487 491 499 503 509 521 523 541 547 557 563   ;

ascrs=zeros(100,1);
for i (1,100,1);
  ascrs[i]=chrs(i+32);
endfor;
ascrs=ascrs|ascrs|ascrs|ascrs|ascrs|ascrs|ascrs ;  @ a hack to handle up to 700 divisions (for 100 or so variables @
@for i (1,rows(@

foo=rndu(rows(ascrs),1);
foo=sortindc(foo);
ascrs=ascrs[foo];    @ this is the scrambled order @


opts=" " $+ upper(opts) $+ " ";

@ replace some stuff in options @
i1=strindx(opts,"IPRIME=",1);
if i1>0 ;
    i2=strindx(opts,"=",i1);
    a1=strsect(opts,1,i1-1);
    opts=a1 $+ " PRIME_START=" + strsect(opts,i2+1,10000);
endif;    

i1=strindx(opts," RANDOM ",1);
if i1>0 ;
    opts=opts $+ " RANDOM= " ;
endif;    


for i3 (1,rows(tlist),1);
  tmp=" " $+ tlist[i3] $+" ";
  i1=strindx(opts,tmp,1);
  if i1>0 ;
       opts=opts $+ " TYPE=" $+ tlist[i3] $+ " " ;
   endif;
endfor;   


@--- done with synonyms replacement  @


doscramble=0;
i1=strindx(opts," SCRAMBLE ",1);
i2=strindx(opts," SCRAMBLED ",1);

if (I2+i1)>0; doscramble=1; endif;

typelist=ones(nvs,1);  @ default of normal rvs @

i1=strindx(opts," TYPE=",1);
if i1>0 ;
   i2=strindx(opts,"=",i1);
   i3=strindx(opts," ",i2);
   atype=strsect(opts,i2+1,i3-(i2+1)) ;
   dtype=missrv(indcv(atype,tlist),1);   @ default is normal @
   typelist=ones(nvs,1)*dtype ;
endif;

i1=strindx(opts," TYPELIST=",1);

if i1>0 ;
   i2=strindx(opts,"=",i1);
   i3=strindx(opts," ",i2);
   tmp=g_strip(strsect(opts,i2+1,i3-(i2+1))) ;
   if strlen(tmp) ne nvs ;
	call grbl2_error("ERROR in halton: length of TYPELIST ne #variables " ,1);
	stop ;
   endif;
   for i (1,nvs,1);
	typelist[i]=missrv(indcv(strsect(tmp,i,1),tlist2),1); 		@ use normal as default (if unrecognized typelist id @
   endfor;
endif;
 

dorandom=0;gotseed=miss(1,1);
i1=strindx(opts," RANDOM=",1);
if i1>0 ;
  dorandom=1 ; 
  i2=strindx(opts,"=",i1);
  i3=strindx(opts," ",i2);
  arand=strsect(opts,i2+1,i3-(i2+1)) ;
  if arand$/="" ;
     gotseed=stof(arand);
     if ismiss(gotseed) ;
             call grbl2_error(" Bad random number specified in g_rnd_halton: " $+arand,1) ; stop ;
     endif;
     rndseed gotseed ;
  endif;    @ something after RANDOM= @
endif;

ndiscard=0;    @ discard from front of halton sequence @
moredraws=0;  @ extra draws for each "Row" (observation) @

i1=strindx(opts," DISCARD=",1);
if i1>0 ;
  i2=strindx(opts,"=",i1);
  i3=strindx(opts," ",i2);
  i1=strsect(opts,i2+1,i3-(i2+1)) ;
  if i1$/="" ;
     ndiscard=stof(i1);
  endif;
endif;


i1=strindx(opts," ADD=",1);
if i1>0 ;
  i2=strindx(opts,"=",i1);
  i3=strindx(opts," ",i2);
  i1=strsect(opts,i2+1,i3-(i2+1)) ;
  if i1$=="." or i1$=="" ; i1="0" ;endif;
   moredraws=stof(i1);
  if moredraws==0 ; moredraws=3 ; endif;
endif;



i1=strindx(opts," PRIMES=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     aprimes=strsect(opts,i2+1,i3-(i2+1)) ;
     plist=stof(aprimes);
     prim=plist;
     if not(scalmiss(indnv(1,plist)));
        call grbl2_error("1.0 specified as a PRIME in G_RND_HALTON: " $+aprimes,1);
         stop;
     endif;
 endif;

i1=strindx(opts," SD=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tsd=stof(tmp) ; endif;
endif;

i1=strindx(opts," MEAN=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tmean=stof(tmp) ; endif;
endif;



i1=strindx(opts," PRIME_START=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     i1=strsect(opts,i2+1,i3-(i2+1)) ;
     iprime=stof(i1);
     if scalmiss(iprime);
       call grbl2_error("Bad start-with-prime-number-index in G_RND_HALTON: " $+ i1 ,1);
       stop;
     endif;
     if iprime< 1 or iprime > rows(prim)  ;
       call grbl2_error("Out of range value of start-with-prime-number-index value (<0 or >" $+
                          g_ftos2(rows(prim),1) $+ ") in G_RND_HALTON: " $+ i1,1) ;
        stop;
     endif;
     prim=prim[iprime:rows(prim)];
endif;

@ now we are ready to go.... @

if rows(prim) lt nvs ;
    call grbl2_error("Not enough primes specified in G_RND_HALTON (must be >= to " $+ g_ftos2(nvs,1),1);
     stop;
endif;


@ todo is the  number of halton numbers needed for each variable,
  PLUS offset to get out of "correlated" starting values @

if ndiscard < prim[nvs] ; ndiscard=ndiscard+prim[nvs]; endif;
diddraws=(moredraws+1)*ndraws  ;
todo=(nrows*diddraws)+ndiscard ;

for hh (1,nvs,1) ;
  s=prim[hh];
  maxscram=(s^8 - ndiscard) ;

  dtype=typelist[hh];
  if s> rows(ascrs);
     call grbl2_error("Error: scramble does not support prime-value of " ,1); stop;
   endif;

  phi={0};
  scramble=ascrs[1:s]  ;

@ each loop will be adding rows(phi)* (s-1) values to X @

  i=1;

  do while 1==1 ;  @ keep adding fractionals until you have enough @
    x=phi;

@ possible improvement: if scrambled, dont bother with y.  @

    for j (1,s-1,1);
       y=phi+(j/s^i);
       if (rows(x)+rows(y)) > (todo+1120) ;
            x=x|y[1:(100+todo-rows(x)),.]    ;
       else;
          x=x|y;
       endif;
       if rows(x) ge todo ; break ; endif;
    endfor;

    phi=x;
    clear x,y ;

   if s>2 and i > 1  and (doscramble==1 and todo<maxscram) ;        @ no scramble for prime=2, or skip if first round @
      yscramble=zeros(rows(scramble)*s+ 100 ,1);
      iat=1;
      for j (1,rows(scramble),1);
        tmp=ascrs[1:s] $+ scramble[j] ;
        iat2=iat+rows(tmp);
        yscramble[iat:iat2-1]=tmp ;
        iat=iat2;
        if iat2 ge todo ; break ; endif;
      endfor;
      scramble=yscramble;
      clear yscramble ;
   endif;

   if rows(phi) ge todo ; break ; endif;

   i=i+1;
  endo;

 i1=ndiscard+1; i2=todo;


  phi=phi[i1:i2];  @drop initial zero and other, possibly "correlated", numbers.@

@ dont scramble if small primes and large sequences (algorithim only supports 8 levels; and besides, you
  dont really need to scramble in those cases @

  if (doscramble==1) and (todo<maxscram)  and (ndiscard<rows(scramble));
      scramble=scramble[(Ndiscard+1):todo];
      i1=sortindc(scramble);
      eps=seqa(1/(rows(phi)+1),1/rows(phi),rows(phi));
      phi=eps[i1];
  endif;


  if dorandom==1 ;      @ "randomize" the halton sequence? @
   arand=rndu(1,1);
   phi=phi+arand;
   phi=(phi.*(phi.<1)) + ((phi-1).*(phi.>1));
  endif;

@ phi is a 0-1 distributed "halton sequence  @

  if dtype==1 ;  @ normal random variable @
     hm1=cdfni(phi-(1e-10));		@ hack to avoid resolusion glitches @
     hm1=hm1.*(hm1 .<= 10)+10.*(hm1 .> 10);          /*truncates the normal dist. at 10 and -10*/
     hm1=hm1.*(hm1 .>= -10)-10.*(hm1 .< -10);
     hm1=reshape(hm1,nrows,diddraws);  @each row is an obs @
  elseif dtype==2 ;   @ extreme valued @
     hm1=-ln(-ln(phi));
     hm1=hm1.*(hm1 .<= 10)+10.*(hm1 .> 10);          /*truncates the extreme value dist. at 10 and -10*/
     hm1=hm1.*(hm1 .>= -10)-10.*(hm1 .< -10);
     hm1=reshape(hm1,nrows,diddraws);  @each row is an obs @

  elseif dtype==4 ;   @ triangle -1 - 1 @
     hm1= (sqrt(2.*phi)-1).*(phi .<= 0.5) + (1-sqrt(2.*(1-phi))) .* (phi .> 0.5); 
     hm1=reshape(hm1,nrows,diddraws);  @each row is an obs @

  elseif dtype==5;   @ uniform -1 - 1 @
     hm1=phi .*2 - 1;  
     hm1=reshape(hm1,nrows,diddraws);  @each row is an obs @

  else;		@ 3: uniform 0-1 @
     hm1=reshape(phi,nrows,diddraws);  @each row is an obs @
  endif;


  if moredraws > 0 ;	
    goo=seqa(1,moredraws,ndraws); @ extract every add'th column ;@
    hm1=hm1[.,goo];
  endif;

  if tsd ne 1 ; hm1=hm1*tsd ; endif;
  if tmean ne 0 ; hm1=hm1+tmean ; endif;
   
  if nvs==1 ; retp(hm1); endif ;

  y=areshape(hm1,nrows|ndraws|1);
  if hh==1 ;
      res=y;
  else ;
      res=aconcat(res,y,1);
  endif;
 
endfor ;        @ get next variable @
retp(res); 
endp ;


/********************/
@ returns a matrix of random numbers generated using modified latin hypercube sampling 

  eps_array=g_rnd_mlhs(nrows,ndraws,nvs,options);

  If nvs>1, this will be a 3d array.
  
  Eg: eps_array[112,43,2] would be the 2nd variable of the 43 draw for the 112th observation

  Options is either a blank string (a " " ), or a string containing a
  space delimited list of options.
  The options are:

     TYPE=atype                :   atype can be NORMAL, EXTREME, or UNIFORM
                                 to return a normal, extreme valued, or uniform random 
                                 variable (respectively). 
                                 By default, NORMAL is used.


    RANDOM=seed            : Set the seed value for the randomized sequence.


   SD=aval                 : Multiply random numbers by aval  (after conversion to normal, etc.).  
                               Thus, for NORMAL variates the returned values will have a 
                               standard deviation of SD. 
   MEAN=aval              :  Add aval to the generated values (after conversion to normal, etc.).
                               Thus, for NORMAL variates the returned values will have a 
                               mean of SD. 
                              Example: if TYPE=UNIFORM, SD=10, MEAN=-5; then a uniform distribution
                                       between -5 and 5 will be returned. In other words, the mean
                                       will be 0!


This is a implementation of the method described in see Hess, Train and Polak
   "On the use of Modified Latin Hypercube Sampling (MLHS) Method in the Estimation of a Mixed 
    Logit model for vehicle choice".

Basically, for the NDRAWS values (for one variable for one observation
 a) Generate a sequence of equally spaced values (0,1/NDRAWS,...,(NDRAWS-1)/NDRAWS)
 b) Add x, a uniform random value, to each number in the sequence (same value added to each);
    where 0 < x  < 1/NDRAWS
 c) Randomly reorder this sequence

A collection of NROWS*NVS of these sequences is used to form the NROWS x NDRAWS x NVS array.


@
proc g_rnd_mlhs(nrows,ndraws,nvs,opts);

local tlist,tsd,tmean,dtype,i1,i2,i3,atype;
local tmp,eps,gotseed;

local jv,e0,s0,jn,eps1,efoo,bvector,arand;


let tlist=NORMAL EXTREME UNIFORM  ;
tsd=1; tmean=0;
opts=" " $+  upper(opts)$+ " " ;
for i3 (1,rows(tlist),1);
  tmp=" " $+ tlist[i3] $+" ";
  i1=strindx(opts,tmp,1);
  if i1>0 ;
       opts=opts $+ " TYPE=" $+ tlist[i3] $+ " " ;
   endif;
endfor;   

dtype=1 ;  @default of normal rvs @
i1=strindx(opts," TYPE=",1);
if i1>0 ;
   i2=strindx(opts,"=",i1);
   i3=strindx(opts," ",i2);
   atype=strsect(opts,i2+1,i3-(i2+1)) ;
   dtype=missrv(indcv(atype,tlist),1);   @ default is normal @
endif;
 

i1=strindx(opts," SD=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tsd=stof(tmp) ; endif;
endif;

i1=strindx(opts," MEAN=",1);
if i1>0 ;
     i2=strindx(opts,"=",i1);
     i3=strindx(opts," ",i2);
     tmp=strsect(opts,i2+1,i3-(i2+1)) ;
     if tmp$/=""; tmean=stof(tmp) ; endif;
endif;

gotseed=miss(1,1);
i1=strindx(opts," RANDOM=",1);
if i1>0 ;
  i2=strindx(opts,"=",i1);
  i3=strindx(opts," ",i2);
  arand=strsect(opts,i2+1,i3-(i2+1)) ;
  if arand$/="" ;
     gotseed=stof(arand);
     if ismiss(gotseed) ;
             call grbl2_error(" Bad random number specified in rnd_mlhs: " $+arand,1) ; stop ; 
     endif;
     rndseed gotseed ;
  endif;    @ something after RANDOM= @
endif;



@ Generate the base vector -- 0.. 1/ndraws ... (ndraws-1)/ndraws  @
bvector=seqa(0,1/ndraws,ndraws);
eps={};
for jv (1,nvs,1);
e0=rndu(nrows,1).* (1/ndraws);
s0=rndu(ndraws,nrows);
for jn (1,nrows,1);
 eps1=bvector + e0[jn];
 
 efoo=sortind(s0[.,jn]);
 eps1=eps1[efoo];
 eps=eps|eps1;
endfor;
endfor;



@ sort into a matrix, or array @
if dtype==1 ;  @ normal random variable @
     eps=cdfni(eps-(1e-10));		@ hack to avoid resolusion glitches @
     eps=eps.*(eps .<= 10)+10.*(eps .> 10);          /*truncates the normal dist. at 10 and -10*/
     eps=eps.*(eps .>= -10)-10.*(eps .< -10);
elseif dtype==2 ;   @ extreme valued @
     eps=-ln(-ln(eps));
     eps=eps.*(eps .<= 10)+10.*(eps .> 10);          /*truncates the extreme value dist. at 10 and -10*/
     eps=eps.*(eps .>= -10)-10.*(eps .< -10);
endif;
if tsd ne 1 ; eps=eps*tsd ; endif;
if tmean ne 0 ; eps=eps+tmean ; endif;

if nvs==1 ;
  eps=reshape(eps,nrows,ndraws);
else ;
  eps=areshape(eps,nrows|ndraws|nvs);
endif;


retp(eps);
endp;






/*******************************/
@ procedure to randomly J items select from a Mx1 list, without replacement, given a 
  Mx1 vector of selection probabilties.
  A selection probablity of 0 means "never choose".
  Note that the list will be normalized.
  A Mx1 0/1 dummy vector is returned, with J values of 1
@

proc g_choose_noreplace(plist,ndo);
local chosen,a,aa,i1,culist;


if ndo>rows(plist);
  retp(ones(rows(plist),1));
endif;

chosen=plist*0;
for m (1,ndo,1);
  plist=plist./sumc(plist);
  culist=cumsumc(plist);
  a=rndu(1,1);
  aa=culist.>a;
  i1=indnv(1,aa);

@clearg hoy ; hoy=aa ;
clearg hoy2 ; hoy2=i1 ;
clearg hoy3 ; hoy3=a ;
clearg hoy4 ; hoy4=culist;@

  chosen[i1]=1;
  plist[i1]=0;
endfor;

retp(chosen);
endp;  


/*******************************************/
@ Compute great circle (as the crow flys) distance between two lat/long pairs @
proc gc_distance(alat1,along1,alat2,along2);
local orlatr,dslatr,lndifr,tmp,gcdist,arcdst;
 along1=abs(along1);
 along2=abs(along2);
  alat2=abs(alat2);
  alat1=abs(alat1);

  if ( ABS(alat1-alat2) Lt 0.001  and ABS(along1-along2) LT 0.001) ;  retp(0);   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 ;
retp(gcdist);
endp;
