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 firefly
{
double gamma,alpha,beta;
double alpha0,alphan;
int m;
int n;
// m number of fireflies
// n number of variables for each firefly 
// I0 light intensity at source
// gamma absorbtion coefficient
// alpha size of the random step
// beta attractiveness
double x[][];
double I0[];
double ymin;
double xmin[];
int imin;
int ngeneration;
public firefly(f_xj ff,int mi,int ni, double xxmin[],double xxmax[],double gammai,double alpha0i,double alphani,int ngenerationi)
{
// initial population set
n=ni;
m=mi;
x=new double[m+1][n];
I0=new double[m+1];
xmin=new double[n];
gamma=gammai;
alpha0=alpha0i;
alphan=alphani;
ngeneration=ngenerationi;
for(int j=0;j<n;j++)
  {xmin[j]=xxmin[j]+Math.random()*(xxmax[j]-xxmin[j]);}
ymin=ff.func(xmin);
  //end of init
int l=0;
double d=0;
double t=0;
while(l<ngeneration)
{ 
// find minimum of population
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]);}
  I0[i]=ff.func(x[i]);
  if(I0[i]<ymin) 
    {ymin=I0[i];imin=i;
    for(int k=0;k<n;k++)
    {xmin[k]=x[i][k];}   
    };
}
for(int j=0;j<n;j++)
  {x[m][j]=xmin[j];}
  I0[m]=ymin;
l++;
double fi=0,fj=0;
for(int i=0;i<=m;i++)
{ //fi=ff.func(x[i]);
  for(int j=0;j<n;j++)
  { 
    //fj=ff.func(x[j]);
    if(I0[i]<I0[j])
    {
    beta=attractiveness(i,j);
    for(int k=0;k<n;k++)
    {   //The best one will remain the rest will change
	    if(i!=imin)
        { x[i][k]=(1.0-beta)*x[i][k]+beta*x[j][k]+alpha*(Math.random()-0.5);} }    
    }
    else
    {
	    for(int k=0;k<n;k++)
        {  if(i!=imin)
	       { x[i][k]=x[i][k]+alpha*(Math.random()-0.5); }   
	    }
    }
  }

}

    //alpha will change by time 
    //(shrink to a smaller value)  
    alpha=alphan+(alpha0-alphan)*Math.exp(-t);
    t=0.1*l;
}
}

public double distance(int i, int j)
{
double d=0;
for(int k=0;k<n;k++)
{d+=(x[i][k]-x[j][k])*(x[i][k]-x[j][k]);}
d=Math.sqrt(d);
return d;	
}

public double attractiveness(int i,int j)
{   double d=distance(i,j);
	return I0[i]/(1+gamma*d*d);
	//return I0[i]*Math.exp(-gamma*d*d);
}

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

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); 
      }


}