import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.event.*;
import javax.swing.*;
import java.io.*;

abstract class f_x
{
  // single function  multi independent variable
  // a single value is returned  (
  // example f=x[0]+sin(x[1])
  // func(x) returns the value of f
  abstract public double func(double x);
}

abstract class f_xj
{ 
  // single function  multi independent variable
  // a single value is returned  (
  // example f=x[0]+sin(x[1])
  // func(x) returns the value of f
  abstract public double func(double x[]);
}

class f1dim extends f_x
{  
   // iki boyutlu fonksiyonu tek boyutlu fonksiyona çeviren
   // pseudo fonksiyon 
   f_xj ff1;
   double pc[];
   double xc[];
   double xt[];
   int n;	 
   public f1dim(f_xj ff,double pi[],double xi[])
   { n=pi.length;
     ff1=ff;	 
	 pc=new double[n];
     xc=new double[n];
     xt=new double[n];     
   	 for(int i=0;i<n;i++)
   	   {pc[i]=pi[i];xc[i]=xi[i];}  
   }

   public double func(double x)
   {
	for (int j=0;j<n;j++) xt[j]=pc[j]+x*xc[j];
	return ff1.func(xt); 	
   }   
}

      class yy extends f_xj
      {  
	  double xi[]; // independent variable data set
      double yi[]; //dependent variable data set
      double a[];  //fit function coefficient set
      int nn;
	  public yy(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 t,double ai[])
      {
      //sample specific function 
	  // saturation pressure as a function of saturation temperature
      a=ai;
	  double T=t+273.15;
      double pi=ai[0]+ai[1]/T+ai[2]/T/T+ai[3]*Math.log10(T)+ai[4]*T
      +ai[5]*T*T+ai[6]*T*T*T+ai[7]*((ai[8]-T)/T)*Math.log10(ai[8]-T);
      //pi=Math.pow(10.0,pi);
      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;
	  }
}

class zz extends yy
{     public zz(String filename,double ia[])
      {super(filename,ia);}
      public double func(double ai[])
      {
      return -super.func(ai);
	  }
}

