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

//  Surface fitting by using moving least square model
//  f(x,y)=a0+a1*x+a2*y+a3*x^2+a4*y^2+a5*xy+a6*x^3+a7*y^3+a8*x^2y
//  +a9*xy^2+a10*x^4+a11*y^4+a12*x^2y^2+a13x^3y+a14*xy^3

abstract class f_xir
{ abstract double func(double x[],int equation_ref);}

//You can define your own fitting function
class fax extends f_xir
{	
	double func(double x[],int i)
	{//A general two dimensional surface polynomial equation
	 // of order (up to) 14 	
    double xx=0.0;
    //density function 		
    if(i==0)       {xx=1.0;           }
    else if (i==1) {xx=x[0];          }
    else if (i==2) {xx=x[1];          }
    else if (i==3) {xx=x[0]*x[0];     }
    else if (i==4) {xx=x[1]*x[1];     }
    else if (i==5) {xx=x[0]*x[1];     }
    else if (i==6) {xx=x[0]*x[0]*x[0];}
    else if (i==7) {xx=x[1]*x[1]*x[1];}
    else if (i==8) {xx=x[0]*x[0]*x[1];}
    else if (i==9) {xx=x[0]*x[1]*x[1];}
    else if (i==10){xx=x[0]*x[0]*x[0]*x[0];}
    else if (i==11){xx=x[1]*x[1]*x[1]*x[1];}
    else if (i==12){xx=x[0]*x[0]*x[1]*x[1];}
    else if (i==13){xx=x[0]*x[0]*x[0]*x[1];}
    else if (i==14){xx=x[0]*x[1]*x[1]*x[1];}
    else if (i==15){xx=x[0]*x[0]*x[0]*x[0]*x[0];}
    else if (i==16){xx=x[1]*x[1]*x[1]*x[1]*x[1];}
    else if (i==17){xx=x[0]*x[1]*x[1]*x[1]*x[1];}
    else if (i==18){xx=x[0]*x[0]*x[0]*x[0]*x[1];}
    else if (i==19){xx=x[0]*x[0]*x[1]*x[1]*x[1];}
    else if (i==20){xx=x[0]*x[0]*x[0]*x[1]*x[1];}    
    return xx;
	}
}    

class NA129H
{
	public static double d(double x[],double xi[])
	{// Euler distance
	double total=0;
	double dx=1;
	for(int i=0;i<xi.length;i++)
	{   dx=x[i]-xi[i];
		total+=dx*dx;}
	return Math.sqrt(total);	
    }
    
    public static double w(double x[],double xi[])
    {
	// Moving least square weight function
	double d=d(x,xi);
	double eps=0.001;
	double alpha=1.0/16.0;
	double ww=Math.exp(-alpha*d*d)/(d*d+eps);
	return ww;
	}
    	
	public static double[][] Transpose(double [][] left)
    { //transpose matrix (if A=a(i,j) Transpose(A)=a(j,i)
    int i,j;
    int n=left.length;
    int m=left[0].length;
    double b[][];
    b=new double[m][n];
    for(i=0;i<n;i++)
      {for(j=0;j<m;j++)  {b[j][i]=left[i][j];} }
    return b;
    }  
    
  public static double[] pivotlugauss(double a[][],double b[])
  { //gauss elimination with partial pivoting
  int n=b.length;	  
  double x[]=new double[n];
  double carpan=0;
  double toplam=0;
  double buyuk;
  double dummy=0;
  //gauss elimination
  int i,j,k,p,ii,jj;
  for(k=0;k<(n-1);k++)
  { //pivoting
	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 elimination process 
    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];
	}
  }
  //back substitution
  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 HEKK(double x[],double xi[][],double yi[],int n)
{ // n dimension surface curve fitting for the vector point x[]
fax f=new fax();
int l=xi.length; //number of independent variables
int i,j,k;
int np1=n+1;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double E[];
E=new double[np1];
for(i=0;i<n+1;i++)
  {  B[i]=0;
	 for(j=0;j<np1;j++)
     {
	 A[i][j]=0.0;    
     for(k=0;k<l;k++) A[i][j]+= w(x,xi[k])*f.func(xi[k],i)*f.func(xi[k],j);
     }
     for(k=0;k<l;k++) B[i]+= w(x,xi[k])*f.func(xi[k],i)*yi[k];
  }  
E=pivotlugauss(A,B);
//X=B/A;
double max=0;
for(i=0;i<n+1;i++)
 if(Math.abs(E[i]) > max) max = Math.abs(E[i]);
for(i=0;i<n+1;i++)
 if((Math.abs(E[i]/max) > 0) && (Math.abs(E[i]/max) < 1.0e-100)) E[i]=0;
double z=funcEKKgenel(E,x);
return z;
}


public static double funcEKKgenel(double e[],double x[])
{
// Moving least square function evaluation
fax f=new fax();
int n=e.length;
double ff=0;
if(n!=0.0)
  { 
  for(int i=n-1;i>=0;i--)
   {ff+=e[i]*f.func(x,i);}
  } 
return ff;
}

public static double[] HEKK(double x[][],double xi[][],double yi[],int n)
{
	// n dimension surface curve fitting for the array of vector points x[]
int m=x.length;
double z[]=new double[m];
for(int i=0;i<m;i++)
{
z[i]= HEKK(x[i],xi,yi,n);
}	
return z;
}

public static double[][] output(double x[][],double xi[][],double yi[],int n)
{  //gives curve fitting values together with the independent value set
int m=x.length;
int nn=x[0].length;
double c[][]=new double[m][nn+1];
double z[]=new double[m];
for(int i=0;i<m;i++)
{for(int j=0;j<nn;j++)
	{c[i][j]=x[i][j];}
}	
for(int i=0;i<m;i++)
{
c[i][nn]= HEKK(x[i],xi,yi,n);
}	
return c;
}
 
  public static void main(String[] args) 
  {
	  /*
  String s1=JOptionPane.showInputDialog("curve fitting input data : ");	  
  double c[][]=Text.readDouble(s1); 
  int n=c[0].length;
  int m=n-1;
  int l=c.length;
  double x[][]=new double[l][m];
  double y[]=new double[l];
  for(int i=0;i<l;i++)
  {for(int j=0;j<m;j++)
   {x[i][j]=c[i][j];}
   y[i]=c[i][m];
  }
  String s2=JOptionPane.showInputDialog("required curve fitting independent variable list data : ");
  double xy[][]=Text.readDouble(s2);
  double z[][]=output(xy,x,y,4);  
  Text.print(z);
  f2a f=new f2a(c,3);
	    surfacePlot p=new surfacePlot(f,-5.0,5.0,-5.0,5.0,-10.0,10.0,"X","Y","Z",150);
        p.plot();
        */
        double a[][]=Text.readDouble("a.txt");
	    Text.print(a);
	    int n=a[0].length;
        int m=n-1;
        System.out.println("n="+n+"m="+m);
        int l=a.length;
	    double x[][]=new double[l][m];
        double y[]=new double[l];
        for(int i=0;i<l;i++)
        { for(int j=0;j<m;j++)
        { System.out.println("i="+i+"j="+j);
           x[i][j]=a[i][j];}
          y[i]=a[i][m];
        }
        Text.print(x);
        Text.printT(y);
	    double z[][]=output(a,x,y,14);  
        Text.print(z);
	    f2a f=new f2a(a,4);
	    surfacePlot p=new surfacePlot(f,-5.0,5.0,-5.0,5.0,-10.0,10.0,"X","Y","Z",50);
        p.plot();
  }
}
