/* Latin Hypercube Search (Multi-dimensional Grid Search) stochastic prototype */
/* Notice:    This set of parameters is only for GAUSS light. You need larger  */
/*            parameters with normal version of GAUSS. Results could vary.     */
new; cls;
xmin={0,0,0};
xmax={10,10,10};
nn=2;            /* more than 10 in normal version of GAUSS required */
nsteps=2;        /* this is not enough */
p=0.2;
fn f(x)=sin(x[1])-cos(x[2])-x[3];
print LHSsearch(xmin,xmax,nn,nsteps,p,&f);


/*
**  LHSCH02.txt - Latin Hypercube Sampling Search(Multi-dimensional Grid Search).
**                Stochastic Sampling version of prototype.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets maximum values of function f(x) between xmin and xmax initially.
**
**  Format:     xstar=LHSsearch(xmin,xmax,nn,nsteps,p,&f);
**
**  Input:      xmin     vector,    nn x 1 vector of initial minimum values to search
**
**              xmax     vector,    nn x 1 vector of initial maximum values to search
**
**              nn       scalar,    number of segments
**
**              nstep    scalar,    number of steps to focus in
**
**              p        scalar,    focus ratio in each step
**
**              &f       pointer to a procedure of objective function f(x) to be maximized
**
**
**  Output:     xstar    vector,    nn x 1 vector of maximum values
**
**  Notice:     It takes lots of memory to run. Light version of GAUSS will not work in most cases.
**              Partition parameter nn should be large enogh. Partial use of this algorithm would
**              be very powerful at the initial stages.
**
**              If you have some trouble to run, start over GAUSS again.
*/
proc LHSsearch(xmin,xmax,nn,nsteps,p,&f);
   local dim,y,range,step,zmax,zmaxj,xstar,i,j,k,index,x,z,d,U;
   local f:proc;
/* indexing */
   dim=rows(xmin);
   y=dimindex(nn,dim);
/* initial step settings */
   range=xmax-xmin;
   step=range/nn;
   zmax=f(xmin); 
   xstar=xmin;
/* main */
   d=1/(nn);
   k=1;
   do while k<=nsteps;
      print;
      j=1;
      do while j<=nn^dim;
         /* LHS before scaling by parameter vector 'step' */
            U=(d*y-d*(y-1)).*rndu(nn^dim,dim)+d*(y-1);
         /*    [    location       ] + [ U(0,1) inside each cube]  */
            x=xmin'+(y[j,.]-1).*step'+U.*step';
         /* calculation f(x) in each hypercube */
            z=zeros(nn^dim,1);
            i=1;
            do while i<=nn^dim;
               z[i]=f(x[i,.]');
               i=i+1;
            endo;
            zmaxj=maxc(z);
         /* uptate xstar if larger solution zmaxj is found */
            if zmax<zmaxj;
               index=maxindc(z);
               xstar=x[index,.]';
               zmax=zmaxj;
            endif;
                  /* just counting */
                  if j%50==0;
                     call sleep(1);
                  endif;
                  print/lz "step" k "   " j "       /" nn^dim;
         j=j+1;
      endo;
   /* update step length for next step to focus in */
      range=range*p;
      xmin=xstar-1/2*range;
      step=range/nn;
      k=k+1;
   endo;
   retp(xstar);
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;