//
//  neldermead.java
//   
//
//  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 neldermead
{
    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 simplex[][] = // [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++)
          simplex[i][FUNC] = func(simplex[i]); 

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

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

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

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

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

            if ((iter % 4*NDIMS) == 0)
            {
                if (simplex[ilo][FUNC] > best - TOL)
                  break; 
                best = simplex[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] += simplex[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] - simplex[ihi][j];
            double fr = func(r); 

            if ((flo <= fr) && (fr < fnhi))  // in zone: accept
            {
                for (int j=0; j<NDIMS; j++)
                  simplex[ihi][j] = r[j]; 
                simplex[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*simplex[ihi][j]; 
                double fe = func(e); 
                if (fe < fr)
                {
                   for (int j=0; j<NDIMS; j++)
                     simplex[ihi][j] = e[j]; 
                   simplex[ihi][FUNC] = fe; 
                   continue;
                }
                else
                {
                   for (int j=0; j<NDIMS; j++)
                     simplex[ihi][j] = r[j];
                   simplex[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*simplex[ihi][j]; 
                double fc = func(c); 
                if (fc <= fr)         
                {
                    for (int j=0; j<NDIMS; j++)
                      simplex[ihi][j] = c[j]; 
                    simplex[ihi][FUNC] = fc; 
                    continue; 
                }
                else   /////// contract
                {
                    for (int i=0; i<NPTS; i++)
                      if (i != ilo)
                      {
                          for (int j=0; j<NDIMS; j++)
                            simplex[i][j] = 0.5*simplex[ilo][j] + 0.5*simplex[i][j];
                          simplex[i][FUNC] = func(simplex[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*simplex[ihi][j]; 
                double fcc = func(cc); 
                if (fcc < fhi)
                {
                    for (int j=0; j<NDIMS; j++)
                      simplex[ihi][j] = cc[j]; 
                    simplex[ihi][FUNC] = fcc; 
                    continue; 
                }
                else    ///////// contract
                {
                    for (int i=0; i<NPTS; i++)
                      if (i != ilo)
                      {
                          for (int j=0; j<NDIMS; j++)
                            simplex[i][j] = 0.5*simplex[ilo][j] + 0.5*simplex[i][j];
                          simplex[i][FUNC] = func(simplex[i]); 
                      }
                }
            }
        }

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

    }

    static double func(double p[])
    {
        ncalls++; 
        return rosen(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; 
    }
}