/* Multi-Candidates Latin Hypercube Search (Multi-dimensional Grid Search)                             */
/* Switching between deterministic and stochastic search with flexible focus ratios. Allow repetition. */
new; cls;
xmin={0,0,0};
xmax={10,10,10};
nn=3;
status={5,1,0};       /* 0 if deterministic, 1 if stochastic search, 2 or more if repeated */
p={0.2,0.1,0.1};      /* the last one is void */
fn f(x)=sin(x[1])-cos(x[2])-x[3];
print rsmLHSsearch(xmin,xmax,nn,status,p,&f);


/*
**  LHSCH04m.txt - Repetition Switching Multi-Candidates Latin Hypercube Search(Multi-dimensional Grid Search).    
**                 (Repetition in stochastic search means 'shake it inside each hypercube'.)
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets maximum values of function f(x) between xmin and xmax initially.
**
**  Format:     xstarstar=rsmLHSsearch(xmin,xmax,nn,status,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
**
**              status   vector,    0 if deterministic, 1 if stochastic search
**                                 (the number of rows of this is the number of steps to focus in)
**
**              p        vector,    focus ratios in each step (the last one is void)
**
**              &f       pointer to a procedure of objective function f(x) to be maximized
**
**
**  Output:  xstarstar   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.
**              If you have some trouble to run, start over GAUSS again.
*/
proc rsmLHSsearch(xmin,xmax,nn,status,p,&f);
   local dim,y1,y2,y3,range0,range,step,xmin0,zmax0,zmax,zmaxj,xstar,i,j,k,index,x,z,d,U,nsteps,c,xstarstar,r;
   local f:proc;
/* indexing */
   dim=rows(xmin);
   y1=dimindex0(nn,dim);      /* location index      */
   y2=dimindex0(nn+1,dim);    /* sub-location index  */
   y3=dimindex0(nn,dim)+1;    /* stochastic index beginning from 1 */
/* initial step settings */
   zmax0=-1e256;
   xmin0=xmin;
   range0=(xmax-xmin)/nn;
/* main */
   nsteps=rows(status);        /* # of steps to focus in */
   d=1/(nn);
   c=1;
   do while c<=nn^dim;
      xmin=xmin0'+y1[c,.].*range0'; xmin=xmin';
      range=range0;
      step=range/nn;
      zmax=f(xmin); 
      xstar=xmin;
      k=1;
      do while k<=nsteps;
         print;
         j=1;
         do while j<=nn^dim;
            if status[k]==0; /* deterministic search */
            /* regular grid including endpoints(allow overlaps) */
               x=xmin'+y1[j,.].*step'+y2.*step'/nn;                  
            /* 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 "candidate" c "/" nn^dim "   step" k "   " j "       /" nn^dim;
            else; /* stochastic search(status[k] times) */ 
               r=1;
               do while r<=status[k];                                /* 'status[k]' times */                      
                  /* LHS before scaling by parameter vector 'step' */
                  U=(d*y3-d*(y3-1)).*rndu(nn^dim,dim)+d*(y3-1);
                  /* [     location       ] + [ U(0,1) inside each cube]  */ 
                  x=xmin'+(y3[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 "candidate" c "/" nn^dim "   step" k "   " j "       /" nn^dim;
                  r=r+1;
               endo;
            endif;            
            j=j+1;
         endo;
      /* update step length for next step to focus in */
         range=range*p[k];                               /* focus ratios are changing */
         xmin=xstar-1/2*range;
         step=range/nn;
         k=k+1;
      endo;
      /* update if larger zmax is found for each candidate */
      if zmax0<zmax;
         xstarstar=xstar;
         zmax0=zmax;
      endif;
      c=c+1;
   endo;
   retp(xstarstar);
endp;


/*
**  dimindex0.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=dimindex0(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (0,1,2,3,...,nn-1) beginning from 0
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     y        matrix ,   (nn^dim) x (dim) of resulting index matrix
**
*/
proc dimindex0(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=zeros(1,dim)|y;           /* insert 0's at the 1st row  */
   retp(y);
endp;