// optmn.java  optimize n parameters to achieve a minimum of a function 
//          the double precision function  y = f(double x[])         

// srchn    does a coarse global search to find initial values 
// provide below, function, f,  for finding minimum or
//                modify this java to pass an external function

class optmn
{
  double f(double x[]) // user defined function for test -1 to 1
  {
      return Math.sin(10.0*x[0]) + Math.sin(15.0*x[1]);
  }

  double optmn(                    // returned minimum value of function  
        int    n,                  // number of parameters 
        double x[],                // initial and returned optimum parameters 
        double yfinal,             // stop when y < yfinal 
        int    maxcnt,             // maximum number of steps 
        double tol,                // relative tolerance on step 
        double step[],             // start and returned step size 
        double xmin[],             // minimum value for each x 
        double xmax[])             // maximum value for each x 
        // double (*f)(double x[]) )  // optional passed by user 
  {
    double dx[] = new double[100]; // changeable step size in x
    double xl[] = new double[100]; // lower x value for test
    double xh[] = new double[100]; // higher x dalue for test
    double xb[] = new double[100]; // best x so far
    double xp[] = new double[100]; // next pivot point
    double dxv[] = new double[100];
    double dv;
    double vxl, vxh, v1 ; // initial, -dx, +dx  value of y 
    double vb, vp;        // best, pivot        
    int best, i, j, k, mstep;
    boolean noimp;
    boolean debug = true;
    int cnt = 0;

    v1 =f(x);                        // initial evaluation 
    noimp = true;                    // check for no improvement 

    while(true) // coming in with  cnt==0  gets return each iteration 
    {
      // v1 and x are available 
      vb = v1;      // best value of y 
      for(i=0; i<n; i++)
      {
        xb[i] = x[i]; // best x parameters 
        xp[i] = x[i]; // pivot x parameters 
      }

      // determine preferred direction, dxv's 
      for(i=0; i<n; i++)
      {
        dx[i] = (xmax[i]-xmin[i])* step[i]; // step size changes 
        if(debug) System.out.println("dx["+i+"]="+dx[i]);
        x[i] = xp[i]-dx[i];     // -dx low  
        if(x[i]<xmin[i]) x[i] = xmin[i];
        vxl=f(x);
        if(vxl<vb)
        {
          vb = vxl;      // best value of y 
          for(j=0; j<n; j++) xb[j] = x[j]; // best x parameters 
        }
        x[i] = xp[i]+dx[i];     // +dx high 
        if(x[i]>xmax[i]) x[i] = xmax[i];
        vxh=f(x);
        if(vxh<vb)
        {
          vb = vxh;      // best value of y 
          for(j=0; j<n; j++) xb[j] = x[j]; // best x parameters 
        }
        x[i] = xb[i];           // back to best 
        dxv[i] = 0.0;
        if     ((vxl > v1 && vxh < v1) || (vxh < v1 && vxh < vxl))
        {
          dxv[i] = +dx[i];
          noimp = false;
        }
        else if((vxl < v1 && vxh > v1) || (vxl < v1 && vxl < vxh))
        {
          dxv[i] = -dx[i];
          noimp = false;
        }
        else    step[i] = step[i]*0.7071;
      } // have each delta 
      if(debug)
      {
        System.out.println("cnt="+cnt+", v1="+v1+", vb="+vb);
        for(j=0; j<n; j++) 
           System.out.println("x["+j+"]="+x[j]+", raw dxv["+j+"]="+dxv[j]);
      }
      cnt = cnt+3;
      if(cnt > maxcnt) return v1;
         
      if(noimp) continue; // step size reduced, end while 

      dv = 0.0;
      for(i=0; i<n; i++)
      dv = dv + dxv[i]*dxv[i];
      dv = Math.sqrt(dv);
      if(dv<1.0e-6)  return v1; // too close to zero 
      for(i=0; i<n; i++)
        dxv[i] = dxv[i]/dv; // preferred direction 
      for(k=0; k<10; k++)
      {
        noimp = true;
        for(i=0; i<n; i++)
        {
          xl[i] = x[i] - dxv[i];
          if(xl[i]<xmin[i]) xl[i] = xmin[i];
          if(xl[i]>xmax[i]) xl[i] = xmax[i];
          xh[i] = x[i] + dxv[i];
          if(xh[i]>xmax[i]) xh[i] = xmax[i];
          if(xh[i]<xmin[i]) xh[i] = xmin[i];
        }
        vxl = f(xl);
        vxh = f(xh);
        if     ((vxl > v1 && vxh < v1) || (vxh < v1 && vxh < vxl))
        {
          noimp = false;
          v1 = vxh;
          vb = vxh;      // best value of y 
          for(j=0; j<n; j++)
          {
            dxv[j] = 2.0*dxv[j];
            x[j] = xh[j];
            xb[j] = x[j]; // best x parameters 
          }
        }
        else if((vxl < v1 && vxh > v1) || (vxl < v1 && vxl < vxh)) 
        {
          noimp = false;
          v1 = vxl;
          vb = vxl;      // best value of y 
          for(j=0; j<n; j++)
          {
            dxv[j] = 2.0*dxv[j];
            x[j] = xl[j];
            xb[j] = x[j]; // best x parameters 
          }  
        }
        else
        {
          for(j=0; j<n; j++)
          {
            dxv[j] = 0.5*dxv[j];
          }  
        } // end if 
        cnt = cnt+3;
        if(debug)
        {
          System.out.println("k="+k+", v1="+v1+", vxl="+vxl+", vxh="+vxh);
          for(j=0; j<n; j++) 
             System.out.println("x["+j+"]="+x[j]+", dxv["+j+"]="+dxv[j]);
        }
        if(cnt>maxcnt) return v1;
      } // loop on expand/contract search along preferred direction 

      // end this iteration, test for finish 
      v1 = vb;
      for(i=0; i<n; i++) x[i] = xb[i];

      if(debug) System.out.println("v1="+v1+", cnt=%"+cnt);
      mstep = 0;
      for(i=0; i<n; i++) // test for all steps below minimum 
      {
        if(step[i]<tol*(xmax[i]-xmin[i])) mstep++;
      }
      if(mstep==n) break; // all parameters close enough 
      if(v1 < yfinal)  break;
      if(cnt > maxcnt) break;
    } // end while 
    return v1;
  } // end optmn 