public class opt
{
// ==================================================
// bir boyutlu lineer olmıyan fonksiyonların minimumu
// ==================================================

//1. golden (altın) oran - Fibonacchi metodu
//________________________________________________________________
public static double golden(f_x f,double a,double b,double epsilon)
{return golden(f,a,b,epsilon,1.0e-10,0);}

public static double golden(f_x f,double a,double b)
{return golden(f,a,b,1.0e-10,1.0e-10,0);}

public static double golden(f_x f,double a,double b,double epsilon,double delta,int print)
{
//  f fonksiyonunun minimumunu bulur
// not  maximum f(x) = minimum (-f(x))
    double r1 = (Math.sqrt(5.0)-1.0)/2.0; // golden ratio
    double r2 = r1*r1;
    double h = b - a;
    double ya = f.func(a);
    double yb = f.func(b);
    double c = a + r2*h;
    double d = a + r1*h;
    double yc = f.func(c);
    double yd = f.func(d);
    int k = 1;
    double dp,dy,p,yp;
    while ((Math.abs(yb-ya)>epsilon) || (h>delta))
    {
      k++;
      if (yc<yd)
        {
        b = d;
        yb = yd;
        d = c;
        yd = yc;
        h = b - a;
        c = a + r2 * h;
        yc = f.func(c);
        }
      else
       {
         a = c;
         ya = yc;
         c = d;
         yc = yd;
         h = b - a;
         d = a + r1 * h;
         yd = f.func(d);
       }//end of if
    }//end of while
    dp = Math.abs(b-a);
    dy = Math.abs(yb-ya);
    p = a;
    yp = ya;
    if (yb<ya)
    {
      p = b;
      yp = yb;
    }
    if(print==1)
    {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
    return p;
}

//2. ikinci derece polinom metodu
//________________________________________________________________ 
public static double polinom2derece(f_x f,double x0,double x1,double x2,double epsilon,double delta,int print)
{  
int maxit=100;    
double f0 = f.func(x0);
    double f1 = f.func(x1);
    double f2 = f.func(x2);
    double f3 = f2;
    double x3 = 0;
    double h = x1 - x0;
    double k=1;
    double dd=Math.abs(f1-f0);
    while ((dd >epsilon) || (h>delta))
    {
      k++;
      x3=(f0*(x1*x1-x2*x2)+f1*(x2*x2-x0*x0)+f2*(x0*x0-x1*x1))/
             (2*f0*(x1-x2)+2.0*f1*(x2-x0)+2.0*f2*(x0-x1));
      f3=f.func(x3);      
      if(x3 >= x0 && x3<x1)
        {x2=x1;f2=f1;x1=x3;f1=f3;}
      else if(x3 >= x1 && x3 <x2) 
        {x0=x1;f0=f1;x1=x3;f1=f3; }
      dd=Math.abs(f1-f0);
      h=Math.abs(x1-x0);
      if(k>maxit) break;          
    }//end of while
    if(print==1)
    {System.out.println("x  = "+x3+"f = "+f3+"hata x ="+h+"hata f(x)=y "+dd); }
    return x3;
}

//2a polinom quadrature metodu
//_________________________________________________________

  public static double quadmin(f_x f,double a,double b,double epsilon)
  {return quadmin(f,a,b,epsilon,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b)
  {return quadmin(f,a,b,1.0e-10,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b,double epsilon,double delta,int print)
  {
      // Minimum of a single variable function by using polynomial quadratures
      double p0 = a;
      int maxj = 20;
      int maxk = 30;
      double big = 1e6;
      double err = 1;
      int k = 1;
      int cond = 0;
      double  h = 1;
      double f1,p1,p2,pmin,y0,y1,y2,ymin,d,hmin,dh,h0,h1,h2;
      double e0,e1,e2;
      double p,dp,dy,yp;
      if (Math.abs(p0)>1e4)
        { h = Math.abs(p0)/1e4;}
      while (k<maxk && err>epsilon && cond !=5 )
      {
        f1 = (f.func(p0 + 0.00001) - f.func(p0 - 0.00001)) / 0.00002;
        if (f1 > 0) {h = -Math.abs(h);}
        p1 = p0 + h;
        p2 = p0 + 2 * h;
        pmin = p0;
        y0 = f.func(p0);
        y1 = f.func(p1);
        y2 = f.func(p2);
        ymin = y0;
        cond = 0;
        int j = 0;
        while (j < maxj && Math.abs(h) > delta && cond == 0)
        {
          if (y0 <= y1)
          {
            p2 = p1;
            y2 = y1;
            h = h / 2;
            p1 = p0 + h;
            y1 = f.func(p1);
          }
          else
          {
            if (y2 < y1)
            {
              p1 = p2;
              y1 = y2;
              h = 2 * h;
              p2 = p0 + 2 * h;
              y2 = f.func(p2);
            }
            else
            {  cond = -1;}
          } //end of if(y0 <= y1)
          j = j + 1;
        if (Math.abs(h) > big || Math.abs(p0) > big) {cond = 5;}
        }//end of while(j<maxj
        if (cond == 5)
        {
              pmin = p1;
              ymin = f.func(p1);
        }
        else
        {
          d = 4 * y1 - 2 * y0 - 2 * y2;
          //Start of a long block:
          if (d < 0) {
            hmin = h * (4 * y1 - 3 * y0 - y2) / d;
          }
          else
          {
            hmin = h / 3;
            cond = 4;
          }
          pmin = p0 + hmin;
          ymin = f.func(pmin);
          dh = Math.abs(h);
          h0 = Math.abs(hmin);
          h1 = Math.abs(hmin - h);
          h2 = Math.abs(hmin - 2 * h);
          if (h0 < dh) {h = h0;}
          if (h1 < dh) {h = h1;}
          if (h2 < dh) {h = h2;}
          if (h == 0) {h = Math.abs(hmin);}
          if (h < delta) {cond = 1;}
          if (Math.abs(h) > big || Math.abs(pmin) > big) {cond = 5;}
          e0 = Math.abs(y0 - ymin);
          e1 = Math.abs(y1 - ymin);
          e2 = Math.abs(y2 - ymin);
          if (e0 != 0 && e0 < err) {err = e0;}
          if (e1 != 0 && e1 < err) {err = e1;}
          if (e2 != 0 && e2 < err)  {err = e2;}
          if (e0 != 0 && e1 == 0 && e2 == 0) {err = 0;}
          if (err < epsilon) {cond = 2;}
          p0 = pmin;
          k = k + 1;
      } // End of the long else block
      if (cond == 2 && h < delta) {cond = 3;}
  }//end of the while
  p = p0;
  dp = h;
  yp = f.func(p);
  dy = err;
  if(print==1)
  {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
  return p;
  }
  

// 3. brent metodu 
//________________________________________________________________
// golden (altın) oran - Fibonacchi metodu ve ikinci derece polinom metodunu birlikte kullanır

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 df(f_x f,double x)
{
// türev alt programı	
// df fonksiyonunun türevini hesaplar
double h0=0.0256;
int i,m;
int n=7;
//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++)
{
T[i][0]=( f.func(x + h[i]) - f.func( x - h[i]))/(2.0*h[i]);
}
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;
}

//4. Türevli brent metodu
//________________________________________________________________
public  static double dbrent(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;
// MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
// {a=d;b=e;c=f}
	int iter;
	boolean ok1,ok2;
	double a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
	d=0.0;
	double fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
	double xmin=bx;
	a=(ax < cx ? ax : cx);
	b=(ax > cx ? ax : cx);
	x=w=v=bx;
	fw=fv=fx=f.func(x);
	dw=dv=dx=df(f,x);
	for (iter=1;iter<=ITMAX;iter++) {
		xm=0.5*(a+b);
		tol1=tol*Math.abs(x)+ZEPS;
		tol2=2.0*tol1;
		if (Math.abs(x-xm) <= (tol2-0.5*(b-a))) {
			xmin=x;
			return xmin;
		}
		if (Math.abs(e) > tol1) {
			d1=2.0*(b-a);
			d2=d1;
			if (dw != dx)  d1=(w-x)*dx/(dx-dw);
			if (dv != dx)  d2=(v-x)*dx/(dx-dv);
			u1=x+d1;
			u2=x+d2;
			ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
			ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
			olde=e;
			e=d;
			if (ok1 || ok2) {
				if (ok1 && ok2)
					d=(Math.abs(d1) < Math.abs(d2) ? d1 : d2);
				else if (ok1)
					d=d1;
				else
					d=d2;
				if (Math.abs(d) <= Math.abs(0.5*olde)) {
					u=x+d;
					if (u-a < tol2 || b-u < tol2)
						d=SIGN(tol1,xm-x);
				} else {
					d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
				}
			} else {
				d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
			}
		} else {
			d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
		}
		if (Math.abs(d) >= tol1) {
			u=x+d;
			fu=f.func(u);
		} else {
			u=x+SIGN(tol1,d);
			fu=f.func(u);
			if (fu > fx) {
				xmin=x;
				return xmin;
			}
		}
		du=df(f,u);
		if (fu <= fx) {
			if (u >= x) a=x; else b=x;
			//MOV3(v,fv,dv, w,fw,dw)
			{v=w;fv=fw;dv=dw;}
			//MOV3(w,fw,dw, x,fx,dx)
			{w=x;fw=fx;dw=dx;}
			//MOV3(x,fx,dx, u,fu,du)
			{x=u;fx=fu;dx=du;}
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				//MOV3(v,fv,dv, w,fw,dw)
				{v=w;fv=fw;dv=dw;}
				//MOV3(w,fw,dw, u,fu,du)
				{w=u;fw=fu;dw=du;}
			} else if (fu < fv || v == x || v == w) {
				//MOV3(v,fv,dv, u,fu,du)
				{v=u;fv=fu;dv=du;}
			}
		}
	}
	System.out.println("DBRENT maximum iterasyon sayısını aştı sonuç geçerli olmıyabilir");
	return xmin;
}

// 5. minumumu brakete alma (sınırlarını bulma) arama metodu
//    altın oran araması kullanır

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;
}
//
// ==================================================
// yardımcı metodlar türevler vs..
// ==================================================

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]=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[][]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]=turev(f,x,ii,jj);
    }
  }
  return dy;
}


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;	
}


