import java.io.*;
import javax.swing.*;

      class fxa extends f_xj
      {  
	  double xi[][]; // m independent variable data set
      double yi[]; //dependent variable data set
      double a[];  //fit function coefficient set
      int nn;//number of data
      int m; //number of idependent variables xi
      double xmax;
  	  public fxa(String filename,double ia[])
      {
        //read the data to curvefit
        //last column is treated as independent data set
        //get the data file and initial fit coefficient when class is defined
        int n=ia.length;
        a=new double[n];
        seta(ia);
        xmax=-1.0e99;
        double b[][]=Text.readDoubleT(filename);
        m=b.length-1;
        yi=b[m];
        nn=yi.length;
        //System.out.println("m="+m);
        xi=new double[m][nn];
        for(int i=0;i<nn;i++) {
	    xi[0][i]=b[0][i];if(xi[0][i]>xmax) {xmax=xi[0][i];}
        xi[1][i]=b[1][i];
	    yi[i]=b[m][i];
        }
        //System.out.println("b=\n"+Matrix.toString(Matrix.T(b)));
        System.out.println("xi=\n"+Matrix.toString(xi));
        System.out.println("yi=\n"+Matrix.toString(yi));      
      nn=yi.length;
      a=ia;
      }
      
      public void seta(double ia[])
      {
	      //assign new fit coefficient set 
	      for(int ii=0;ii<nn;ii++)
              a[ii]=ia[ii];
	  }

	  public double[] geta()
      {
	     // return fit coefficient set
	     return a;
	  }	  
       double Ps(double x[],double ai[])
      {
      // sample specific function m dimensional independent variable
      // independent variables : x[0],x[1]...x[m-1]
      // x[0] t
      // x[1] kalınlık
      // ai[0] a
      // 
      a=ai;
      double t=x[0];
      double h=x[1];
      double a=ai[0]+ai[1]*h;
      double k=ai[2]+ai[3]*h;
      double pi=a*Math.exp(-k*t);
      return pi; //kPa
      }
	        
      public double func(double ai[])
      {
      double ff=0;
      double w;
      double yy;
      double xx[]=new double[xi.length];
       for(int i=0;i<nn;i++)
	    {   for(int j=0;j<m;j++) {xx[j]=xi[j][i];} 
	    		    w=Ps(xx,ai);
	    yy=yi[i];
        //System.out.println("x="+xi[i]+"w="+w+"yy="+yy+"ff="+ff);
	    w-=yy; 
        ff+=w*w;
        //System.out.println("x="+xi[i]+"w="+w+"yy="+yy+"ff="+ff); 
	    }
      return ff;
	  }
}

