// FSCO5H Newton Raphson metoduyla
// birden fazla değişkenli denklem sisteminin köklerini bulmak
// ve optimizasyon Hessian matris H sayısal olarak hesaplanmış

import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
//============= Tanımı gereken fonksiyon ================

class fxa 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];
ff=3.0*x[0]*x[0]-4.0*x[0]*x[1]+2.0*x[1]*x[1]-x[0]-x[1];
return ff; //minimum değil maxsimum istendiğinden - ile çarptık
}
}


class fxb extends fi_xi
{
f_xj fxx;	

fxb(f_xj fxi) 
{fxx=fxi;}
 	
public double[] func(double x[])
{
// Türev hesabı
double h=0.0013421432165;
double hh=1.0/h;
int i,m;
int n=x.length;
double f1;
double f2;
double f3;
double f4;
double x1[];
x1=new double[n];
double x2[];
x2=new double[n];
double x3[];
x3=new double[n];
double x4[];
x4=new double[n];
double fx[]=new double[n];

for(int j=0;j<n;j++)
{
for(i=0;i<n;i++)
{
x1[i]=x[i];
x2[i]=x[i];
x3[i]=x[i];
x4[i]=x[i];
}
x1[j]+=2.0*h;
x2[j]+=h;
x3[j]-=h;
x4[j]-=2.0*h;

f1=fxx.func(x1);
f2=fxx.func(x2);
f3=fxx.func(x3);
f4=fxx.func(x4);
fx[j]=(-f1+8.0*f2-8.0*f3+f4)/12.0*hh;
}
return fx;
} 

}