// ==================================================
// çok boyutlu lineer olmıyan fonksiyonların minimumu
// ==================================================
//

// 1. Newton metodu
//________________________________________________________________
public static double[] newtond( f_xj 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-10;
int n=x.length;
double b[];
b=new double[n];
double dy[][]=H(f,x);
//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++)
  //  {   //Hessian matrisi
  //      dy[ii][jj]=turev(f,x,ii,jj);
  //  }
  //System.out.println(Matrix.toString(dy));
  //}
    b=Matrix.multiply(Matrix.divide(df(f,x),dy),-ti);
    x=Matrix.add(x,b);
    System.out.println("i= "+i+"\n" );
    for(int kk=0;kk<x.length;kk++)
        System.out.println(x[kk]+" ");
}
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 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;
}
//2. Nelder ve Mead metodu (simpleks arama metodu)
//________________________________________________________________
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,500,1.0e-20,0);}

          

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;   	   
}
// 3. en dik yamac metodu
//________________________________________________________________
public static double[] endikyamac( f_xj f,double x[])
{
// en dik yamac metodu
// linmin ve brent metodunukullanır	
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-12;
int n=x.length;
double b[];
b=new double[n];
double lambda;
double dy[]=new double[n];
//dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    dy=turev(f,x);
    lambda=linmin(f,x,dy);
    b=Matrix.multiply(dy,lambda);
    x=Matrix.add(x,b);
    System.out.println("i= "+i+"\n"+Matrix.toString(x));
}
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;
}

