// Broyden - Fletcher-Goldfarb-Shanno (BFGS) metodu 
// Metod Referansı : Intruduction to Optimum Design, Second Edition
// Jasbir S. Arora, Elsevier, ISBN: 0-12-064155-0 
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;

class f1 extends f_xj
{
public double func(double x[]) 
{
//çözümü istenen fonksiyon	
double ff;
ff=3.0*x[0]*x[0]-4.0*x[0]*x[1]+2.0*x[1]*x[1]-x[0]-x[1];
//ff=5.0*x[0]*x[0]+2.0*x[0]*x[1]+x[1]*x[1]+7.0;
return ff; //minimum değil maxsimum istendiğinden - ile çarptık
}
}

class df1 extends fi_xi
{
public double[] func(double x[])
{
//çözümü istenen fonksiyonun türev fonksiyonları		
double ff[]=new double[2];
//ff[0]=10.0*x[0]+2.0*x[1];
//ff[1]=2.0*x[0]+2.0*x[1];
ff[0]=6*x[0]-4.0*x[1]-1;
ff[1]=-4.0*x[0]+4.0*x[1]-1;
return ff;
} 
}



public class OPO12I
{
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[][] V_VT_y_s(double [] left,double [] y,double s[])
{
//multiplys a vector transpose with a vector	
int n=left.length;
double aa[][]=new double[n][n];
double ss=y_s(y,s);
double x1=y_s(y,s);
for(int i=0;i<n;i++)
  {
  for(int j=0;j<n;j++)
    {aa[i][j]=left[i]*left[j]/ss;}
  }
return aa;  

}

public static double y_s(double [] y,double [] s)
{
//scalar multiply two vectors	
int n=y.length;
double aa=0;
for(int i=0;i<n;i++)
  {aa+=y[i]*s[i];} 
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-10;
    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-12;
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;
}

// Broyden - Fletcher-Goldfarb-Shanno (BFGS) metodu 
//
// sadece  fonksiyon verilmiş türev sayısal hesaplanıyor
//________________________________________________________________

public static double[] BFGS( f_xj f,double x[])
{		
int i,j,k;
int nmax=100;
double gamma=1.0;
double tolerance=1.0e-6 ;
int n=x.length;
double alpha;
double s[]=new double[n];
double c[] =new double[n];
double c1[] =new double[n];
double d[] =new double[n];
double y[] =new double[n];
double alpha_d;
double H[][]=I(n);
double H1[][]=I(n);
//double H[][]={{1,0},{0,0}};
double D[][]=new double[n][n];
double E[][]=new double[n][n];
int nn=15;
j=0;k=0; 
c=turev(f,x);
for(i=0;i<c.length;i++)
{d[i]=-c[i];} //first iteration : steepest descent
gamma=norm(c);
while( k++ < nmax && gamma > tolerance )
{   
    alpha=linmin(f,x,d);
    for(i=0;i<n;i++)
    {s[i]=alpha*d[i];
     x[i]+=s[i]; 
    }   
    c1=turev(f,x);    
    for(i=0;i<n;i++)
    {y[i]=c1[i]-c[i];}    
    D=V_VT_y_s(y,y,s);
    E=V_VT_y_s(c,c,d);
    for(i=0;i<n;i++)
    {for(j=0;j<n;j++)
    {H[i][j]=H[i][j]+D[i][j]+E[i][j];
     H1[i][j]=H[i][j];
    }
    }
    for(i=0;i<n;i++)
    {c[i]=c1[i];
     c1[i]=-c1[i];
    }     
    d=pivotlugauss(H1,c1); 
    gamma=norm(c);
}
if(k >= 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;
}

// fonksiyon ve türev fonksiyonu verilmiş
//_______________________________________________________________

public static double[] BFGS( f_xj f,fi_xi df,double x[])
{		
int i,j,k;
int nmax=100;
double gamma=1.0;
double tolerance=1.0e-6 ;
int n=x.length;
double alpha;
double s[]=new double[n];
double c[] =new double[n];
double c1[] =new double[n];
double d[] =new double[n];
double y[] =new double[n];
double alpha_d;
double H[][]=I(n);
double H1[][]=I(n);
//double H[][]={{1,0},{0,0}};
double D[][]=new double[n][n];
double E[][]=new double[n][n];
int nn=15;
j=0;k=0; 
c=df.func(x);
for(i=0;i<c.length;i++)
{d[i]=-c[i];} //first iteration : steepest descent
gamma=norm(c);
while( k++ < nmax && gamma > tolerance )
{   
    alpha=linmin(f,x,d);
    for(i=0;i<n;i++)
    {s[i]=alpha*d[i];
     x[i]+=s[i]; 
    }   
    c1=df.func(x);    
    for(i=0;i<n;i++)
    {y[i]=c1[i]-c[i];}    
    D=V_VT_y_s(y,y,s);
    E=V_VT_y_s(c,c,d);
    for(i=0;i<n;i++)
    {for(j=0;j<n;j++)
    {H[i][j]=H[i][j]+D[i][j]+E[i][j];
     H1[i][j]=H[i][j];
    }
    }
    for(i=0;i<n;i++)
    {c[i]=c1[i];
     c1[i]=-c1[i];
    }     
    d=pivotlugauss(H1,c1); 
    gamma=norm(c);
}
if(k >= 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[] multiply(double[] left,double right)
{
//multiplying a vector with a constant
int i;
int n=left.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
  {
  b[i]=right*left[i];
  }
return b;
}


public static double[] divide(double[] left,double[][] right)
{
return pivotlugauss(right,left);
}

  //en kucuk kareler metodu	
  public static double[] pivotlugauss(double a[][],double b[])
  { //kısmi pivotlu gauss eleme yöntemi
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss eleme
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivotlama
	p=k;
	buyuk=Math.abs(a[k][k]);
	for(ii=k+1;ii<n;ii++)
	{  dummy=Math.abs(a[ii][k]);
	   if(dummy > buyuk) {buyuk=dummy;p=ii;}
	}
	if(p!=k)
	{ for(jj=k;jj<n;jj++)
	  { dummy=a[p][jj];
	    a[p][jj]=a[k][jj];
	    a[k][jj]=dummy;
	  }
	  dummy=b[p];
	  b[p]=b[k];
	  b[k]=dummy;
	}
	//gauss elemeyi çözme  
    for(i=k+1;i<n;i++)
    {  carpan=a[i][k]/a[k][k];
       a[i][k]=0;
       for(j=k+1;j<n;j++)
       {   a[i][j]-=carpan*a[k][j]; }
           b[i]   =b[i]   -carpan*b[k];
	}
  }
  //geriye doğru yerine koyma
  x[n-1]=b[n-1]/a[n-1][n-1];
  for(i=n-2;i>=0;i--)
  {
    toplam=0;
    for(j=i+1;j<n;j++)
    {  toplam+=a[i][j]*x[j];}
    x[i]=(b[i]-toplam)/a[i][i];
  }
  return x;
}

public static double[] add(double[] left,double[] right)
{
//addition of two vectors
int n1=left.length;
int n2=right.length;
int nMax;
int i;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[];
b=new double[nMax];
for(i=0;i<n1;i++)
  {
    b[i]=b[i]+left[i];
  }
for(i=0;i<n2;i++)
  {
    b[i]=b[i]+right[i];
  }
return b;
//end of vector addition method
}

public static double[][] add(double[][] left,double[][] right)
{
//addition of two matrices
int n1=left.length;
int m1=left[0].length;
int n2=right.length;
int m2=right[0].length;
int nMax,mMax;
int i,j;
if(m1>=m2) mMax=m1;
else       mMax=m2;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[][];
b=new double[nMax][mMax];
for(i=0;i<n1;i++)
  {
  for(j=0;j<m1;j++)
    {
    b[i][j]=left[i][j]+right[i][j];
    }
  }
return b;
//end of matrix addition method
}


	 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);
     f1 f_x=new f1();
     //df1 df_x=new df1();
     //double [] r1= BFGS(f_x,df_x,x0);
     double [] r1= BFGS(f_x,x0);
     s=" çözüm seti BFGS: \n"+Matrix.toStringT(r1);
     String s2="Broyden - Fletcher-Goldfarb-Shanno (BFGS) çok değişkenli fonksiyon optimizasyonu:";
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     }      
}
