/* Improved Latin Hypercube */
new; cls;
nn=10;
dim=3;
x=iLHS(nn,dim);

library pgraph;
graphset;
_plctrl=-1;
pqgwin auto;
call hist(x[.,1],19);
call hist(x[.,2],19);
call hist(x[.,3],19);
xy(x[.,1],x[.,2]);
xy(x[.,2],x[.,3]);
xy(x[.,1],x[.,3]);
xyz(x[.,1],x[.,2],x[.,3]);


/*
**  iLHS.txt - Improved Latin Hypercube U(0,1).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Improved Latin Hypercube numbers in a very easy way. New algorithm.
**              Basically, it shows a little vibration around iLMP pattern in each dimension.
**
**  Format:     x=iLHS(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     x        matrix ,   (nn^dim) x (dim) of resulting U(0,1) numbers
**
**  Notice:     It takes lots of memory to run. Light version of GAUSS will not work in most cases.
**
**              This procedure uses procedure 'dimindex1' inside along with 'rerankindx'.
*/
proc iLHS(nn,dim);
   local indexvec,d,U,i,index;
   indexvec=dimindex1(nn,dim);
/* divide between 0 and 1 into (nn^dim) segments and sample from each segment */
   d=1/(nn^dim);
   U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   do while NOT(U>0);                                                      /* check if U contains 0's */
      U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   endo;
/* randomize the order as a row */
   x=zeros(nn^dim,dim);
   i=1;
   do while i<=nn^dim-1;
      index=ceil((nn^dim-(i-1))*rndu(1,1));
      x[i,.]=U[index,.];
      if index==1;
         U=U[2:rows(U),.];
      elseif index==rows(U);
         U=U[1:rows(U)-1,.];
      else;
         U=U[1:(index-1) (index+1):rows(U),.];
      endif;
      i=i+1;
   endo;
   x[nn^dim,.]=U;
   retp(x);
endp;


/*
**  dimindex1.txt - Dimension Indexing by re-partitioning the same index numbers
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets improved index numbers for given dimension in a very easy way.
**
**  Format:     x=dimindex1(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     x        matrix ,   (nn^dim) x (dim) of resulting index matrix
**
**  Notice:     It takes lots of memory to run. Light version of GAUSS will not work in most cases.
**              Here, index starts from 1 as in GAUSS.
**
**  This algorithm can be easily implemented in EXCEL, Stata, etc... No more hustle on quasi-random sequences.
*/
proc dimindex1(nn,dim);
   local x,b,n,z,y,j;
/* convert x(base 10) to y(base b) */
   x=seqa(1,1,nn^dim-1); 
   b=nn;
   n=maxc(log(x)/log(b))+1;
   z=reshape(b,n,rows(x));
   y=rev(recserrc(x,z))';
/* adjustments */
   y=y+1;                     /* shift all elements by 1    */
   y=ones(1,dim)|y;           /* insert 1's at the 1st row  */
/* re-indexing each number */
   x=zeros(nn^dim,dim);
   j=1;
   do while j<=dim;
      x[.,j]=rerankindx(y[.,j]);
      j=j+1;
   endo;
   retp(x);
endp;

proc rerankindx(x);
   local k,y,i,j;
   k=1;
   y=zeros(rows(x),cols(x));
   j=minc(x);
   do while j<=maxc(x);
      i=1;
      do while i<=rows(x);
         if x[i]==j;
            y[i]=k;
            k=k+1;
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   retp(y);
endp;