/* All Search(Primitive Multi-dimensional Grid Search) */
new; cls;
xmin={1,0,0};
xmax={10,10,0.02};
d=0.01;
fn f(x)=sin(x[1])-cos(x[2])-x[3];
print Asearch(xmin,xmax,d,&f);


/*
**  Asearch.txt - ALL Search(Primitive Multi-dimensional Grid Search).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets values at the maximum of f(x) exactly between xmin and xmax
**              to the precision. Evaluates all points.
**
**  Format:     xstar=Asearch(xmin,xmax,d,&f);
**
**  Input:      xmin     vector,    dim x 1 vector of minimum values to search
**
**              xmax     vector,    dim x 1 vector of maximum values to search
**
**              d        scalar,    precision as decimal number
**
**              &f       pointer to a procedure of objective function f(x) to be maximized
**
**
**  Output:     xstar    vector,    dim x 1 vector of values at the maximum of the function
**
**  Notice:     Overflow-free. Light or normal version does not matter. Only time matters.
**              This program is very light. "Stupid Search."
*/

proc Asearch(xmin,xmax,d,&f);
   local dim,a,range,NN,x,zmax,zmaxj,xstar;
   local f:proc;
   dim=rows(xmin);
   a=1/d;
   range=xmax-xmin;
   NN=a.*range+1;
   x=zeros(dim,1);
   zmax=f(x);
   xstar=xmin;
   do while not x==(NN-1);
      x=recXindex(x,NN);
      zmaxj=f(xmin+d.*x);
      if zmax<zmaxj;
         xstar=xmin+d.*x;
         zmax=zmaxj;
      endif;
/*
      print (xmin+d.*x)';
*/
   endo;
   retp(xstar);
endp;


/*
**  recXindex.txt - Recursive Flexible N-base Dimension Indexing.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Gets flexible N-base index numbers for given dimension recursively 
**              without overflow problem.
**
**  Format:     x=recXindex(x,NN);
**
**  Input:      x        vector,    dim x 1 of current flexible N-base index vector
**
**              NN       vector,    dim x 1 of base numbers
**
**  Output:     x        vector,    dim x 1 of next N-base index vector
**
**  Notice:     You could start from x={0,0,0,0,0}; for example.
*/
proc recXindex(x,NN);
   local dim,j;
   dim=rows(x);
   x[dim]=x[dim]+1;
   j=dim;
   do until j==0;
      if x[j]==NN[j];
         x[j]=0;
         x[j-1]=x[j-1]+1;
      endif;
      j=j-1;
   endo;
   retp(x);
endp;