  double srchn(                         // returned minimum value of function  
           int    n,                    // number of parameters 
           double x[],                  // returned best found 

           int    maxcnt,               // maximum steps each dimension 
           double xmin[],               // minimum value for each x 
           double xmax[])               // maximum value for each x 
           // double (*f)(double x[]) ) // optional passed function 
  {
    double dx[] = new double[100];    // step size for each parameter 
    double step[] = new double[100];  // stem number for each parameter 
    double xbest[] = new double[100]; // best values of x 
    double ybest;      // best value of y 
    int i, j, k, mymax;
    boolean debug = false; // lots of printout
    double yy;
  
    // check input 
    if(maxcnt<3) maxcnt=3; // two ends and middle of each parameter 
    if(n>6 && maxcnt>3) maxcnt=3;
    if(n==6 && maxcnt>4) maxcnt=4;
    if(n==5 && maxcnt>10) maxcnt=10;
    if(n==4 && maxcnt>20) maxcnt=20;
    if(n==3 && maxcnt>200) maxcnt=200;
  
    mymax = 1;
    for(i=0; i<n; i++)
    {
      mymax = mymax*maxcnt;
      dx[i] = (xmax[i]-xmin[i])/(double)(maxcnt-1);
      x[i] = xmin[i]; // start of sweep 
      xbest[i] = x[i];
      step[i] = 0; // dx[i]*(double)step[i]  mymax 
      if(debug) 
        System.out.println("dx["+i+"]="+dx[i]+", x["+i+"]="+x[i]);
    }
    ybest = f(x); // first evaluation 
    if(debug)
      System.out.println("srchn evaluation at minimums, maxcnt="+maxcnt+
                         ", mymax="+mymax+", ybest="+ybest);
    if(debug)
    {
      System.out.println("index     minimum      maximum \n");
      for(k=0; k<n; k++) System.out.println(k+", "+xmin[k]+", "+xmax[k]);
    }
    j = 0;
    for(i=0; i<mymax-1; i++)
    {
      for(j=0; j<n; j++)
      {
        step[j]++;
        if(step[j]==maxcnt)
        {
          step[j]=0;
          x[j] = xmin[j]+dx[j]*(double)step[j];
        }
        else
        { 
          x[j] = xmin[j]+dx[j]*(double)step[j];
          break;
        }
      
      } // end for j 

      if(debug) for(k=0; k<n; k++) 
        System.out.println("i="+i+", j="+j+", x["+k+"]="+x[k]);
    
      yy = f(x);
      if(debug) System.out.println("ybest="+ybest+",  yy="+yy);
      if(yy<ybest)
      {
        for(k=0; k<n; k++) xbest[k] = x[k];
        ybest = yy;
      }
    }
    for(k=0; k<n; k++) x[k] = xbest[k];
    return ybest;
  } // end srchn 
} // end class optmn