import java.util.*;
import java.lang.*;
import java.text.*;

public class Fit{
    private double Data[];
    private String Method, Order;
    private int Degree, N, Ncards, nx, Iorder, Imethod;
    private double x[][];
    private double y[];
    private boolean Constant, debug;
    private Hashtable AllowedMethods, AllowedOrders;
    private Util util = new Util();
    private int Nparms;   //Nb. of smoothing parameters 
    private double z[];   //Smoothing parameters
    private double Y[];   //Fit estimates of y
    private double b[];   //X[T]y (b of matrix formulation Ax=b)
    private double rms, Rsquared, MeanSquare, StandardDeviation,
                   correlation;
    public boolean GAUSS_SINGULAR=false;
    public boolean WRONG_Nb_DATA=false;
    public boolean INSUFFICIENT_DATA=false;

    Fit(double Data[], int N, String Method, String Order, int Degree,
        boolean Constant, boolean debug){
        this.N = N;
        this.Data = Data;
        this.Method = Method;
        this.Order = Order;
        this.Degree = Degree;
        this.Constant = Constant;
        this.debug = debug;
        SetupOptions(); function_fit();

    }        
    public void SetupOptions()
    {
        AllowedMethods = new Hashtable();
        AllowedMethods.put("polynome", new Integer(0) );
        AllowedMethods.put("linear",   new Integer(1) );
        AllowedOrders = new Hashtable();
        AllowedOrders.put("xy", new Integer(0) );
        AllowedOrders.put("yx", new Integer(1) );
    }
    public int getMethod()
    {
        Integer n = (Integer)AllowedMethods.get(Method);
        if(n == null){
            System.out.println("Illegal method " +Method +" - EXIT -");
            System.exit(0);
        }
        return(n.intValue());
    }
    public int getOrder()
    {
        Integer n = (Integer)AllowedOrders.get(Order);
        if(n == null){
            System.out.println("Illegal order " +Order +" - EXIT -");
            System.exit(0);
        }
        return(n.intValue());
    }
    public void function_fit()
    {
        int i, j, k, m;
        String StrMethod=null;
        
        Imethod = getMethod();
        Iorder  = getOrder();
        switch(Imethod){
        case 0:       //Polynome
	      nx = 1;
            break;
        case 1:       //Linear
	      nx = Degree;
	      break;
        default:
	      break;
        }
        k = nx + 1;  //Number of values per data card 
        if(N%k != 0){
            WRONG_Nb_DATA = true;
            return;
//	      System.out.print("Illegal number of data - EXIT - ");
//          System.out.print("from function_fit\n");
//	      System.exit(0);
        }
        Ncards = N / k;   //Number of data cards 
        if(Ncards <2){
            INSUFFICIENT_DATA = true;
            return;
        }
//      x = new double[Ncards][]; // These are also valid assignments...
//      for(i=0; i<Ncards; i++)x[i] = new double[nx];
        x = new double[Ncards][nx];
        y = new double[Ncards];
        Y = new double[Ncards];
        for(i=0; i<Ncards; i++){
	      j = i * k;
	      switch(Iorder){
	      case 0:    //input order: x's followed by y
	          y[i] = Data[j+nx];
                for(m=0; m<nx; m++)x[i][m] = Data[j+m];
	          break;
	      case 1:         //input order: y followed by x's 
	          y[i] = Data[j];
                for(m=0; m<nx; m++)x[i][m] = Data[j+1+m];
	          break;
	      default:
	          break;
            }
        }
        z = new double[Degree+1];
        Nparms = Degree + 1;
        if(!Constant)Nparms--;
        switch(Imethod){
        case 0:     //Polynome
	      polynomial_fit();
            if(GAUSS_SINGULAR)return;
	      Yestimates_polynomial();
	      break;
        case 1:     //Linear
	      linear_fit();
            if(GAUSS_SINGULAR)return;
	      Yestimates_linear(); 
	      break;
        default:
	      break;
        }
        if(Degree == 1)CorrelationCoef();
        statistics();
        rms = root_mean_square_error();
        return;          
    }
    public void polynomial_fit()
    {
        int i, j, k, l, n;
        double zero=0.0e0;
        double Xmean, Ymean;

 /*        
        //Subtract from points x it's mean value.
        //This may improve the accuracy of "difficult" cases.
        for(Xmean=zero, i=0; i<Ncards; i++)Xmean += x[i][0];
        Xmean /= Ncards;
        for(Ymean=zero, i=0; i<Ncards; i++)Ymean += y[i];
        Ymean /= Ncards;
        for(i=0; i<Ncards; i++){
            x[i][0] -= Xmean;
 //           y[i]    -= Ymean;
        }
System.out.println(Xmean +" " +Ymean);
*/
        j = 2*Degree + 1;
        double s[] = new double[j];
        b = new double[j];

        s[0] = (double)Ncards;
        for(i=1; i<j; i++){
	      for(s[i]=zero, k=0; k<Ncards; k++){
	          s[i] += Math.pow(x[k][0], (double)i);
	      }
        }
        n = Degree + 1; k = 0; l=0;
        if(!Constant){n = Degree; k = 2; l=1;}
        for(i=0; i<j-k; i++)b[i] = s[k+i];
        double a[][] = new double[n][n];
        for(i=0; i<n; i++){
	      for(j=0; j<=i; j++){
	          a[i][j] = b[i+j];
	          a[j][i] = a[i][j];
	      }
        }
        j = Degree + 1;
        for(s[0]=zero, k=0; k<Ncards; k++)s[0] += y[k];
        for(i=1; i<j; i++){
	      for(s[i]=zero, k=0; k<Ncards; k++){
                s[i] += y[k] * Math.pow(x[k][0], (double)i);
	      }
        }
        for(i=0; i<n; i++)b[i] = s[l+i];
        if(debug){
//            util.CenterString("Linear equation",true);
//            util.output_double_vec(n, b, false, 8);
//            System.out.print("\n");
//            util.printmat(a, s, n, n, false);
        }
        util.AX_Bgauss(n, a, b, z);
        if(util.GAUSS_SINGULAR)GAUSS_SINGULAR = true;
        return;
    }
    public void Yestimates_polynomial()
    {
        int i, k, n;
        double zero=0.0e0;
	
        if(Constant){
	      n = Degree+1;
	      for(i=0; i<Ncards; i++){
	          for(Y[i]=z[0], k=1; k<n; k++){
		        Y[i] += z[k] * Math.pow(x[i][0], (double)k);
	          }
	      }
        }else{
	      n = Degree;
	      for(i=0; i<Ncards; i++){
	          for(Y[i]=zero, k=1; k<=n; k++){
		        Y[i] += z[k-1] * Math.pow(x[i][0], (double)k);
	          }
	      }
        }
        return;
    }
    public void linear_fit()
    {
        int i, j, k, n;
        b = new double[Degree+1];
        double a[][] = new double[Degree+1][Degree+1];
        double zero=0.0e0;

        n = Degree;
        for(i=0; i<n; i++){
	      for(j=0; j<=i; j++){
	          for(a[i][j]=zero, k=0; k<Ncards; k++){
		        a[i][j] += x[k][j] * x[k][i];
	          }
	          a[j][i] = a[i][j];
	      }
            for(b[i]=zero, k=0; k<Ncards; k++){
	          b[i] += x[k][i] * y[k];
            }
        }
        if(Constant){
	      for(i=0; i<n; i++){
	          for(a[n][i]=zero, k=0; k<Ncards; k++){
		        a[n][i] += x[k][i];
	          }
	          a[i][n] = a[n][i];
	      }
	      a[n][n] = (double)Ncards;
	      for(b[n]=zero, k=0; k<Ncards; k++)b[n] += y[k];
	      n++;
        }
        if(debug){
//            util.CenterString("Linear equation",true);
//            util.output_double_vec(n, b, false, 8);
//            System.out.print("\n");
//            util.printmat(a, b, n, n, false);
        }
        util.AX_Bgauss(n, a, b, z);
        if(util.GAUSS_SINGULAR)GAUSS_SINGULAR = true;
        return;
    }
    public void Yestimates_linear()
    {
        int i, k;
        double zero=0.0e0;

        for(i=0; i<Ncards; i++){
	      for(Y[i]=zero, k=0; k<Degree; k++){
	          Y[i] += z[k] * x[i][k];
	      }
	      if(Constant)Y[i] += z[Degree];
        }
        return;
    }
    public void statistics()
    {
        int i, j;
        double zero=0.0e0, one=1.0e0, Bterm, yy, SS, denom,
               yAv, yAvSquared, delta=1.0e-13;

        for(Bterm=zero, i=0; i<Nparms; i++)Bterm += b[i]*z[i];
        for(yy=zero, i=0; i<Ncards; i++)yy += y[i]*y[i];
        SS = yy - Bterm;
        for(yAv=zero, i=0; i<Ncards; i++)yAv += y[i];
        yAv /= Ncards;
        denom = yy;
        if(Constant)denom -= Ncards * yAv * yAv;
        Rsquared = one - SS / denom;
        MeanSquare = SS / (double)(Ncards - Nparms);
        if(Math.abs(MeanSquare) < delta)MeanSquare = zero;
        StandardDeviation = Math.sqrt(MeanSquare);
        return;
    }
    public void CorrelationCoef()
    {
        int i, k;
        double zero=0.0e0, one=1.0e0, Yav, Sy2, Sy, Sxy,
               prod, denom;

        double Xav[] = new double[Degree];
        double Sx2[] = new double[Degree];
        double  Sx[] = new double[Degree];

        //Averages
        for(i=0; i<Degree; i++){
            for(Xav[i]=zero, k=0; k<Ncards; k++){
                Xav[i] += x[k][i];
            }
            Xav[i] /= (double)Ncards;
        }
        for(Yav=zero, i=0; i<Ncards; i++){
            Yav += y[i];
        }
        Yav /= (double)Ncards;
        //S2 terms
        for(i=0; i<Degree; i++){
            for(Sx2[i]=zero, k=0; k<Ncards; k++){
                Sx2[i] += (x[k][i] - Xav[i]) * (x[k][i] - Xav[i]);
            }
            Sx2[i] /= (double)Ncards;
        }
        for(Sy2=zero, i=0; i<Ncards; i++){
            Sy2 += (y[i] - Yav) * (y[i] - Yav);
        }
        Sy2 /= (double)Ncards;
        //Sxy
        for(Sxy=zero, i=0; i<Ncards; i++){
            for(prod=one, k=0; k<Degree; k++){
                prod *= (x[i][k] - Xav[k]);
            }
            prod *= (y[i] - Yav);
            Sxy += prod;
        }
        Sxy /= (double)Ncards;
        //S terms
        for(i=0; i<Degree; i++){
            Sx[i] = Math.sqrt(Sx2[i]);
        }
        Sy = Math.sqrt(Sy2);
        //Correlation coefficient
        for(denom=one, i=0; i<Degree; i++){
            denom *= Sx[i];
        }
        denom *= Sy;
        correlation = Sxy / denom;
    }
    public double root_mean_square_error()
    {
        int i;
        double zero=0.0e0, rms;

        for(rms=zero, i=0; i<Ncards; i++){
            rms += (Y[i] - y[i]) * (Y[i] - y[i]);
        }
        rms /= (double)Ncards;
        return( Math.sqrt(rms) );
    }
    public void putNparms(int n){
        Nparms = n;
        return;
    }
    //Return the (linear or polynomial) degree
    public int getDegree(){
        return(Degree);
    }
    //Return the "compute constant" boolean variable
    public boolean isConstant(){
        return(Constant);
    }
    //Return Nb. of smoothing parameters
    public int getNparms(){
        return(Nparms);
    }
    //Return smoothing parameters
    public double[] getParameters(){
        return(z);
    }
    //Return Nb. of input cards
    public int getN(){
        return(Ncards);
    }
    //Return Nb. of x's/card
    public int getNx(){
        return(nx);
    }
    //Return input x array
    public double[][] getX(){
        return(x);
    }
    //Return input y vector
    public double[] getY(){
        return(y);
    }
    //Return output (fitted estimates) y vector
    public double[] getYY(){
        return(Y);
    }
    //Return rms statistic
    public double getRMS(){
        return(rms);
    }
    //Return correlation statistic
    public double getCorrelation(){
        return(correlation);
    }
    //Return R-squared statistic
    public double getRsquared(){
        return(Rsquared);
    }
    //Return residual mean square statistic
    public double getMeanSquare(){
        return(MeanSquare);
    }
    //Return residual standard deviation statistic
    public double getStandardDeviation(){
        return(StandardDeviation);
    }
}