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[]) 
{
//function to minimize
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];
return ff; 
}
}

public class NA43C
{

public static double VT_V(double [] left)
{
//multiplies 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[] substract(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 substraction method
}

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[] df(f_xj f,double xxi[])
{
// f(x0,x1,x2,...xn) derivative of the function
// [df/dx0, df/dx1,df/dx2,....,df/dxn]
// df/dxj  j=0...x.length
// 
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)
{ //second derivative
  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 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)
{
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("Maximum iteration in BRENT is exceeded");
	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)
{return mnbrak(f,ax,(ax-dx),(ax+dx));}

public  static double[] mnbrak(f_x f,double ax)
{return mnbrak(f,ax,(ax-0.1),(ax+0.1));}

public  static double[] mnbrak(f_x f,double ax,double bx,double cx)
{
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];
}
// 4. Polak Ribiere Method 
//________________________________________________________________
public static double[] polak_ribiere( f_xj f,double x[])
{	
int i,j,ii,jj,k;
int nmax=500;
double gamma=1.0;
double tolerance=1.0e-7;
int n=x.length;
double alpha;
double g[]=new double[n];
double ge[]=new double[n];
double d[]=new double[n];
double beta;
double alpha_d;
double Q[]=new double[n];
int nn=15;
i=0;
j=0;k=0;
while( j < nmax && gamma > tolerance )
{
    g=df(f,x);
    gamma=norm(g);
    if(k<1) {for(i=0;i<n;i++){d[i]=-g[i];};ge=g;}
    else    
    {   // Polack-Ribiere conjugate gradient definition
	    beta=VT_X(substract(ge,g),g)/VT_V(ge);
     for(i=0;i<n;i++){Q[i]=beta*d[i];};
    }
    for(i=0;i<n;i++){d[i]=-g[i]+Q[i];}
    alpha=linmin(f,x,d);
    for(i=0;i<n;i++){x[i]+=alpha*d[i];}
    ge=g;
    j++;k++;
    if(k>=nn) k=0; //re-apply a step of steepest descent
}
if(i >= nmax) JOptionPane.showMessageDialog(null,"Warning : Maximum number of iteration exceeded \n"+
  " result may not be valid","MAXIMUM ITERATION WARNING",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="enter optimum point first estimation values with one empty space in between: ";   
     double [] x0=verigir(s);
     f1 f_x=new f1();
     double [] r1= polak_ribiere(f_x,x0);
     s=" minimum values : \n"+Matrix.toStringT(r1);
     String s2="Polak-Ribiere multivariable non linear optimisation:";
     JOptionPane.showMessageDialog(null,s,s2,JOptionPane.PLAIN_MESSAGE);
     System.exit(0);
     }    
}
