/*****************************************************************
   JamaFit.java
      uses the NIST JAMA
   http://math.nist.gov/javanumerics/jama
//  matrix library to fit a least squares parabola 
//      y(x) = b0 + b1 x + b2 x^2
//  
//  The JAMA library must be in  same directory as program
//      else modify CLASSPATH variable to include  JAMA directory
//      uses  Matrix.class,   see Matrix.java
//  or the JAMA HTML documentation for how it works.
//   
//*****************************************************************/
    
    
import Jama.*;  
import java.io.*;
    
public class JamaFit      
{   
  public static void main(String[] argv)  
{   
    double []   x = {1., 1.05, 1.15, 1.32,  1.51, 1.68, 1.92};
    double []   y = {0.52, 0.73, 1.08, 1.44, 1.39, 1.46, 1.58}; 
    double [] sig = {0.1, 0.1, 0.2, 0.3, 0.2, 0.1, 0.1};
    int Nd = 7;      // Number of data points
    double sig2,  s,sx,sxx,sy,sxxx,sxxxx,sxy,sxxy, rhl;
  //----------------------------------------------------
  // Create 3 3 A matrix, 3x1 X vector and 3x1 B vector
  //----------------------------------------------------
        
    double [][] Sx = new double[3][3];
    double [] [] Sy = new double[3][1]; 
    
  // Generate matrix elements

     s= sx = sxx = sy = sxxx = sxxxx = sxy = sxy = sxxy = 0;
     for (int i=0; i <= Nd-1; i++)
    {
     sig2   = sig[i]*sig[i];
         s += 1./sig2;
     sx    += x[i]/sig2;
     sy    += y[i]/sig2;
     rhl    = x[i]*x[i];
     sxx   += rhl/sig2;
     sxxy  += rhl*y[i]/sig2;
     sxy   += x[i]*y[i]/sig2;
     sxxx  += rhl*x[i]/sig2;
     sxxxx += rhl*rhl/sig2;
    }
    
  // Load arrays Sx & Sy

     Sx[0][0] = s;
     Sx[0][1] = Sx[1][0] = sx;
     Sx[0][2] = Sx[2][0] = Sx[1][1] = sxx;
     Sx[1][2] = Sx[2][1] = sxxx;
     Sx[2][2] = sxxxx;
     Sy[0][0] = sy;
     Sy[1] [0]= sxy;
     Sy[2][0] = sxxy;
  //-------------------------------------------------------------
  // Form matrices for Matrix Class from Java 2D arrays
  //-------------------------------------------------------------
     Matrix MatSx = new Matrix(Sx);
//     Matrix MatSy = new Matrix(Sy);
//-------------------------------------------------------------
//  Form Matrix Sy with Jama set
     Matrix MatSy = new Matrix(3,1);
     MatSy.set(0,0,sy);
     MatSy.set(1,0,sxy);
     MatSy.set(2,0,sxxy);
  //-------------------------------------------------------------
  // Solve matrix equation for B via inverse,  test inverse
  //-------------------------------------------------------------
    Matrix B = MatSx.inverse().times(MatSy);
    Matrix Itest = MatSx.inverse().times(MatSx);
//  Print out matrices using Jama's print 
    System.out.print( "B Matrix via inverse" );
    B.print (16,14);
    System.out.print( "MatSx.inverse().times(MatSx) " );
    Itest.print (16,14);
//-------------------------------------------------------------
  // Solve matrix equation directly for  B
  //-------------------------------------------------------------
    B = MatSx.solve(MatSy);
//  Print out matrices using Jama's print 
    System.out.print( "B Matrix via direct" );
    B.print (16,14);
  //-----------------------------------------------------
  // Extract using Jama get & Print parabola coefficients
  //-----------------------------------------------------
 
    System.out.println( "FitParabola2 Final Results" ); 
    System.out.println( "                                 " );
    System.out.println( "y(x) = b0 + b1 x + b2 x^2           " );
    System.out.println( "                                 " );
    System.out.println("   b0 = "  + B.get(0,0)              );
    System.out.println("   b1 = "  + B.get(1,0)              );
    System.out.println("   b2 = "  + B.get(2,0)              );
    System.out.println( "                                 " );
 
 //test fit 
    for (int i=0; i <= Nd-1; i++)
    {
      s=B.get(0,0)+B.get(1,0)*x[i]+B.get(2,0)*x[i]*x[i];
      System.out.println("   i, xi, yi, yfit = "  
			 +i+ ", "+  x[i] +", " +y[i]+", "  + s );
    }
}   
   
}   // End of file