/* Scrambled Improved LHS-Sobol Sequence */
/* Notice:Very slow. Wait for a while... */
new; cls;
nn=5;
dim=3;
x=siLHSSobol(nn,dim);

library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
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]);



/*
**  siLHSHalton.txt - Scrambled Improved LHS-Sobol Sequence.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Scrambled Improved LHS-Sobol numbers in a very easy way.
**
**  Format:     xx=siLHSSobol(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     xx       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 'dimindex3' with 'rerankindx' and 'Sobol' inside.
*/
proc siLHSSobol(nn,dim);
   local y,x,bb,SO,xx,j,i;
   {y,x,bb}=dimindex3(nn,dim); 
/*  y: location parameter of Hypercube            */
/*  x: ordering parameter in each quasi sequence  */
/* bb: dimension parameter in each quasi sequence */
   SO=Sobol((nn^dim),dim);
   xx=zeros(nn^dim,dim);
   j=1;
   do while j<=dim;
      i=1;
      do while i<=nn^dim;
         if rndu(1,1)<0.5;
            xx[i,j]=SO[x[i,j],bb[i,j]]/nn+(y[i,j]-1)/nn;   /* Now, use bb as random dimension index */
            /*      [    in each cell   ]+[ location  ]    */
         else;
             xx[i,j]=(1/nn-SO[x[i,j],bb[i,j]]/nn)+(y[i,j]-1)/nn; 
            /*      [         in each cell      ]+[ location  ]    */
         endif;
         i=i+1;
      endo;
      j=j+1;
   endo;
   retp(xx);
endp;

proc Sobol(nn,dim);
   local x,lastq,Sseed,maxdim,seed,v,p,m,i,j,k,l,jj,newv,inc,rec,tempseed;
   if dim<2 or dim>40;
      errorlog "ERROR: Dimension number must be an integer where 2 <= dim <= 40.";
      retp(".");
   endif;
   x=zeros(nn,dim);
   lastq=zeros(1,dim);
   Sseed=0;
   maxdim=40;
   seed=2;  /* exclude seed=1 corresponding to zero vector at the beginning */
   do while seed<=nn+1;
      let v[40,30]=
          1  0  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  0  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  0  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5  0  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  1  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  7  0  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  9  9  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7 13  3  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5 11 27  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  3  5  1 15  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7  3 29  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  7 21  0  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  9 23 37  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  5 19 33  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3 13 11  7  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7 13 25  5  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  5 11  7 11  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  3 13 39  0  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1 15 17 63 13  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  1  5  5  1 27 33  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  3 25 17 115 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3 15 29 15 41  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  7  3 23 79  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  9 31 29 17  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5 13 11  3 29  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  9  5 21 119 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  1 23 13 75  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3 11 27 31 73  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7  7 19 25 105 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

          1  3  5  5 21  9  7  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1 15  5 49 59  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  1  1  1 33 65  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  5 15 17 19 21  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  7 11 13 29  3  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  7  5  7 11 113 0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  5  3 15 19 61  0  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  1  1  9 27 89  7  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  1  3  7 31 15 45 23  0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
          1  3  3  9  9 25 107 39 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
         ;
      p={1,   3,   7,  11,  13,  19,  25,  37,  59,  47,
        61,  55,  41,  67,  97,  91, 109, 103, 115, 131,
        193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
        213, 191, 253, 203, 211, 239, 247, 285, 369, 299};
      i=1;
      v[i,1:30]=ones(1,30);
      i=2;
      do while i<=maxdim;
         m=sumc(p[i].>=(2^seqa(1,1,9)));
         j=p[i];
         inc=zeros(m,1);
         k=m;
         do while k>=1;
            jj=floor(j/2);
            inc[k]=(j/=2*jj);
            j=jj;
            k=k-1;
         endo;
         j=m+1;
         do while j<=30;
            newv=v[i,j-m];
            l=1;
            k=1;
            do while k<=m;
               l=2*l;
               if inc[k];
                  newv=eor(newv,l*v[i,j-k]);
               endif;
               k=k+1;
            endo;
            v[i,j]=newv;
            j=j+1;
         endo;
         i=i+1;
      endo;
      l=1;
      j=30-1;
      do while j>=1;
         l=2*l;
         v[1:dim,j]=v[1:dim,j]*l;
         j=j-1;
      endo;
      rec=1/(2*l);
      if seed==Sseed+1;
         l=low0(seed);
      elseif seed<=Sseed;
         Sseed=0;
         l=1;
         lastq[1:dim]=0;
         tempseed=Sseed;
         do while tempseed<=seed-1;
            l=low0(tempseed);
            i=1;
            do while i<=dim;
               lastq[i]=eor(lastq[i],v[i,l]);
               i=i+1;
            endo;
            tempseed=tempseed+1;
         endo;
         l=low0(seed);
      elseif Sseed+1<seed;
         tempseed=Sseed+1;
         do while tempseed<=seed-1;
            l=low0(tempseed);
            i=1;
            do while i<=dim;
               lastq[i]=eor(lastq[i],v[i,l]);
               i=i+1;
            endo;
            tempseed=tempseed+1;
         endo;
         l=low0(seed);
      endif;
      i=1;
      do while i<=dim;
         x[seed-1,i]=lastq[i]*rec;
         lastq[i]=eor(lastq[i],v[i,l]);
         i=i+1;
      endo;
      Sseed=seed;
      seed=seed+1;
   endo;
   retp(x);
