// OPO12F Davidon-Fletcher-Powell optimizasyonu
// Metod Referansı : Numerical Analysis, 8th edition Richard L. Burden, J. Douglas Faires 3
// sayfa 617
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
import java.io.*;

      class fa extends f_xj
      {  
	  double xi[]; // independent variable data set
      double yi[]; //dependent variable data set
      double a[];  //fit function coefficient set
      int nn;
  	  public fa(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 Ps(double x,double ai[])
      {
      //sample specific function 
	  // saturation pressure as a function of saturation temperature
      a=ai;
      double pi=ai[0]*(Math.exp(-ai[1]*x))+ai[2]*(Math.exp(ai[3]*x*x));
      return pi; //kPa
      }      
      
      public double func(double ai[])
      {
      double ff=0;
      double w;
      double yy;
       for(int i=0;i<nn;i++)
	    {w=Ps(xi[i],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 OPO12F2
{
public static double VT_V(double [] left)
{
//multiplys a vector transpose with a vector	
int n=left.length;
double tot=0;
for(int i=0;i<n;i++)
  {
  tot+=left[i]*left[i];
  }
return tot;  
}

public static double[][] V_VT(double [] left)
{
//multiplys a vector transpose with a vector	
int n=left.length;
double aa[][]=new double[n][n];
for(int i=0;i<n;i++)
  {
  for(int j=0;j<n;j++)
    {
	aa[i][j]=left[i]*left[j];
	}
  }
return aa;  
}

public static double  VT_X(double [] left,double [] right)
{
//multiplys a vector transpose with a vector	
int n=left.length;
double tot=0;
for(int i=0;i<n;i++)
  {
  tot+=left[i]*right[i];
  }
return tot;  
}

public static double norm(double v[])
{
// vector norm
double total=0;;
for(int i=0;i<v.length;i++)
  {total+=v[i]*v[i];}
return Math.sqrt(total);
}

public static double turev(f_xj f,double x[],int denklem_ref,int x_ref)
{
// 
// verilen fonksiyonun ikinci türevlerinin matrisi
// fonksiyon f in denklem_ref sayılı fonksiyonunun(sayılar 0 dan başlar)
// x[x_ref] değişkenine göre türevi (sayılar 0 dan başlar) 
// df_denklem_ref(x)/d_x_ref
// bu metod newtond metodu içinde kullanılmak içindir.
double h0=0.256808;
int i,m;
int n=7;
double f1[];
f1=new double[x.length];
double f2[];
f2=new double[x.length];
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//turev of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i]*1.0094847;
f1=df(f,x1);
f2=df(f,x2);
T[i][0]=( f1[denklem_ref] - f2[denklem_ref])/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}
public static double[] df(f_xj f,double xxi[])
{
//Gredyen vektörü
// f(x0,x1,x2,...xn) founksiyonunun türev vektörü
// [df/dx0, df/dx1,df/dx2,....,df/dxn]
// df/dxj  j=0...x.length
//
// bu fonksiyon bilgisayar tarafından hesaplanmaktadır
// kullanıcı tarafından tanımlanan f(x) fonksiyonunun
// sayısal olarak alınan türevidir.	

// 
double a[]=new double[xxi.length];
for(int i=0;i<xxi.length;i++)
   {
   a[i]=turev(f,xxi,i);
   }
return a;   	   
}
public static double turev(f_xj f,double xxi[],int x_ref)
{
// df/dxj
double h0=0.256808;
int i,m;
int n=7;
double x1[];
x1=new double[xxi.length];
double x2[];
x2=new double[xxi.length];
for(i=0;i<xxi.length;i++)
{
x1[i]=xxi[i];
x2[i]=xxi[i];
}
//turev of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
//h[i]=h0*Math.pow(r,i);
h0/=2.0;
h[i]=h0;
}
//first turev (difference formula)
for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
T[i][0]=(f.func(x1)-f.func(x2))/(2.0*h[i]);
x1[x_ref]=xxi[x_ref];
x2[x_ref]=xxi[x_ref];
}

for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] -
  h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i] - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}
public static double[] turev(f_xj f_deriv,double x[])
{
// df/dxj  j=0...x.length
// This method calculates turev of a function with more than one independent variable. 
// Accuracy of method can be adjusted by changing variables h0 and n
// function input should be in the form given in abstract class
// f_xj,j=0...x.length = df/dx(x_ref)
double a[]=new double[x.length];
for(int x_ref=0;x_ref<x.length;x_ref++)
   {
   a[x_ref]=turev(f_deriv,x,x_ref);
   }
return a;   	   
}

public static double SIGN(double a,double b) {return (b > 0.0 ? Math.abs(a) : -Math.abs(a));}
public static double MAX(double a,double b) {return (a > b ? a : b);}

	 
public static double linmin(f_xj f,double p[],double xi[])
{
	double tol=2.0e-4;
    int n=p.length;	
    double pcom[],xicom[];
	int j;
	double xx,xmin,bx,ax;
	f1dim f1=new f1dim(f,p,xi);
	ax=0.0;
	xx=1.0;
	bx=2.0;
	double aa[]=mnbrak(f1,ax,xx,bx);
	ax=aa[0];
	xx=aa[1];
	bx=aa[2];
	xmin=brent(f1,ax,xx,bx,tol);
	return xmin;
}

public static double[][] linminiter(f_xj f,double p[],double xi[])
{   // yeni iterasyon noktasını verir
    int n=p.length;
    double xmin=linmin(f,p,xi);
    double aa[][]=new double[2][n];
	for (int j=0;j<n;j++) 
	{   xi[j] *= xmin;
		p[j] += xi[j];
		aa[0][j]=p[j];
		aa[1][j]=xi[j];
	}
    return aa;	
}
public  static double[] brentf(f_x f,double ax,double bx,double cx,double tol)
{
// ax,bx,cx üç adet tahmin noktası (bxax ile cx arasında yer almalıdır)
// f : fonksiyon (f_x abstract sınıfında tanımlanmış)
// tol tolerans	
int ITMAX=100;
double CGOLD=(3.0-Math.sqrt(5))/2.0; //altın oran
double ZEPS=1.0e-10;
double xmin;
double aa[]=new double[2];
// SHFT(double a,double b,double c,d) (a)=(b);(b)=(c);(c)=(d);
	int iter;
	double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
	double e=0.0;
	d=0.0;
	a=((ax < cx) ? ax : cx);
	b=((ax > cx) ? ax : cx);
	x=w=v=bx;
	fw=fv=fx=f.func(x);
	for (iter=1;iter<=ITMAX;iter++) {
		xm=0.5*(a+b);
		tol2=2.0*(tol1=tol*Math.abs(x)+ZEPS);
		if (Math.abs(x-xm) <= (tol2-0.5*(b-a))) {
			xmin=x;
	        aa[0]=xmin;
	        aa[1]=fx;
	        return aa;
		}
		if (Math.abs(e) > tol1) {
			r=(x-w)*(fx-fv);
			q=(x-v)*(fx-fw);
			p=(x-v)*q-(x-w)*r;
			q=2.0*(q-r);
			if (q > 0.0) p = -p;
			q=Math.abs(q);
			etemp=e;
			e=d;
			if (Math.abs(p) >= Math.abs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
				d=CGOLD*(e=(x >= xm ? a-x : b-x));
			else {
				d=p/q;
				u=x+d;
				if (u-a < tol2 || b-u < tol2)
					d=SIGN(tol1,xm-x);
			}
		} else {
			d=CGOLD*(e=(x >= xm ? a-x : b-x));
		}
		u=(Math.abs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
		fu=f.func(u);
		if (fu <= fx) {
			if (u >= x) a=x; else b=x;
			{v=w;w=x;x=u;}
			{fv=fw;fw=fx;fx=fu;}
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				v=w;
				w=u;
				fv=fw;
				fw=fu;
			} else if (fu <= fv || v == x || v == w) {
				v=u;
				fv=fu;
			}
		}
	}
	System.out.println("BRENT metodunda maksimum iterasyon sayısı aşıldı");
	xmin=x;
	//minumum değer a[0] ile fonksiyon değeri a[1] ile geridöndürülmektedir
	aa[0]=xmin;
	aa[1]=fx;
	return aa;
}
public  static double[] mnbrak(f_x f,double ax,double dx)
{
// dx ikinci derece polinom için gerekli iki ek noktayı oluşturmak için 
// ana noktadan seçilen uzaklık	
return mnbrak(f,ax,(ax-dx),(ax+dx));
}

public  static double[] mnbrak(f_x f,double ax)
{
// dx ikinci derece polinom için gerekli iki ek noktayı oluşturmak için 
// ana noktadan seçiln uzaklık	
return mnbrak(f,ax,(ax-0.1),(ax+0.1));
}

public  static double[] mnbrak(f_x f,double ax,double bx,double cx)
{
// ax,bx,cx üç adet tahmin noktası verildiğinde içinde minimum bulunan ax,bx,cx üç adetnoktayı bize arama
// sonucu bulur 
// f : fonksiyon (f_x abstract sınıfında tanımlanmış)
// tol tolerans	
int ITMAX=100;
double xmin;
double GOLD=(Math.sqrt(5.0)+1.0)/2.0;;
double GLIMIT=100.0;
double TINY=1.0e-20;
double fa,fb,fc;
fa=0;
fb=0;
fc=0;
	double ulim,u,r,q,fu,dum;
	double aa[]=new double[3];
	fa=f.func(ax);
	fb=f.func(bx);
	if (fb > fa) {
		//SHFT(dum,*ax,*bx,dum)
		{dum=ax;ax=bx;bx=dum;}
		//SHFT(dum,*fb,*fa,dum)
		{dum=fb;fb=fa;fa=dum;}
	} 
	cx=(bx)+GOLD*(bx-ax);
	fc=f.func(cx);
		aa[0]=ax;
		aa[1]=bx;
		aa[2]=cx;
	while (fb > fc) {
		r=(bx-ax)*(fb-fc);
		q=(bx-cx)*(fb-fa);
		u=(bx)-((bx-cx)*q-(bx-ax)*r)/
			(2.0*SIGN(MAX(Math.abs(q-r),TINY),q-r));
		ulim=(bx)+GLIMIT*(cx-bx);
		if ((bx-u)*(u-cx) > 0.0) {
			fu=f.func(u);
			if (fu < fc) {
				ax=bx;
				bx=u;
				fa=fb;
				fb=fu;
				return aa;
			} else if (fu > fb) {
				cx=u;
				fc=fu;
				return aa;
			}
			u=(cx)+GOLD*(cx-bx);
			fu=f.func(u);
		} else if ((cx-u)*(u-ulim) > 0.0) {
			fu=f.func(u);
			if (fu < fc) {
				//SHFT(bx,cx,u,cx+GOLD*(cx-bx))
				{bx=cx;cx=u;u=cx+GOLD*(cx-bx);}
				//SHFT(fb,fc,fu,f.func(u))
				{fb=fc;fc=fu;fu=f.func(u);}
			}
		} else if ((u-ulim)*(ulim-cx) >= 0.0) {
			u=ulim;
			fu=f.func(u);
		} else {
			u=cx+GOLD*(cx-bx);
			fu=f.func(u);
		}
		//SHFT(ax,bx,cx,u)
		{ax=bx;bx=cx;cx=u;}
		//SHFT(fa,fb,fc,fu)
		{fa=fb;fb=fc;fc=fu;}
		aa[0]=ax;
		aa[1]=bx;
		aa[2]=cx;
	}
		return aa;
}	
public  static double brent(f_x f,double ax,double bx,double cx,double tol)
{
 double aa[]=brentf(f,ax,bx,cx,tol);
 return aa[0];
}

public static double[][] I(int n)
{
//unit matrix
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=1.0;
return b;
}
// Davidon - Fletcher-Powell metodu 
//________________________________________________________________
public static double[] davidon_fletcher_powell( f_xj f,double x[])
{		
int i,j,k;
int nmax=500;
double gamma=1.0;
double tolerance=1.0e-12 ;
int n=x.length;
double alpha;
double dx[]=new double[n];
double g[] =new double[n];
double ge[]=new double[n];
double dg[]=new double[n];
double d[] =new double[n];
double r[] =new double[n];
double rdg,dxdg;
double beta;
double alpha_d;
double Q[][]=I(n);
double Q1[][]=new double[n][n];
double Q2[][]=new double[n][n];
int nn=15;
i=0;
j=0;k=0; 
g=turev(f,x);
while( j++ < nmax && gamma > tolerance )
{
    for(int ii=0;ii<n;ii++)
    {
     d[ii]=0.0;
     for(int kk=0;kk<n;kk++)
	    {d[ii]=d[ii]-Q[ii][kk]*g[kk];}
    }
    System.out.println("j="+j+"d="+Matrix.toString(d));
    alpha=linmin(f,x,d);  
    System.out.println("alfa="+alpha);
    if(alpha==0) break;
    for(i=0;i<n;i++){dx[i]=alpha*d[i];x[i]+=dx[i];ge[i]=g[i];}
    System.out.println("j="+j+"dx=\n"+Matrix.toString(dx));
    System.out.println("j="+j+"x=\n"+Matrix.toString(x));  
    g=turev(f,x);
    System.out.println("j="+j+"g=\n"+Matrix.toString(g));
    for(i=0;i<n;i++){dg[i]=g[i]-ge[i];}
    gamma=norm(g);
    for(int ii=0;ii<n;ii++)
    {r[ii]=0.0;
	    for(int kk=0;kk<n;kk++)
	    {r[ii]+=Q[ii][kk]*dg[kk];}
    }
    rdg=VT_X(r,dg);
    dxdg=VT_X(dx,dg); 
    Q1=V_VT(r);
    Q2=V_VT(dx);
    for(int ii=0;ii<n;ii++)
      for(int jj=0;jj<n;jj++)
        Q[ii][jj]=Q[ii][jj]-Q1[ii][jj]/rdg+Q2[ii][jj]/dxdg;
    gamma=norm(dx);
    System.out.println("j="+j+"Q=\n"+Matrix.toString(Q));
j++;k++;
if(k>=nn) {k=0;Q=I(n);}
}
if(j >= nmax) JOptionPane.showMessageDialog(null,"Uyarı maximum iterasyon sayısı aşıldı \n"+
  " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE);
return x;
}
public static double[] davidon_fletcher_powell_V1( f_xj f,double x[])
{		
int i,j,k;
int nmax=500;
double gamma=1.0;
double tolerance=1.0e-12 ;
int n=x.length;
double alpha;
double dx[]=new double[n];
double g[] =new double[n];
double ge[]=new double[n];
double dg[]=new double[n];
double d[] =new double[n];
double r[] =new double[n];
double rdg,dxdg;
double beta;
double alpha_d;
double Q[][]=I(n);
double Q1[][]=new double[n][n];
double Q2[][]=new double[n][n];
int nn=15;
i=0;
j=0;k=0; 
g=turev(f,x);
while( j++ < nmax && gamma > tolerance )
{
    
    for(int ii=0;ii<n;ii++)
    {
     d[ii]=0.0;
     for(int kk=0;kk<n;kk++)
	    {d[ii]=d[ii]-Q[ii][kk]*g[kk];}
    }
    alpha=linmin(f,x,d);    
    if(alpha==0) break;
    for(i=0;i<n;i++){dx[i]=alpha*d[i];x[i]+=dx[i];ge[i]=g[i];}
    g=turev(f,x);
    for(i=0;i<n;i++){dg[i]=g[i]-ge[i];}
    gamma=norm(g);
    for(int ii=0;ii<n;ii++)
    {  r[ii]=0.0;
	    for(int kk=0;kk<n;kk++)
	    {r[ii]+=Q[ii][kk]*dg[kk];}
    }
    rdg=VT_X(r,dg);
    dxdg=VT_X(dx,dg); 
    beta=dxdg/rdg;
    Q1=V_VT(r);
    Q2=V_VT(dx);
    for(int ii=0;ii<n;ii++)
      for(int jj=0;jj<n;jj++)
        Q[ii][jj]=beta*(Q[ii][jj]-Q1[ii][jj]/rdg)+Q2[ii][jj]/dxdg;
    gamma=norm(dx);
j++;k++;
if(k>=nn) {k=0;Q=I(n);}
}
if(j >= nmax) JOptionPane.showMessageDialog(null,"Uyarı maximum iterasyon sayısı aşıldı \n"+
  " çözüm geçerli olmıyabilir","MAKSİMUM ITERASYON SAYISI UYARISI",JOptionPane.WARNING_MESSAGE);
return x;
}
	 public static double[] verigir(String s)
     {
     String s1=JOptionPane.showInputDialog(s);        
     StringTokenizer token=new StringTokenizer(s1);
     int n=token.countTokens()-1;
     int m=n+1;
     double a[]=new double[m];
     int j=0;           
     while(token.hasMoreTokens())
     {
     Double ax=new Double(token.nextToken());
     a[j++]=ax.doubleValue();
     }
     return a; 
     }
	 public static void main (String args[])  
     {
	 String s="kök tahmin değerlerini arada bir boşluk bırakarak giriniz: ";
     double [] x0=verigir(s);
     fa f_x=new fa("a.txt",x0);
     double [] r1= davidon_fletcher_powell(f_x,x0);
     double [] r2= davidon_fletcher_powell_V1(f_x,x0);
     s=" çözüm seti DFP : \n"+Matrix.toStringT(r1);
     s+=" çözüm seti düzeltilmiş DFP : \n"+Matrix.toStringT(r2);
     String s2="Davidon - Fletcher - Powell çok değişkenli fonksiyon optimizasyonu:";
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     }      
}
