/* Threefold Dimension Indexing */
new; cls;
nn=3;
dim=3;
{y,x,bb}=dimindex3(nn,dim);
print/rz "1st Index y=" y;
print/rz "2nd Index x=" x;
print/rz "3rd Index bb=" bb;


/*
**  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;