/* Latin Hypercube of regular definition up to dim=1,2,3,..,M */
new; cls;
nn=3000;
dim=3;
x=LHss1(nn,dim);

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



proc LHss1(nn,dim);
   local x,i;
   x=U01ss(nn);
   if dim/=1;
      i=1;
      do while i<=dim-1;
         x=x~U01ss(nn);
         i=i+1;
      endo;
   endif;
   retp(x);
endp;

proc U01ss(nn);
   local d,i,U,x,index;
/* divide between 0 and 1 into n segments and sample from each segment */
   d=1/nn;
   U=zeros(nn,1);
   i=1;
   do while i<=nn;
      U[i]=(d*i-d*(i-1))*rndu(1,1)+d*(i-1);
      i=i+1;
   endo;
/* check if the first element is zero or not */
   do while U[1]==0;
      U[1]=d*rndu(1,1);
   endo;
/* randomize the order */
   x=zeros(nn,1);
   i=1;
   do while i<=nn-1;
      index=ceil((nn-(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]=U;
   retp(x);
endp;