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

import java.io.*;
import java.util.*;
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
  	  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);
        double b[][]=Text.readDoubleT(filename);
        
        //System.out.println("m="+m);
        xi=b[0];
        yi=b[1];
        nn=yi.length;
        //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;
	  }	  

      public  double Ps(double t)
      {//ilk tanım Ps formülü	
      double T=t+273.15;
      double ln10=Math.log(10.0);
      double A[]={40.69889*ln10,-2.36254e3*ln10,-13.06883,7.616005e-3*ln10,0.2342564,3.761111e2};
      double p1=A[0]+A[1]/T+A[2]*Math.log(T)+A[3]*T+A[4]*((A[5]-T)/T)*Math.log(A[5]-T);
      return Math.exp(p1);
      }

      public double dPsdT(double t)
      {
      double T=t+273.15;
      double ln10=Math.log(10.0);
      double A[]={40.69889*ln10,-2.36254e3*ln10,-13.06883,7.616005e-3*ln10,0.2342564,3.761111e2};
      double p1=A[0]+A[1]/T+A[2]*Math.log(T)+A[3]*T+A[4]*((A[5]-T)/T)*Math.log(A[5]-T);
      p1=Math.exp(p1);
      double p2=-A[1]/(T*T)+A[2]/T+A[3]-A[4]/T-A[4]*(-A[5]/(T*T)-1)*Math.log(A[5]-T);
      p2=p2*p1;
      return p2;	
      }

      public double ortalama(double x[])
      {//ortalama değerleri hesaplar
      int N=x.length;
      double tot=0;
      for(int i=0;i<N;i++) tot+=x[i];
      return tot/N;  
      }

      public double s(double x[],int j,int m)
      {// j çevresinde +/- m/2 veri bölgesinde 
       // yerel standart sapma değerlerini hesaplar
       int N=x.length;
       int m1,m2,m3,m4;
       m1=m/2;
       m2=j-m1;
       m3=j+m1;
       if(m2<0) {m2=0;m3=m;if(m3>N){m3=N;}}
       else if(m3>N) {m3=N;m2=m3-m;if(m2<0){m2=0;}} 
       double tot=0;
       m4=m3-m2;
       double y[]=new double[m];
       for(int i=m2;i<m3;i++) {y[i-m2]=x[i];}
       double yort=ortalama(y);
       double yi=0;
       for(int i=0;i<m;i++)
       {yi=y[i]-yort;
        tot+=yi*yi;
       }
       return Math.sqrt(tot/(m-1));  
       }

       public double W(int i,int m)
       {   //weight function
	   double dd=s(xi,i,m)+dPsdT(xi[i])*s(yi,i,m);
	   return 1/(dd*dd);
       }

	  double Ps(double t,double ai[])
      {
	  // t derece C 
      // sample specific function m dimensional independent variable
      // independent variables : x[0],x[1]...x[m-1]
      double T=t+273.15;
      double b0=ai[0];
      double b1=ai[1];
      double a4=ai[2];
      double a5=ai[3];
      double a6=ai[4];
      double a7=ai[5];
      double Tc=374.23;
      double Tt=169.861;
      double Pt=0.42; //kPa
      double Pc=4056;
      double A=(Tc-Tt)/Tt;
      double tt=(T-Tt)/(Tc-Tt);
      double R=8.314471;
      double s=b0/R;
      double r=b1/(R*Tt)-s;
      double q=(1+A*tt);
      double rat=A*r*tt/q;
      //System.out.println("R="+R+"s="+s+"r="+r+"q="+q+"rat="+rat);
      double P0=1.0-Pt/(Pc-Pt)*(1-Math.pow(q,s)*Math.exp(rat));
      double tf=1-tt;
      double teta=0.199;
      double P1=2-a4*tf+a5*Math.pow(tf,(2.0-teta))+ a6*tf*tf*tf+a7*tf*tf*tf*tf;
      double N=87.0*Tt/Tc;
      double pi=Math.pow((Math.pow(P0,N)+Math.pow(P1,N)),(1.0/N));
      //System.out.println("P0="+P0+"P1="+P1+"pi="+pi+"N="+N);
      return (pi-1.0)*(Pc-Pt)+Pt; //kPa
      }
	        
      public double func(double ai[])
      {
      double ff=0;
      double dP;
      double yy;
      double xx[]=new double[xi.length];
      double Wi;
       for(int i=0;i<nn;i++)
	    {   
	    dP=Ps(xi[i],ai);
	    yy=yi[i];
        //System.out.println("t="+xi[i]+"dP="+dP+"yy="+yy+"Pst="+Ps(xi[i]));
	    dP-=yy; 
	    Wi=W(i,30);
        ff+=Wi*dP*dP;
        //System.out.println("x="+xi[i]+"yy="+yy+"dP="+dP+"ff="+ff+"Wi="+Wi); 
	    }
      return ff;
	  }
}



public class SCO11H1
{
   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,5000,1.0e-20,0);}

	
	public static void main(String args[]) throws IOException
    {
    //String in_name=JOptionPane.showInputDialog(" veri dosyasının isminigiriniz : ");
    String in_name="R134a_Ps.txt";
    double a[];
    a=new double[6];
    a[0] = 50;//-45.37032;
    a[1] = 26233.885;
    a[2] = 4.164859;
    double a4=a[2];
    a[3] = -0.11599104+0.29506258*a4*a4-0.00021222*a4*a4*a4*a4*a4;
    a[4] = -0.01546028+0.0897816*a4*a4-0.05322199*a4*a4*a4;
    a[5] = -0.05725757-0.06817687*a4+0.00047188*a4*a4*a4*a4*a4;  
    double da[];
    da=new double[6];
    da[0] =100;
    da[1] =100000;
    da[2] =100;
    da[3] =100;
    da[4] =100;
    da[5] =100;
    fxa f=new fxa(in_name,a);
    int i=0;
    double t=100;
    System.out.println("t="+t+"p="+f.Ps(t,a)+"F="+f.func(a));
    
    double p[]=nelder(f,a,da,2000,1e-30,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);
    String s2="Nelder-Mead en küçük kareler lineer olmiyan eğri uydurma : ";
    JOptionPane.showMessageDialog(null,s1,s2,JOptionPane.PLAIN_MESSAGE);
    }

}
