/* Latin Mid Points */
new; cls;
nn=10;
dim=3;
x=LMP(nn,dim);

library pgraph;
graphset;
_plctrl=-1;
xyz(x[.,1],x[.,2],x[.,3]);



/*
**  LMP.txt - 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=LMP(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 is NOT randomized. If needed, see and use LHmp.txt instead.
**
**              This procedure uses procedure 'dimindex' inside.
*/
proc LMP(nn,dim);
   local indexvec,x;
   indexvec=dimindex(nn,dim);
/* representative mid-point values */
   x=(indexvec-0.5)/nn;
   retp(x);
endp;


/*
**  dimindex.txt - Dimension Indexing.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets index numbers for given dimension in a very easy way.
**
**  Format:     y=dimindex(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 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.
*/
proc dimindex(nn,dim);
   local x,b,n,z,y;
/* 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  */
   retp(y);
endp;