
//
//  neldermead.java
abstract class f_xi
{
  // multifunction multi independent variable
  // a single value is returned indiced to equation_ref
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x,1) returns the value of f[1]
  // func(x,0) returns the value of f[0]

  abstract double func(double x[],int equation_ref);
  // single function multi independent variable
  // a single value is returned indiced to equation_ref
  // example f=x[0]+sin(x[1])
  abstract double func(double x[],int equation_ref);
}

class f extends f_xi
{
  double func(double p[])
 // Rosenbrock banana, two dimensions
 {
     double r0 = 10.0*(p[1] - SQR(p[0]));
     double r1 = 1.0 - p[0];
     return r0*r0 + r1*r1;
 }

}
//
//
//  simplex minimizer
//  Nelder & Mead 1965 Computer J, v.7, 308-313.
//  Lagarias et al 1998 SIAM J.Optim. p.112
//  M.Lampton UCB SSL 2004
//
//  I use an augmented simplex that keeps funcval with each vertex.

import java.io.*;
import java.util.*;
import java.text.DecimalFormat;


public class nelder
{
    static final int NDIMS = 2;
    static final int NPTS = 3;
    static final int FUNC = NDIMS;
    static final int MAXITER = 200;
    static int ncalls = 0;
    static final double TOL = 1E-6;

    public static void main (String args[])
    {
        descend();
        System.out.println("\nThat's all folks");
    }

    static void descend()
    {

        ////// set up the starting simplex //////////////////

        double p[][] = // [row][col] = [whichvx][coord,FUNC]
        {{ 10.0, 10.0, 0.0},
         { 10.0, 11.0, 0.0},
         { 11.0, 10.0, 0.0}};

        double best = 1E99;

        //////////////// initialize the funcvals ////////////////

        for (int i=0; i<NPTS; i++)
          p[i][FUNC] = func(p[i]);

        System.out.println("ncalls = "+fwi(ncalls,6));
        int iter=0;

        for (iter=1; iter<MAXITER; iter++)
        {
            /////////// identify lo, nhi, hi points //////////////

            double flo = p[0][FUNC];
            double fhi = flo;
            int  ilo=0, ihi=0, inhi = -1; // -1 means missing
            for (int i=1; i<NPTS; i++)
            {
               if (p[i][FUNC] < flo)
                 {flo=p[i][FUNC]; ilo=i;}
               if (p[i][FUNC] > fhi)
                 {fhi=p[i][FUNC]; ihi=i;}
            }
            double fnhi = flo;
            inhi = ilo;
            for (int i=0; i<NPTS; i++)
              if ((i != ihi) && (p[i][FUNC] > fnhi))
                {fnhi=p[i][FUNC]; inhi=i;}

            for (int j=0; j<=NDIMS; j++)
              System.out.print(fwd(p[ilo][j], 18, 9));
            System.out.println();

            ////////// exit criterion //////////////

            if ((iter % 4*NDIMS) == 0)
            {
                if (p[ilo][FUNC] > best - TOL)
                  break;
                best = p[ilo][FUNC];
            }

            ///// compute ave[] vector excluding highest vertex //////

            double ave[] = new double[NDIMS];
            for (int j=0; j<NDIMS; j++)
              ave[j] = 0;
            for (int i=0; i<NPTS; i++)
              if (i != ihi)
                for (int j=0; j<NDIMS; j++)
                   ave[j] += p[i][j];
            for (int j=0; j<NDIMS; j++)
               ave[j] /= (NPTS-1);


            ///////// try reflect ////////////////

            double r[] = new double[NDIMS];
            for (int j=0; j<NDIMS; j++)
              r[j] = 2*ave[j] - p[ihi][j];
            double fr = func(r);

            if ((flo <= fr) && (fr < fnhi))  // in zone: accept
            {
                for (int j=0; j<NDIMS; j++)
                  p[ihi][j] = r[j];
                p[ihi][FUNC] = fr;
                continue;
            }

            if (fr < flo)  //// below zone; try expand, else accept
            {
                double e[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  e[j] = 3*ave[j] - 2*p[ihi][j];
                double fe = func(e);
                if (fe < fr)
                {
                   for (int j=0; j<NDIMS; j++)
                     p[ihi][j] = e[j];
                   p[ihi][FUNC] = fe;
                   continue;
                }
                else
                {
                   for (int j=0; j<NDIMS; j++)
                     p[ihi][j] = r[j];
                   p[ihi][FUNC] = fr;
                   continue;
                }
            }

            ///////////// above midzone, try contractions:

            if (fr < fhi)  /// try outside contraction
            {
                double c[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  c[j] = 1.5*ave[j] - 0.5*p[ihi][j];
                double fc = func(c);
                if (fc <= fr)
                {
                    for (int j=0; j<NDIMS; j++)
                      p[ihi][j] = c[j];
                    p[ihi][FUNC] = fc;
                    continue;
                }
                else   /////// contract
                {
                    for (int i=0; i<NPTS; i++)
                      if (i != ilo)
                      {
                          for (int j=0; j<NDIMS; j++)
                            p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                          p[i][FUNC] = func(p[i]);
                      }
                    continue;
                }
            }

            if (fr >= fhi)   /// over the top; try inside contraction
            {
                double cc[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                double fcc = func(cc);
                if (fcc < fhi)
                {
                    for (int j=0; j<NDIMS; j++)
                      p[ihi][j] = cc[j];
                    p[ihi][FUNC] = fcc;
                    continue;
                }
                else    ///////// contract
                {
                    for (int i=0; i<NPTS; i++)
                      if (i != ilo)
                      {
                          for (int j=0; j<NDIMS; j++)
                            p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                          p[i][FUNC] = func(p[i]);
                      }
                }
            }
        }

        System.out.println("ncalls, iters, Best ="+fwi(ncalls,6)+fwi(iter,6)+fwd(best,16,9));

    }

    static double func(double p[])
    {
        ncalls++;
        return parab(p);
    }

    static double rosen(double p[])
    // Rosenbrock banana, two dimensions
    {
        double r0 = 10.0*(p[1] - SQR(p[0]));
        double r1 = 1.0 - p[0];
        return SQR(r0) + SQR(r1);
    }

    static double parab(double p[])
    // simple paraboloid
    {
        return SQR(p[0]-2) + SQR(p[1]-20);
    }


    /////////////////////////////////utilities ////////////////////

    static double SQR(double x)
    {
        return x*x;
    }

    static String fwi(int n, int w)
    // converts an int to a string with given width.
    {
        String s = Integer.toString(n);
        while (s.length() < w)
           s = " " + s;
        return s;
    }


    static String fwd(double x, int w, int d)
    // converts a double to a string with given width and decimals.
    {
        java.text.DecimalFormat df = new DecimalFormat();
        df.setMaximumFractionDigits(d);
        df.setMinimumFractionDigits(d);
        df.setGroupingUsed(false);
        String s = df.format(x);
        while (s.length() < w)
          s = " " + s;
        if (s.length() > w)
        {
            s = "";
            for (int i=0; i<w; i++)
              s = s + "-";
        }
        return s;
    }
}
