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

abstract class f_xj
{
  // single function  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[]);
}

public class particle_swarm
{
int m;
int n;
double alpha,beta,beta0,teta,gamma;

// m number of population
// n number of variables for each individual 
double x[][];  //location for each population member
double y[];    //function value of each population member
double yxmin;  //current best minimum
double ygmin;  //global best minimum
double xmin[]; //current best x
double gmin[]; //global best x
double v[][];  //velocity vector
double xxmin[];
double xxmax[];
int ngeneration;  //number of generations in iteration
f_xj ff;

public particle_swarm(f_xj ffi,int mi, double xxmini[],double xxmaxi[],double alphai,double betai,int ngenerationi)
{
xxmin=xxmini;
xxmax=xxmaxi;
teta=0.5;
gamma=0.5;
ff=ffi;
//initial population set
n=xxmin.length;
m=mi;
alpha=alphai;
beta0=betai;
beta=beta0;
ngeneration=ngenerationi;
x=new double[m][n];
v=new double[m][n];
y=new double[m];
xmin=new double[n];
gmin=new double[n];

ngeneration=ngenerationi;
yxmin=1.0e50;
for(int i=0;i<m;i++)
{for(int j=0;j<n;j++)
  {x[i][j]=xxmin[j]+Math.random()*(xxmax[j]-xxmin[j]);
   v[i][j]=Math.random()*(xxmax[j]-xxmin[j])*0.01;
  }

  y[i]=ff.func(x[i]);
  if(y[i]<yxmin) 
    {yxmin=y[i];
     ygmin=y[i];
    for(int k=0;k<n;k++)
    {xmin[k]=x[i][k];
     gmin[k]=x[i][k];
    } //end for(int k=0;k<n;k++)  
    }//end if(y[i]<yxmin) 
}
//end of init
int t=0;

while(t<ngeneration)
{	
for(int i=0;i<m;i++)
{for(int j=0;j<n;j++)
  {x[i][j]=x[i][j]+v[i][j];
   if(x[i][j]<xxmin[j])  x[i][j]=xxmin[j]+Math.random()*(xxmax[j]-xxmin[j]);
   if(x[i][j]>xxmax[j])  x[i][j]=xxmax[j]-Math.random()*(xxmax[j]-xxmin[j]);
  }
   y[i]=ff.func(x[i]);
   if(y[i]<yxmin)  
   {yxmin=y[i];
    for(int k=0;k<n;k++)
    {xmin[k]=x[i][k];} 
   }
   if(y[i]<ygmin)  
   {ygmin=y[i];
   for(int k=0;k<n;k++)
   {gmin[k]=x[i][k];}
   }
   for(int j=0;j<n;j++)
   {v[i][j]=teta*v[i][j]+alpha*Math.random()*(gmin[j]-xmin[j])+beta*Math.random()*(xmin[j]-x[i][j]);}
}//end for int i=..	    
t++;
beta=beta0*Math.pow(gamma,0.001*t);
}

}


public String toString()
{ String s="ymin = "+ff.func(xmin)+"\nxmin = ";
  for(int k=0;k<n;k++)
  {s+=xmin[k]+" ";}
  s+="\ny global min = "+ff.func(gmin)+"\nxmin = ";
  for(int k=0;k<n;k++)
  {s+=gmin[k]+" ";}
  return s;
}

public String toString(double x[])
{ String s="ymin = "+ff.func(x)+"\nxmin = ";
  for(int k=0;k<n;k++)
  {s+=x[k]+" ";}
  return s;
}

public double[] xmin()
{ 
  return xmin;
}
public double[] gmin()
{ 
  return gmin;
}

public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance,int printlist)
{
	int i,j;
double x[][]=new double[a.length+1][a.length];
double p[][]=new double[a.length+1][a.length+1];
for(i=0;i<x.length;i++)
   {for(j=0;j<x[0].length;j++)
     {if(i==j){x[i][j]=a[j]+da[j];p[i][j]=x[i][j];}
      else    {x[i][j]=a[j];p[i][j]=x[i][j]; }
     }
     p[i][j] = fnelder.func(p[i]);
   }


   // Inlet variable definitions
   // fnelder : abstract multivariable function f(x)
   // x : independent variable set of n+1 simplex elements
   // maxiteration : maximum iteration number
   // tolerance : 
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// construct the starting simplex //////////////////
            //double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// calculate the first function values for the simplex ////////////////
 
            int iter=0;
            
            for (iter=1; iter<maxiteration; iter++)
            {
	            ///////////  define lo, nhi, hi (low high next_to_high //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                for (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 (i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                ////////// exit criteria //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    // calculate the avarage (including maximum value)
                    pavg=0;
                    for(i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(iter);
                       for (j=0; j<=NDIMS; j++)
                          { System.out.print(p[ilo][j]+" ");}
                       System.out.println("");
                    }
                    for(i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// calculate avarage without maximum value //////

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


                ///////// reflect ////////////////

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

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

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

                ///////////// shrink:

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

                if (fr >= fhi)   /// 
                {
                    double cc[] = new double[NDIMS];
                    for (j=0; j<NDIMS; j++)
                      cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                    double fcc = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// 
                    {
                        for (i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                    }
                }
            }
       return z;
}
      public static double[] nelder(f_xj fnelder,double a[],double da[],double tolerance)
      {return nelder(fnelder,a,da,500,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[])
      {return nelder(fnelder,a,da,500,1.0e-10,0);}
      
      public static double[] nelder(f_xj fnelder,double a[])
      {
	          double [] da=new double[a.length];
              for(int i=0;i<a.length;i++) da[i]=0.1*a[i];
              return nelder(fnelder,a,da); 
      }


}