public class NA126
{//2. Nelder ve Mead method
//________________________________________________________________
public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance,int printlist)
{
double x[][]=new double[a.length+1][a.length];
for(int i=0;i<x.length;i++)
   {for(int j=0;j<x[0].length;j++)
     {if(i==j){x[i][j]=a[j]+da[j];}
      else    {x[i][j]=a[j]; }
     }
   }
   //  Nelder mead multidimensional simplex minimization method
   //  Reference : Nelder & Mead 1965 Computer J, v.7, 308-313.
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// create 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 (int i=0; i<NPTS; i++)
            {
              for (int j=0; j<NDIMS; j++)
                  {p[i][j] = x[i][j];}
              p[i][NDIMS] = fnelder.func(p[i]);
            }
            int iter=0;
            for (iter=1; iter<maxiteration; iter++)
            {
                /////////// define lowest, next to highest and highest points lo, nhi, hi  //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                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;}
                ////////// output criteria //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    // standart error criteria of yi : it should be smaller than the tolerance value
                    // calculate the mean value (including biggest
                    pavg=0;
                    for(int i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(iter);
                       for (int j=0; j<=NDIMS; j++)
                          { System.out.print(p[ilo][j]+" ");}
                       System.out.println("");
                    }
                    for(int i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (int j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// calculate ave[] vector excluding the biggest //////

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


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

                double r[] = new double[NDIMS];
                for (int 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 (int j=0; j<NDIMS; j++)
                      p[ihi][j] = r[j];
                    p[ihi][FUNC] = fr;
                    continue;
                }

                if (fr < flo)  //// enlarge, 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 = fnelder.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;
                    }
                }

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

                if (fr < fhi)  
                {
                    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 = fnelder.func(c);
                    if (fc <= fr)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   ///////  shrink:
                    {
                        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] = fnelder.func(p[i]);
                          }
                        continue;
                    }
                }

                if (fr >= fhi)   /// 
                {
                    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 = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// 
                    {
                        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] = fnelder.func(p[i]);
                          }
                    }
                }
            }
       return z;
}
      public static double[] nelder(f_xj fnelder,double a[],double da[],double tolerance)
      {return nelder(fnelder,a,da,300,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[])
      {return nelder(fnelder,a,da,500,1.0e-20,0);}
      
      public static double[] nelder(f_xj fnelder,double x0[])
      {  double [] dx0=new double[x0.length];
         for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
         return nelder(fnelder,x0,dx0,500,1.0e-15,0);
      }
      
      public static double[] nelder(f_xj fnelder,double x0[],double tolerance)
      {
        double [] dx0=new double[x0.length];
        for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
                      return nelder(fnelder,x0,dx0,tolerance);
      }   
      
      public static double[] nelder(f_xj fnelder,double x0[],int maxiteration,double tolerance)
      {
        double [] dx0=new double[x0.length];
        for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
        return nelder(fnelder,x0,dx0,maxiteration,tolerance,0);
      }
      
    public static void main(String args[]) throws IOException
    {
    //String in_name=JOptionPane.showInputDialog(" veri dosyasının isminigiriniz : ");
    String in_name="kalinlik1.txt";
    double a[];
    double x[][],y[]; //curve fitting value
    double xu1[][],yu1[]; ////curve fitting value
    double xu2[][],yu2[];
    double xu3[][],yu3[];
    a=new double[4];
    a[0] = 1.0; 
    a[1] = 1.0;
    a[2] = 1.0;
    a[3] = 1.0;
    fxa f=new fxa(in_name,a);
    int i=0,j=0;
    double p[]=nelder(f,a);
    f=new fxa(in_name,p);
    double tmax=f.xmax;
    double dt=0.1; //time step
    int n=(int)(tmax/dt);
    System.out.println("length="+f.nn+"tmax="+tmax+"n="+n);
    xu1=new double[f.m][n];
    yu1=new double[n];
    xu2=new double[f.m][n];
    yu2=new double[n];
    xu3=new double[f.m][n];
    yu3=new double[n];
    double z[]=new double[2];
    for(i=0;i<n;i++)
    {
    xu1[0][i]=dt*i;
    xu1[1][i]=6;
    z[0]=xu1[0][i];
    z[1]=xu1[1][i];
    yu1[i]=f.Ps(z,p);
    System.out.println("x0="+xu1[0][i]+"x1="+xu1[1][i]+"y="+yu1[i]);
    }  
    for(i=0;i<n;i++)
    {
    xu2[0][i]=dt*i;
    xu2[1][i]=8;
    z[0]=xu2[0][i];
    z[1]=xu2[1][i];
    yu2[i]=f.Ps(z,p);
    System.out.println("x0="+xu2[0][i]+"x1="+xu2[1][i]+"y="+yu2[i]);
    }
    for(i=0;i<n;i++)
    {
    xu3[0][i]=dt*i;
    xu3[1][i]=10;
    z[0]=xu3[0][i];
    z[1]=xu3[1][i];
    yu3[i]=f.Ps(z,p);
    System.out.println("x0="+xu3[0][i]+"x1="+xu3[1][i]+"y="+yu3[i]);
    }    
    Plot pp=new Plot(f.xi[0],f.yi);
    pp.setPlabel("Moistre Ratio - time plot with thickness parameter");
    pp.setXlabel("t, time hour");
    pp.setYlabel("MR Mositure Ratio ");
    pp.setPlotType(0,24); 
    pp.addData(xu1[0],yu1);
    pp.setPlotType(1,0);
    pp.addData(xu2[0],yu2);
    pp.setPlotType(1,1);
    pp.addData(xu3[0],yu3);
    pp.setPlotType(1,2);
    pp.plot();

    String s1=" Optimisation value  : \n"+Matrix.toStringT(p)+"\n";
    s1+="Value of the original function = "+f.func(a)+"\n";
    s1+="Function value = "+f.func(p);
    s1+="standard deviation= "+Math.sqrt(f.func(p)/(f.nn-1));
    String s2="Nelder-Mead non linear least square curve fitting : ";
    System.out.println(s1);
    }
}
