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

class SVDLeastsquare
{

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[][] T(double [][] left)
{return Transpose(left);}

public static double[][] I(int n)
{
//unit matrix
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=1.0;
return b;
}
//Converts vector to diagonal values of matrix
public static double[][] S(double[] left,int m)
{
int n=left.length;
double b[][];
b=new double[m][n];
for(int i=0;i<n;i++)
    b[i][i]=left[i];   
return b;
}

//Converts double to diagonal values of matrix
public static double[][] S(double left,int n)
{
double b[][];
b=new double[n][n];
for(int i=0;i<n;i++)
    b[i][i]=left;
return b;
}

//takes invert value of a diagonal matrix
public static double[][] invS(double[][] S)
{
int m=S.length;
int n=S[0].length;
double b[][];
b=new double[m][n];
for(int i=0;i<n;i++)
    {  if(S[i][i]==0) b[i][i]=0; 
	   else b[i][i]=1.0/S[i][i];
	}
return b;
}

//square root of a diagonal matrix
public static double[][] sqrtS(double[] left,int m)
{
int n=left.length;
double b[][];
b=new double[m][n];
for(int i=0;i<n;i++)
    b[i][i]=Math.sqrt(left[i]);
return b;
}
//multiply two matrices
public static double[][] multiply(double[][] left,double[][] right)
{
//multiplication of two matrices
int ii,jj,i,j,k;
int m1=left[0].length;
int n1=left.length;
int m2=right[0].length;
int n2=right.length;
double[][] b;
b=new double[n1][m2];
 if(m1 !=n2 )
 {
 System.out.println("inner matrix dimensions must agree"+"m1="+m1+"n2="+n2);
 return b;
 }
 for(i=0;i<n1;i++)
   {
   for(j=0;j<m2;j++)
     {
     for(k=0;k<m1;k++)
	    {b[i][j]+=left[i][k]*right[k][j];
	    }
     }
   }
return b;
//end of multiply of two matrices
}

//multiply a matrix with a double value
public static double[][] multiply(double[][] left,double right)
{
//multiplying a matrix with a constant
int i,j;
int n=left.length;
int m=left[0].length;
double b[][];
b=new double[n][m];
for(i=0;i<n;i++)
  {
  for(j=0;j<m;j++)
    b[i][j]=left[i][j]*right;
  }
return b;
//end of multiplying a matrix with a constant double
}

public static double[] multiply(double[][] left,double[] right)
{
//multiplication of one matrix with one vector
int ii,jj,i,j,k;
int n1=left.length;
int m1=left[0].length;

int m2=right.length;
double[] b;
b=new double[m2];
 for(i=0;i<n1;i++)
   {
     b[i]=0;
     for(k=0;k<m1;k++)
	 {b[i]+=left[i][k]*right[k];
	 }
   }
return b;
//end of multiply of a matrix and a vector
}

//substract two matrices
public static double[][] substract(double[][] left,double[][] right)
{
//addition of two matrices
int n1=left.length;
int m1=left[0].length;
int n2=right.length;
int m2=right[0].length;
int nMax,mMax;
int i,j;
if(m1>=m2) mMax=m1;
else       mMax=m2;
if(n1>=n2) nMax=n1;
else       nMax=n2;
double b[][];
b=new double[nMax][mMax];
for(i=0;i<n1;i++)
  {
  for(j=0;j<m1;j++)
    {
    b[i][j]=b[i][j]+left[i][j];
    }
  }
for(i=0;i<n2;i++)
  {
  for(j=0;j<m2;j++)
    {
    b[i][j]=b[i][j]-right[i][j];
    }
  }
return b;
//end of matrix substraction method
}
 
public static double[] SVD_LSQ(double xi[],double yi[],int n_leastsquare)
{ // Singular value decomposition least square approximation
  double U[][];
  double V[][];
  double S[][];
  double invS[][];
  double A[][]=A(xi,n_leastsquare);
  int m=A.length;
  int m1=A[0].length;
  System.out.println("A="+Matrix.toString(A));
  double ATA[][]=multiply(Transpose(A),A);
  System.out.println("ATA="+Matrix.toString(ATA));
  //double aa[][]=Householder(ATA);
  double cc[][];
  cc=eigen_jacobi(ATA);//eigenQL(aa);
   System.out.println("c="+Matrix.toString(cc));
  int n=cc[0].length;
  int n1=cc.length;
  double dummy[]=new double[n1];
  for(int i=0;i<n;i++)
  {for(int k=i+1;k<n;k++)
    {if(cc[0][i]<cc[0][k])
      { for(int j=0;j<n1;j++)
        { dummy[j]=cc[j][k];
          cc[j][k]=cc[j][i];
          cc[j][i]=dummy[j];
        }
      }
    }
  }
  double x[]=new double[n];
  for(int i=0;i<n;i++)
  {x[i]=cc[0][i];}
  S=sqrtS(x,m1);
  System.out.println("S="+Matrix.toString(S));
  invS=invS(S);
  System.out.println("invS="+Matrix.toString(invS));
  double xx[][]=new double[n][n];
  V=new double[n][n];
  U=new double[n][n];
  double L[]=new double[n];
 
  for(int i=0;i<n;i++)
  { for(int j=0;j<n;j++)
    {xx[i][j]=cc[i+1][j];}
  }

  double bb[]=new double[n];  
  for(int j=0;j<n;j++)
  { for(int i=0;i<n;i++)
    {L[j]+=xx[i][j]*xx[i][j];}
    L[j]=Math.sqrt(L[j]);
  }
  for(int i=0;i<n;i++)
  { for(int j=0;j<n;j++)
    {V[i][j]=xx[i][j]/L[i];}
  }
  System.out.println("VT="+Matrix.toString(T(V)));
  U=multiply(A,multiply(V,invS)); 
   System.out.println("U="+Matrix.toString(U));
  double c[]=multiply(T(U),yi);
   System.out.println("c="+Matrix.toString(c)); 
   double z[]=new double[n];
   for(int i=0;i<n;i++) z[i]=c[i]/S[i][i];
   System.out.println("z="+Matrix.toString(z)); 
   double a[]=multiply(V,z);
   System.out.println("a="+Matrix.toString(a)); 
   return a;
}
public static double func(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=e.length;
double ff;
if(n!=0.0)
  { ff=e[n-1];
  for(int i=n-2;i>=0;i--)
   { ff=ff*x+e[i]; }
  }
else
  ff=0;
return ff;
}

public static double error(double x[],double y[],double e[])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 for(k=0;k<n;k++)
 {
 total+=(y[k]-func(e,x[k]))*(y[k]-func(e,x[k]));
 }
total=Math.sqrt(total);
return total;
}