public class OPO12E3
{
// Newton Raphson Metodunu kullanarak denklem sisteminin köklerini bulmak	
// Bu metod denklem sistemi çözümü için Matrix.java programını
// kullanmaktadır.

public static double[] nelder_newton(f_xj fnelder,double a[])
{
	//Nelder-Meads ve Fletcher- Reeves metodlarının kombine kullanımı
double x0[]=nelder(fnelder,a,0.1);
fxb f2 =new fxb(fnelder);
return newtond(f2,x0);
}

public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance,int printlist)
{
	int i,j;
double x[][]=new double[a.length+1][a.length];
double p[][]=new double[a.length+1][a.length+1];
for(i=0;i<x.length;i++)
   {for(j=0;j<x[0].length;j++)
     {if(i==j){x[i][j]=a[j]+da[j];p[i][j]=x[i][j];}
      else    {x[i][j]=a[j];p[i][j]=x[i][j]; }
     }
     p[i][j] = fnelder.func(p[i]);
   }
   //System.out.println("x=\n"+Matrix.toString(x));
   // 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]);
              System.out.println("i="+i+"p=\n"+Matrix.toString(p));
            }
            */
            int iter=0;
            
            for (iter=1; iter<maxiteration; iter++)
            {
                System.out.println("iter="+iter+"p=\n"+Matrix.toString(p));
	            ///////////  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 (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 (i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                 System.out.println("büyük y i="+fhi+" "+ihi+" ikinci büyük y i="+fnhi+" "+inhi+" küçük y i="+flo+" "+ilo);    
                ////////// çı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(i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(iter);
                       for (j=0; j<=NDIMS; j++)
                          { System.out.print(p[ilo][j]+" ");}
                       System.out.println("");
                    }
                    for(i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    if(sterr < tolerance)
                      { for (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 (j=0; j<NDIMS; j++)
                  ave[j] = 0;
                for (i=0; i<NPTS; i++)
                  if (i != ihi)
                    for (j=0; j<NDIMS; j++)
                       ave[j] += p[i][j];
                for (j=0; j<NDIMS; j++)
                   ave[j] /= (NPTS-1);


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

                double r[] = new double[NDIMS];
                for (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 (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 (j=0; j<NDIMS; j++)
                      e[j] = 3*ave[j] - 2*p[ihi][j];
                    double fe = fnelder.func(e);
                    if (fe < fr)
                    {
                       for (j=0; j<NDIMS; j++)
                         p[ihi][j] = e[j];
                       p[ihi][FUNC] = fe;
                       continue;
                    }
                    else
                    {
                       for (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 (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 (j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   /////// daralt
                    {
                        for (i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (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 (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 (j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// 
                    {
                        for (i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (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,500,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[])
      {return nelder(fnelder,a,da,500,1.0e-15,0);}
      
      public static double[] nelder(f_xj fnelder,double x0[])
      { 
        double [] dx0=new double[x0.length];
        for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
	      return nelder(fnelder,x0,dx0,500,1.0e-15,0);
	  }
      
	  public static double[] nelder(f_xj fnelder,double x0[],double tolerance)
      { 
        double [] dx0=new double[x0.length];
        for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
	      return nelder(fnelder,x0,dx0,tolerance);
	  }
      
	  public static double[] nelder(f_xj fnelder,double x0[],int maxiteration,double tolerance)
      { 
        double [] dx0=new double[x0.length];
        for(int i=0;i<x0.length;i++) dx0[i]=0.2*x0[i];
	      return nelder(fnelder,x0,dx0,maxiteration,tolerance,0);
	  }
	  
public static double[] newtond( fi_xi f,double x[])
{
// lineer olmıyan denklem siteminin kökleri
// fj(xi)=0 lineer olmıyan denklem sistemi için çözüm matrisi 
// türev matrisi dfj/dxi turev alt programıyla hesaplanır.
// Newton-Raphson metodu ile lineer olmıyan denklem sistemi çözümü 
//ti :önem katsayısı
//x bağımsız değişken vektörü
//y bağımlı değişken vektörü
//dy bağımlı değişkenin türev fonksiyonları vektörü
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-11;
int n=x.length;
double b[];
b=new double[n];
double dy[][];
dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
  for(ii=0;ii<n;ii++)
  {
    for(jj=0;jj<n;jj++)
    {
        dy[ii][jj]=dfx(f,x,ii,jj);
    }
  }
    b=multiply(divide(dy(f,x),H(f,x)),-ti);
    x=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 dfx(fi_xi f,double x[],int denklem_ref,int x_ref)
{ // Türev hesabı
double h=0.00011354523465;
double hh=1.0/h;
int i,m;
int n=x.length;
double f1[];
f1=new double[n];
double f2[];
f2=new double[n];
double f3[];
f3=new double[n];
double f4[];
f4=new double[n];

double x1[];
x1=new double[n];
double x2[];
x2=new double[n];
double x3[];
x3=new double[n];
double x4[];
x4=new double[n];

for(i=0;i<n;i++)
{
x1[i]=x[i];
x2[i]=x[i];
x3[i]=x[i];
x4[i]=x[i];
}
x1[x_ref]+=2.0*h;
x2[x_ref]+=h;
x3[x_ref]-=h;
x4[x_ref]-=2.0*h;

f1=f.func(x1);
f2=f.func(x2);
f3=f.func(x3);
f4=f.func(x4);
double fx=(-f1[denklem_ref]+8.0*f2[denklem_ref]-8.0*f3[denklem_ref]+f4[denklem_ref])/12.0*hh;
//System.out.println("f("+x_ref+","+denklem_ref+") = "+fx);
return fx;
}

public static double[][] H(fi_xi 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]=dfx(f,x,ii,jj);
    }
  }
  System.out.println(Matrix.toString(dy));
  return dy;
}

public static double[] dy(fi_xi f,double x[])
{
//birinci türev 
  int n=x.length;	
  double deltay[]=f.func(x);
  System.out.println(Matrix.toStringT(deltay));
  return deltay;
}

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[] verigir()
{
    String s=JOptionPane.showInputDialog("kök tahmin değerlerini arada bir boşluk bırakarak giriniz: ");        
    StringTokenizer token=new StringTokenizer(s);
    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[])  
     {
	 fxa f1 =new fxa();
     double [] x0=verigir();
     String s2="Newton-Raphson kök bulma : ";
     double [] r1= nelder_newton(f1,x0);
     String s=" kök değerleri  : \n"+Matrix.toStringT(r1);
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     } 
}
