/* 2D Mid-point version of Latin Hypercube (nn^2 x 2) */
new; cls; 
nn=50;
x=LH2mp(nn);

library pgraph;
graphset;
_plctrl=-1;
xtics(0,1,0.1,0);
ytics(0,1,0.1,0);
xy(x[.,1],x[.,2]);



/*
**  LH2mp.txt - 2D Latin Mid Points U(0,1).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Latin Mid-point numbers in a very easy way. New algorithm.
**
**  Format:     x=LH2mp(nn);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**  Output:     x        matrix ,   (nn^2) x 2 of resulting U(0,1) numbers
**
**  Notice:     It takes lots of memory to run.
*/
proc LH2mp(nn);
   local x,b,n,z,y,indexvec,i,index,U;
/* convert x(base 10) to y(base b) */
   x=seqa(1,1,nn^2-1); 
   b=nn;
   n=maxc(log(x)/log(b))+1;
   z=reshape(b,n,rows(x));
   y=rev(recserrc(x,z))';
   indexvec=y+1;                     /* shift all elements by 1    */
   indexvec=ones(1,2)|indexvec;      /* insert 1's at the 1st row  */
/* representative mid-point values */
   U=(indexvec-0.5)/nn;
/* randomize the order as a row */
   x=zeros(nn^2,2);
   i=1;
   do while i<=nn^2-1;
      index=ceil((nn^2-(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^2,.]=U;
   retp(x);
endp;