// Newton metodu ile optimizasyon
// 
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=2.0*x[0]*x[1]+2.0*x[0]-x[0]*x[0]-2.0*x[1]*x[1];
return -ff; //minimum değil maxsimum istendiğinden - ile çarptık
}
}
class f2 extends fi_xi
{
public double[] func(double x[]) 
{
//çözümü istenen fonksiyon	
double ff[]={-(2.0*x[1]+2.0-2.0*x[0]),-(2.0*x[0]-4.0*x[1])};
return ff;
}
}

class f3 extends fij_xi
{
public double[][] func(double x[]) 
{
//çözümü istenen fonksiyon	
double ff[][]={{2,-2},{-2,-4}};
return ff; //minimum değil maxsimum istendiğinden - ile çarptık
}
}

public class OPO12D
{
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[] 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]=dfdx(f,xxi,i);
   }
return a;   	   
}


public static double dfdx(f_xj f,double xxi[],int xref)
{ int i;
	double h=1.0e-3;
  double hh=1.0/h; 
  double xp2h[]=new double[xxi.length];
  double xph[]=new double[xxi.length];
  double xm2h[]=new double[xxi.length];
  double xmh[]=new double[xxi.length];
  double fp2h,fph,fm2h,fmh;
  for(i=0;i<xxi.length;i++)
  {
  if(i!=xref)
    {xp2h[i]=xxi[i];xph[i]=xxi[i];xm2h[i]=xxi[i];xmh[i]=xxi[i];}
  else {xp2h[i]=xxi[i]+2.0*h;xph[i]=xxi[i]+h;xm2h[i]=xxi[i]-2.0*h;xmh[i]=xxi[i]-h;}  
  }
  fp2h=f.func(xp2h);
  fph=f.func(xph);
  fm2h=f.func(xm2h);
  fmh=f.func(xmh);
  return (-fp2h+8.0*fph-8.0*fmh+fm2h)/12.0*hh;
}
public static double d2fdx2(f_xj f,double xxi[],int xiref,int xjref)
{ //ikincitürev matrisi
  int i;
  double h=1.0e-4;
  double hh=1.0/h;
  int n=xxi.length; 
  double ff;
  double xp2h[]=new double[xxi.length];
  double xph[]=new double[xxi.length];
  double xm2h[]=new double[xxi.length];
  double xmh[]=new double[xxi.length];
  double fp2h,fph,fm2h,fmh;
     for(int j=0;j<n;j++)
     {if(j!=xjref) {xp2h[j]=xxi[j];xph[j]=xxi[j];xm2h[j]=xxi[j];xmh[j]=xxi[j];}
      else {xp2h[j]=xxi[j]+2.0*h;xph[j]=xxi[j]+h;xm2h[j]=xxi[j]-2.0*h;xmh[j]=xxi[j]-h;}  
     }  
     fp2h=dfdx(f,xp2h,xiref);
     fph=dfdx(f,xph,xiref);
     fm2h=dfdx(f,xm2h,xiref);
     fmh=dfdx(f,xmh,xiref);
     ff=(-fp2h+8.0*fph-8.0*fmh+fm2h)/12.0*hh;
  return ff;
}

public static double[][]H(f_xj f,double x[])
{
  int n=x.length;	
  double dy[][]=new double[n][n];;
  for(int ii=0;ii<n;ii++)
  {
    for(int jj=0;jj<n;jj++)
    {   //Hessian matrisi
        dy[ii][jj]=d2fdx2(f,x,ii,jj);
    }
  }
  return dy;
}
// 4. modifiye Newton metodu (en dik yamaç ve newton metodlarını birleştirir)
//________________________________________________________________

public static double[] newton( f_xj f,fi_xi df,fij_xi d2f,double x[])
{
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-10;
int n=x.length;
double b[];
b=new double[n];
i=0;
for(i=0;i<n;i++)
{b[i]=1.0;}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    b=Matrix.multiply(Matrix.divide(df.func(x),d2f.func(x)),-ti);
    x=Matrix.add(x,b);
}
if(i >= 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[] dampednewton(f_xj f,fi_xi df,fij_xi d2f,double x[])
{
//  Damped Newton metodu
//  Newton-Raphson metodu ile en dik yamaç metodunu birlikte kullanır.
// linmin ve brent metodunukullanır	
double ti=1.0;
int i,ii,jj;
int nmax=800;
double tolerance=1.0e-14;
int n=x.length;
double b[];
b=new double[n];
double lambda;
double g[]=new double[n];
double G[][]=d2f.func(x);//Hessian matrix
double g1,g2,g3,g4,g5;
//dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    g=df.func(x);
    
    g1=Matrix.VT_V(g);
    g2=Matrix.VT_G_V(g,G);
    g1*=g1;
    g3=g1/g2;
    g4=Matrix.VT_G_V(g,Matrix.inverse(G));
    lambda=linmin(f,x,g);
    if(g4>Math.max(g3,0))
          b=Matrix.multiply(Matrix.divide(g,G),lambda);
    else
          b=Matrix.multiply(g,lambda);
    x=Matrix.add(x,b);
}
if(i >= 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[] 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[] 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 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 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;	
}
// veri giriş çıkış metodu
//________________________________________________________________
     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=new f1();
     f2 df=new f2();
     f3 d2f=new f3();
     double r1[]=dampednewton(f,df,d2f,x0);
     s=" çözüm seti  : \n"+Matrix.toStringT(r1);
     String s2="Newton metodu ile lineer olmıyan çok değişkenli fonksiyon optimizasyonu:";
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}
