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

      class yy extends f_xj
      {  
	  double xi[]; // independent variable data set
      double yi[]; //dependent variable data set
      double a[];  //fit function coefficient set
      int nn;
  	  public yy(String filename,double ia[])
      {
        //read the data to curvefit
        //get the data file and initial fit coefficient when class is defined
	    xi=new double[500];
        yi=new double[500];
        int n=ia.length;
        a=new double[n];
        seta(ia);
        int i=-1;
        try{
	    BufferedReader fin=new BufferedReader(new FileReader(filename));
        try {
        while(fin != null)
        {
        i++;
        xi[i]=Text.readDouble(fin);
        yi[i]=Text.readDouble(fin);
        }
        } catch(EOFException e_eof) {System.out.println("end of file"); }
        }catch(IOException e_io) {System.out.println("dosya bulunamadı"); }        
      nn=i;
      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 Cp(double T,double ai[])
      {
      //sample specific function 
	  // saturation pressure as a function of saturation temperature
	  // a = ai[0]
	  // k = ai[1]
	  // b = ai[2]
      a=ai;
      double A1=Math.sinh(ai[2]/T);
      double A2=Math.cosh(ai[4]/T);
      double cp=ai[0]+ai[1]*ai[2]/T/(A1*A1)+a[3]*a[4]/T/(A2*A2);
      return cp; //kPa
      }      
      
      public double func(double ai[])
      {
      double ff=0;
      double w;
      double yy;
       for(int i=0;i<nn;i++)
	    {w= Cp(xi[i],ai)-yi[i];
        ff+=w*w;
	    }
      return ff;
	  }
}

public class SCO11H
{//2. Nelder ve Mead metodu (simpleks arama metodu)
//________________________________________________________________
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 çok boyutlu simplex minimizasyon metodu
   //  Nelder & Mead 1965 Computer J, v.7, 308-313.

   // Giriş değişkenleri tanımlaması
   // fnelder : abstract çok boyutlu fonksiyon f(x)
   // x : for n boyutlu n+1 simlex noktasını içerir bağımsız değişken seti
   // maxiteration : maximum iterasyon sayısı
   // tolerance : 
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// başlangıç simplexini oluştur //////////////////
            double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// iilk fonksiyon değerlerini hesapla ////////////////
            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++)
            {
                ///////////  lo, nhi, hi noktalarını tanımla //////////////
                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;}
                ////////// çıkış kriteri //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    //yi nin standart hata kriteri set değerinden (tolerance)
                    // küçük olmalı
                    // ortalama değeri hesapla (en büyük değer de dahil olmak üzere)
                    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];
                }

                ///// ave[] vektorünü en büyük değeri hariç tutarak hesapla //////

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


                ///////// yansıt ////////////////

                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)  //// genişlet, else kabul et
                {
                    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;
                    }
                }

                ///////////// daralt:

                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   /////// daralt
                    {
                        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 void main(String args[]) throws IOException
    {
    String in_name=JOptionPane.showInputDialog(" veri dosyasının ismini giriniz : ");
    //String in_name="D1.txt";
    double x[],y[]; //uydurulan veri
    double xu[],yu[]; //uydurulan veri
    double a[];
    a=new double[5];
    a[0] = 1.0; 
    a[1] = 1.0;
    a[2] = 1.0; 
    a[3] = 1.0; 
    a[4] = 1.0;   
    double da[];
    da=new double[5];
    da[0] =0.3;
    da[1] =0.3;
    da[2] =0.3;
    da[3] =0.3;
    da[4] =0.3;
    yy f=new yy(in_name,a);
    int i=0;
    double p[]=nelder(f,a,da,1000,1e-15,1);
    String s1=" optimizasyon değeri  : \n"+Matrix.toStringT(p)+"\n";
    s1+="orijinal fonksiyon değeri = "+f.func(a)+"\n";
    s1+="fonksiyon değeri = "+f.func(p)+"\n";
    s1+="standart sapma= "+Math.sqrt(f.func(p)/(f.nn-1));
    f=new yy(in_name,p);
    i=0;
    double tmax=f.xi[f.nn-1];
    double dt=0.01; //ZAMAN ARALIĞI
    int n=(int)(tmax/dt);
    System.out.println("length="+f.nn+"tmax="+tmax+"n="+n);
    xu=new double[n];
    yu=new double[n];
    x=new double[f.nn];
    y=new double[f.nn];
    for(i=0;i<f.nn;i++)
    {x[i]=f.xi[i];y[i]=f.yi[i];}
    double t=0;
    for(i=0;i<n;i++)
    {t=dt*i;//System.out.println("t="+t+"ti="+f.xi[i]+"Pi="+f.yi[i]+"P="+f.Ps(t,p))
    xu[i]=t;
    yu[i]=f.Ps(t,p);
    System.out.println("x="+xu[i]+"y="+yu[i]);
    };  
    Plot pp=new Plot(xu,yu);
    pp.setPlabel("lineer olmıyan 2 boyutlu Nelder-mead optimizasyon metodu ile eğri uydurma");
    pp.setXlabel("x");
    pp.setYlabel("y");
    pp.setPlotType(0,0); 
    pp.addData(x,y);
    pp.setPlotType(1,24);// üçgen olarak çiz
    pp.plot();
    String s2="Nelder-Mead en küçük kareler lineer olmiyan eğri uydurma : ";
    //JOptionPane.showMessageDialog(null,s1,s2,JOptionPane.PLAIN_MESSAGE);
    System.out.println(s1);
    }

}