// 4. modifiye Newton metodu (en dik yamaç ve newton metodlarını birleştirir)
//________________________________________________________________
public static double[] dampednewton( f_xj f,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=500;
double tolerance=1.0e-12;
int n=x.length;
double b[];
b=new double[n];
double lambda;
double g[]=new double[n];
double G[][]=H(f,x);
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=turev(f,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;
}

// 5. Powell metodu
//________________________________________________________________
public static double SQR(double a) {return a*a;}

public static double[] powell(f_xj f,double p[])
{  
	// 	
	int n=p.length;
	double xi[][]=Matrix.I(n);
	double fret;
	double ftol=1e-12;
	double sqrarg;
	int ITMAX=200;
	int i,ibig,j;
	double t,fptt,fp,del;
	double [] pt,ptt,xit;
	pt=new double[n];
	ptt=new double[n];
	xit=new double[n];
	double aa[][]=new double[2][n];
	int iter;
	fret=f.func(p);
	for (j=0;j<n;j++) pt[j]=p[j];
	for (iter=0;iter<=ITMAX;iter++) 
	{
		fp=fret;
		ibig=0;
		del=0.0;
		for (i=0;i<n;i++) {
			for (j=0;j<n;j++) xit[j]=xi[j][i];
			fptt=fret;
			aa=linminiter(f,p,xit);
			p=aa[0];
			xit=aa[1];
			fret=f.func(xit);
			if (Math.abs(fptt-fret) > del) {
				del=Math.abs(fptt-fret);
				ibig=i;
			}
		}
		if (2.0*Math.abs(fp-fret) <= ftol*(Math.abs(fp)+Math.abs(fret))) {
			break;
		}
		if (iter == ITMAX) System.out.println("Too many iterations in routine POWELL");
		for (j=0;j<n;j++) {
			ptt[j]=2.0*p[j]-pt[j];
			xit[j]=p[j]-pt[j];
			pt[j]=p[j];
		}
		fptt=f.func(ptt);
		if (fptt < fp) {
			t=2.0*(fp-2.0*(fret)+fptt)*SQR(fp-(fret)-del)-del*SQR(fp-fptt);
			if (t < 0.0) {
				aa=linminiter(f,p,xit);
			    p=aa[0];
			    xit=aa[1];
			    fret=f.func(xit);
				for (j=0;j<n;j++) xi[j][ibig]=xit[j];
			}
		}
	}
return p;
}

// 7.  Davidon - Fletcher - Powell metodu
//________________________________________________________________


public static double[] dfpmin( f_xj f,double p[])
{
	double ftol=1.0e-10;
	int n=p.length;
	int ITMAX=200;
    double EPS=1.0e-10;
	int j,i,its;
	double fp,fae,fad,fac;
	double hessin[][]=new double[n][n];
    double xi[]=new double[n];
	double g[]=new double[n];
	double dg[]=new double[n];
	double hdg[]=new double[n];
	fp=f.func(p);
    double aa[][]=new double[2][n];
	g=turev(f,p);
	int iter=0;
	double fret;
	for (i=0;i<n;i++) 
	{
		for (j=0;j<n;j++) hessin[i][j]=0.0;
		hessin[i][i]=1.0;
		xi[i] = -g[i];
	}
	for (its=0;its<ITMAX;its++) {
		iter=its;
		aa=linminiter(f,p,xi);
		p=aa[0];
		xi=aa[1];
		fret=f.func(xi);
		if (2.0*Math.abs(fret-fp) <= ftol*(Math.abs(fret)+Math.abs(fp)+EPS)) {
			return p;
		}
		fp=fret;
		for (i=0;i<n;i++) dg[i]=g[i];
		fret=f.func(p);
		g=turev(f,p);
		for (i=0;i<n;i++) dg[i]=g[i]-dg[i];
		for (i=0;i<n;i++) {
			hdg[i]=0.0;
			for (j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
		}
		fac=fae=0.0;
		for (i=0;i<n;i++) {
			fac += dg[i]*xi[i];
			fae += dg[i]*hdg[i];
		}
		fac=1.0/fac;
		fad=1.0/fae;
		for (i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
		for (i=0;i<n;i++)
			for (j=0;j<n;j++)
				hessin[i][j] += fac*xi[i]*xi[j]
					-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
		for (i=0;i<n;i++) {
			xi[i]=0.0;
			for (j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
		}
	}
	System.out.println("Too many iterations in DFPMIN");
	return p;
}

public static double[] davidon_fletcher_powell( f_xj f,double x[])
{
//  Davidon - Fletcher - Powell metodu
//  
// linmin ve brent metodunu kullanır	
double ti=-1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-12;
int n=x.length;
double b[];
b=new double[n];
double lambda;
double g[]=new double[n];
double gold[]=new double[n];
double xold[]=new double[n];
double H[][]=new double[n][n];
double Hold[][]=Matrix.I(n);
double delta[]=new double[n];
double gamma[]=new double[n];
double g1[][]=new double[n][n];
double g3[][]=new double[n][n];
double g2,g4,g5;
double aa[][]=new double[2][n];
i=0;
    //ilk iterasyon olarak en dik yamaç metodunu kullan
    g=turev(f,x);
    for(i=0;i<n;i++)
    {
    b[i]=1.0;
    xold[i]=x[i];
    gold[i]=g[i];
    delta[i]=0.0;
    gamma[i]=0.0;
    }
    aa=linminiter(f,x,g);
	x=aa[0];
	b=aa[1];    
    System.out.println("Horiginal=\n"+Matrix.toString(Matrix.inverse(H(f,x))));
while( i++ < nmax && Matrix.abs(b) > tolerance )
{   
	g=turev(f,x);
    for(i=0;i<n;i++)
    { 
    delta[i]=x[i]-xold[i];
    gamma[i]=g[i]-gold[i];
    }
    g1=Matrix.V_VT(delta);//iki boyut
    g2=Matrix.VT_X(g,gamma);
    g1=Matrix.divide(g1,g2);//iki boyut
    g3=Matrix.V_VT(Matrix.multiply(Hold,gamma));
    //g3=Matrix.multiply(g3,Matrix.Transpose(g3));
    g4=Matrix.VT_G_V(gamma,Hold);
    g3=Matrix.divide(g3,g4);
    H=Matrix.substract(Matrix.add(Hold,g1),g3);
    System.out.println("H=\n"+Matrix.toString(H));
    lambda=-linmin(f,x,g);
    b=Matrix.multiply(Matrix.multiply(H,g),lambda);
    x=Matrix.add(x,b);
    for(i=0;i<n;i++)
    {
       for(int j=0;j<n;j++)
         {Hold[i][j]=H[i][j];}
       xold[i]=x[i];
       gold[i]=g[i];  
    } 
}
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;
}

// 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; 
     }
}