public static double[][]  funcSVD_LSQ(double xi[],double yi[],int polinomkatsayisi,int aradegersayisi)
 {
	//aradegersayisi:  x--o--o--x--o--o--x zincirinde x deneysel noktalar ise
	// ara değer sayısı 2 dir
	 int n=xi.length;
	int nn=(n-1)*(aradegersayisi+1)+1;
	double z[][]=new double[2][nn];
    double E[]=SVD_LSQ(xi,yi,polinomkatsayisi);	
    System.out.println("coefficients :\n"+Matrix.toStringT(E));
    double dx=0;
    int k=0;
    int i;
    for(i=0;i<(n-1);i++)
       {z[0][k]=xi[i];z[1][k]=func(E,z[0][k]);k++;
       for(int j=0;j<aradegersayisi;j++) 
       {dx=(xi[i+1]-xi[i])/((double)aradegersayisi+1.0);
       z[0][k]=z[0][k-1]+dx;z[1][k]=func(E,z[0][k]);k++;}
       }
       z[0][k]=xi[i];z[1][k]=func(E,z[0][k]);
    return z;
 } 
  public static void print(String s)
  {JOptionPane.showMessageDialog(null,s,
     "Singular value decomposition : ",JOptionPane.PLAIN_MESSAGE);}
  
  public static double[][] eigen_jacobi(double A[][])
  {	  //Eigen value calculation
   double eps=1e-15;
   double err=1;
   int m=100;
   int n=A.length;
   double P[][]=new double[n][n];   
   double R[][]=new double[n][n];   
   int p=n,q=n;
   double max=-9.99e100;
   double fi;
   P=I(n);
   for(int i=0;i<n;i++)
   {for(int j=i+1;j<n;j++)
       {if(Math.abs(A[i][j]) > max) {max=Math.abs(A[i][j]);p=i;q=j;}}
   }  
   while(Math.abs(A[p][q]) > eps )
   {
    fi=0.5*Math.atan2(2.0*A[p][q],(A[p][p]-A[q][q]));
    R=I(n);
    R[p][p]=Math.cos(fi);
    R[q][q]=Math.cos(fi);
    R[p][q]=-Math.sin(fi);
    R[q][p]=Math.sin(fi);
    P=multiply(P,R);
    A=multiply(T(R),multiply(A,R));
    max=-9.99e100;
    for(int i=0;i<n;i++)
   {for(int j=i+1;j<n;j++){if(Math.abs(A[i][j])>max) {max=Math.abs(A[i][j]);p=i;q=j;}}}  
    }
   double lambda[]=new double[n];
   for(int i=0;i<n;i++) {lambda[i]=A[i][i];}
   double PP[][]=new double[n+1][n];
   PP[0]=lambda;
   for(int i=0;i<n;i++)
   {for(int j=0;j<n;j++)
     {PP[i+1][j]=P[i][j];}
   }
   return PP;
  }
  
  public static double[] gausswithpartialpivot(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++)
  { //partial 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;
	}
	//  
    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];
	}
  }
  // backward 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[][] A(double xi[],int n_leastsquare)
{
int n=xi.length;
int m=n_leastsquare+1;
int l=xi.length;
int i,j,k;
double A[][];
A=new double[n][m];
double B[];
B=new double[n];
double X[];
X=new double[n];
double cc=1.0;
for(i=0;i<n;i++)
  { 
     cc=1.0;
	 for(j=0;j<m;j++)
     {
	 A[i][j]=cc;
	 cc*=xi[i];    
     }
  }  
return A;
}

  public static void main(String args[]) throws IOException
  {
  double x[];
  double y[];
  String s1=JOptionPane.showInputDialog("Data file name (a.txt)");
  // JFileChooser fc=new JFileChooser();
  // if (fc.showOpenDialog(null) == JFileChooser.APPROVE_OPTION) {File file = fc.getSelectedFile();s1=file.getName(); }
  double a[][]=Text.readDoubleT(s1);
  x=a[0];
  y=a[1];
  int n=Integer.parseInt(JOptionPane.showInputDialog("polinom derecesi (n)"));
  double z[][]=funcSVD_LSQ(x,y,n,4);
  System.out.println("Polynomial Least Square Curve fitting\n"+Matrix.toStringT(z)); 
  Plot pp=new Plot(a[0],a[1]);
  pp.setPlabel("Polynomial Singular Value Decomposition Least Square Curve fitting");
  pp.setXlabel("x");
  pp.setYlabel("y=f(x)");
  pp.setPlotType(0,22); 
  pp.addData(z[0],z[1]);
  pp.setPlotType(1,0); 
  pp.setGrid(1,1);
  pp.plot();
  }
}