endp;


proc low0(x);
   local n,z,y;
   if x==0;
      retp(0);
   else;
      n=maxc(log(x)/log(2))+1;
      z=reshape(2,n,rows(x));
      y=rev(recserrc(x,z));
      y=0|y;
      y=rev(y);
      y=seqa(1,1,rows(y)).*(y.==0);
      y=delif(y,y.==0);
      retp(y[1]);
   endif;
endp;

proc eor(a,b);
   local a2,b2,dim,y,x;
   a2=dectobin(a);
   b2=dectobin(b);
   dim=maxc(rows(a2)|rows(b2));
   if rows(a2)>rows(b2);
      b2=zeros(rows(a2)-rows(b2),1)|b2;
   elseif rows(a2)<rows(b2);
      a2=zeros(rows(b2)-rows(a2),1)|a2;
   endif;
   y=a2+b2;
   y=y.*(y./=2);
   x=bintodec(y);
   retp(x);
endp;

proc dectobin(x);           /* dec to bin */
   local n,z,y;
   if x==0;
      y=0;
   else;
      n=maxc(log(x)/log(2))+1;
      z=reshape(2,n,rows(x));
      y=rev(recserrc(x,z)); 
   endif;
   retp(y);
endp;

proc bintodec(y);           /* bin to dec */
   if rows(y)==1;
      retp(y);
   else;
      retp(polyeval(2,y));
   endif;
endp;


/*
**  dimindex3.txt - Threefold Swapped Dimension Indexing
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets threefold swapped index numbers for given dimension in a very easy way.
**
**  Format:     {y,x,bb}=dimindex3(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     y        matrix,   (nn^dim) x (dim) of 1st resulting index matrix
**
**              x        matrix,   (nn^dim) x (dim) of 2nd resulting index matrix for each group
**
**              bb       matrix    (nn^dim) x (dim) of 3rd swapped column numbers for each group   
**
**  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.
**
**              I on purpose program this in a primitive way. We could just build up indeces if
**              you use the latest version of  GAUSS or that kind of software.
**
**  This algorithm can be easily implemented in EXCEL, Stata, etc...
*/
proc(3)=dimindex3(nn,dim);
   local x,b,n,z,y,j,bb,i,rindex;
/* 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;
/* swapping column numbers by random ordering */
   bb=zeros(nn^dim,dim);
   i=1;
   do while i<=((nn^dim)/nn);
      rindex=rankindx(rndu(dim,1),1);             /* random order of column numbers */
      bb[nn*(i-1)+1:nn*(i-1)+nn,.]=rindex'.*ones(nn,dim);
      i=i+1;
   endo;
   retp(y,x,bb);
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;