// least_square_fit_4d.java  tailorable code, provide your input and setup 

// The purpose of this package is to provide a reliable and convenient
// means for fitting existing data by a few coefficients. The companion
// package check_fit provides the means to use the coefficients for
// interpolation and limited extrapolation.
// 
// This package implements the least square fit. 
// 
// The problem is stated as follows :
//   Given measured data for values of Y based on values of X1,X2 and X3. e.g.
// 
//    Y_actual         X1      X2     X3     X4
//    --------       -----   -----  -----  -----
//     32.5           1.0     2.5    3.7    4.2
//      7.2           2.0     2.5    3.6    2.5
//      6.9           3.0     2.7    3.5    1.9
//     22.4           2.2     2.1    3.1    4.7
//     10.4           1.5     2.0    2.6    0.9
//     11.3           1.6     2.0    3.1    5.1
// 
//  Find a, b, ... such that   Y_approximate =  a + b * X1 + c * X2 + d * X3 +
//  e * X4 + f * X1*X2 + g *X1*X3 ... k * X1*X2*X3*X4 ... z * X4*X4*X4*X4 
//  Find the 69 coefficients, including the constant
//  such that the sum of (Y_actual - Y_approximate) squared is minimized.
//  This is minimizing the root mean squire error.
// 
// The method for determining the coefficients a, b and c follows directly
// form the problem definition and mathematical analysis. (See more below)
// 
// Y is called the dependent variable and X1 .. Xn the independent variables.
// The procedures below implements a few special cases and the general case.
//    The number of independent variables can vary.
//    The approximation equation may use powers of the independent variables
//    The user may create additional independent variables e.g. X2 = SIN(X1)
//    with the restriction that the independent variables are linearly
//    independent.  e.g.  Xi not equal  p Xj + q  for all i,j,p,q
// 
// 
// 
// The mathematical derivation of the least square fit is as follows :
// 
// Given data for the independent variable Y in terms of the dependent
// variables S,T,U and V  consider that there exists a function F
// such that     Y = F(S,T,U,V)
// The problem is to find coefficients a,b,c and d such that
//            Y_approximate = a * S + b * T + c * U + d * V
// and such that the sum of ( Y - Y_approximate ) squared is minimized.
// 
// Note: a, b, c, d are scalars. S, T, U, V, Y, Y_approximate are vectors.
// 
// To find the minimum of  SUM( Y - Y_approximate ) ** 2
// the derivatives must be taken with respect to a,b,c and d and
// all must equal zero simultaneously. The steps follow :
// 
//  SUM( Y - Y_approximate ) ** 2 = SUM( Y - a*S - b*T - c*U - d*V ) ** 2
// 
// d/da =  -2 * S * SUM( Y - A*S - B*T - C*U - D*V )
// d/db =  -2 * T * SUM( Y - A*S - B*T - C*U - D*V )
// d/dc =  -2 * U * SUM( Y - A*S - B*T - C*U - D*V )
// d/dd =  -2 * V * SUM( Y - A*S - B*T - C*U - D*V )
// 
// Setting each of the above equal to zero and putting constant term on left
//    the -2 is factored out,
//    the independent variable is moved inside the summation
// 
//  SUM( a*S*S + b*S*T + c*S*U + d*S*V = S*Y )
//  SUM( a*T*S + b*T*T + c*T*U + d*T*V = T*Y )
//  SUM( a*U*S + b*U*T + c*U*U + d*U*V = U*Y )
//  SUM( a*V*S + b*V*T + c*V*U + d*V*V = V*Y )
// 
// Distributing the SUM inside yields
// 
//  a * SUM(S*S) + b * SUM(S*T) + c * SUM(S*U) + d * SUM(S*V) = SUM(S*Y)
//  a * SUM(T*S) + b * SUM(T*T) + c * SUM(T*U) + d * SUM(T*V) = SUM(T*Y)
//  a * SUM(U*S) + b * SUM(U*T) + c * SUM(U*U) + d * SUM(U*V) = SUM(U*Y)
//  a * SUM(V*S) + b * SUM(V*T) + c * SUM(V*U) + d * SUM(V*V) = SUM(V*Y)
// 
// To find the coefficients a,b,c and d solve the linear system of equations
// 
//    | SUM(S*S)  SUM(S*T)  SUM(S*U)  SUM(S*V) |   | a |   | SUM(S*Y) |
//    | SUM(T*S)  SUM(T*T)  SUM(T*U)  SUM(T*V) | x | b | = | SUM(T*Y) |
//    | SUM(U*S)  SUM(U*T)  SUM(U*U)  SUM(U*V) |   | c |   | SUM(U*Y) |
//    | SUM(V*S)  SUM(V*T)  SUM(V*U)  SUM(V*V) |   | d |   | SUM(V*Y) |
// 
// Some observations :
//     S,T,U and V must be linearly independent.
//     There must be more data sets (Y, S, T, U, V) than variables.
//     The analysis did not depend on the number of independent variables
//     A polynomial fit results from the substitutions S=1, T=X, U=X**2, V=X**3
//     The general case for any order polynomial follows, fit_pn.
//     Any substitution such as three variables to various powers may be used.

public class least_square_fit_4d
{
  double rms_err, avg_err, max_err;
  int idata = 0; // reset after each use of the data set
  int Apwr[] = new int[200]; // powers of each variable in a term
  int Bpwr[] = new int[200];
  int Cpwr[] = new int[200];
  int Dpwr[] = new int[200];
  int debug = 0;
  

  least_square_fit_4d()
  {
    int n;
    int mm=4; // 4 dimensions
    int nnn[] = new int[1]; // number of terms 
  
    System.out.println("least_square_fit_4d.java");
  
    // sample polynomial least square fit, 3th power 
    {
      n=3; // need constant term and powers 1,2,3
      System.out.println("fit u "+(n)+" degree polynomial in 4 dimensions");
      gen_4d_powers(n, mm, nnn, Apwr, Bpwr, Cpwr, Dpwr);
      double A[][] = new double[nnn[0]][nnn[0]];
      double C[] = new double[nnn[0]];
      double Y[] = new double[nnn[0]];
      fit_4d(n, mm, nnn, A, Y, C);
      check_4d(n, mm, nnn, C);
      System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                         ", max_err="+max_err);
      System.out.println(" ");
    }
    {
      n=4; // need constant term and powers 1,2,3,4 
      System.out.println("fit u "+(n)+" degree polynomial in 4 dimensions");
      gen_4d_powers(n, mm, nnn, Apwr, Bpwr, Cpwr, Dpwr);
      double A[][] = new double[nnn[0]][nnn[0]];
      double C[] = new double[nnn[0]];
      double Y[] = new double[nnn[0]];
      fit_4d(n, mm, nnn, A, Y, C);
      check_4d(n, mm, nnn, C);
      System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                         ", max_err="+max_err);
      System.out.println(" ");
    }
  } // end least_square_fit_4d

    int data_set4d(double abcdu[], int ijkm[])
  {
    int i=ijkm[0];
    int j=ijkm[1];
    int k=ijkm[2];
    int m=ijkm[3]; // would have used l, ell, yet too close to i and 1
    int n = 4;
    double u, a, b, c, d;

    i++;
    if(i>n) {i=0; j++;}
    if(j>n) {j=0; k++;}
    if(k>n) {k=0; m++;}
    if(m>n) {i=0; j=0; k=0; m=0; return 0;}
    a = (double)i;
    b = (double)j;
    c = (double)k;
    d = (double)m;
    u = 1.0 + 2.0*a + 3.0*b + 4.0*c + 5.0*d +
        6.0*a*a + 7.0*a*b + 8.0*a*c + 9.0*a*d +
        10.0*b*b + 11.0*b*c + 12.0*b*d + 13.0*c*c +
        14.0*c*d + 15.0*d*d + 16.0*a*a*a + 17.0*a*a*b +
        18.0*a*a*c + 19.0*a*a*d + 20.0*a*b*b + 21.0*a*b*c +
        22.0*a*b*d + 23.0*a*c*c + 24.0*a*c*d + 25.0*a*d*d +
        26.0*b*b*b + 27.0*b*b*c + 28.0*b*b*d + 29.0*b*c*c +
        30.0*b*c*d + 31.0*c*c*c + 32.0*c*c*d + 33.0*c*d*d +
        34.0*d*d*d + 35.0*a*a*a*a + 36.0*a*a*a*b +
        37.0*a*a*a*c + 38.0*a*a*a*d + 39.0*a*a*b*b +
        40.0*a*a*b*c + 41.0*a*a*b*d + 42.0*a*a*c*c +
        43.0*a*a*c*d + 44.0*a*a*d*d + 45.0*a*b*b*b +
        46.0*a*b*b*c + 47.0*a*b*b*d + 48.0*a*b*c*c +
	49.0*a*b*c*d + 50.0*a*a*d*d + 51.0*a*c*c*c +
        52.0*a*c*c*d + 53.0*a*c*d*d + 54.0*a*d*d*d +
        55.0*b*b*b*b + 56.0*b*b*b*c + 57.0*b*b*b*d +
        58.0*b*b*c*c + 59.0*b*b*c*d + 60.0*b*b*d*d +
        61.0*b*c*c*c + 62.0*b*c*c*d + 63.0*b*c*d*d +
        64.0*b*d*d*d + 65.0*c*c*c*c + 66.0*c*c*c*d +
        67.0*c*c*d*d + 68.0*c*d*d*d + 69.0*d*d*d*d;

    abcdu[0]=a;
    abcdu[1]=b;
    abcdu[2]=c;
    abcdu[3]=d;
    abcdu[4]=u;
    ijkm[0]=i;
    ijkm[1]=j;
    ijkm[2]=k;
    ijkm[3]=m;
    return 1;
  } // end data_set4d

  void fit_4d(int n, int mm, int nnn[], double A[][], double Y[], double C[])
  {
    int i, j, k, nn;
    double Av[] = new double[10]; // at least n
    double Bv[] = new double[10]; // powers of variables
    double Cv[] = new double[10];
    double Dv[] = new double[10];
    double abcdu[] = new double[5];
    int ijkm[] = {0,0,0,0};
    double term_i, term_j;
    nn = nnn[0];
    for(i=0; i<nn; i++)
    {
      for(j=0; j<nn; j++)
      {
        A[i][j] = 0.0;
      }
      Y[i] = 0.0;
    }
    while(data_set4d(abcdu, ijkm)==1)
    {
      Av[0] = 1.0;
      Bv[0] = 1.0;
      Cv[0] = 1.0;
      Dv[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*abcdu[0]; // a
        Bv[i] = Bv[i-1]*abcdu[1]; // b
	Cv[i] = Cv[i-1]*abcdu[2]; // c
	Dv[i] = Dv[i-1]*abcdu[3]; // d
      }
      for(i=0; i<nn; i++)
      {
	term_i = Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]];
	for(j=0; j<nn; j++)
	{
	  term_j = Av[Apwr[j]] * Bv[Bpwr[j]] * Cv[Cpwr[j]] * Dv[Dpwr[j]];
	  A[i][j] = A[i][j] + term_i * term_j;
	}
	Y[i] = Y[i] + abcdu[4] * term_i;
      }
    }
    simeq(A, Y, C);
    for(i=0; i<nn; i++) System.out.println("C["+i+"]="+C[i]);
  } // end fit_4d 

  void check_4d(int n, int mm, int nnn[], double X[])
  {
    double a1, b1, c1, d1, u1, ua, diff;
    double sumsq = 0.0;
    double sum = 0.0;
    double maxe = 0.0;
    double amin=0.0, amax=0.0, bmin=0.0, bmax=0.0, cmin=0.0;
    double cmax=0.0, dmin=0.0, dmax=0.0, umin=0.0, umax=0.0;
    double Av[] = new double[10]; // at least n
    double Bv[] = new double[10]; // powers of variables
    double Cv[] = new double[10];
    double Dv[] = new double[10];
    int i, k;
    int nn = nnn[0];
    double abcdu[] = new double[5];
    int ijkm[] = {0,0,0,0};
  
    k = 0;
    while(data_set4d(abcdu, ijkm)==1)
    {
      a1=abcdu[0];
      b1=abcdu[1];
      c1=abcdu[2];
      d1=abcdu[3];
      u1=abcdu[4];
      Av[0] = 1.0;
      Bv[0] = 1.0;
      Cv[0] = 1.0;
      Dv[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*a1;
	Bv[i] = Bv[i-1]*b1;
	Cv[i] = Cv[i-1]*c1;
	Dv[i] = Dv[i-1]*d1;
      }

      if(k==0)
      {
	amin=a1; amax=a1;
	bmin=b1; bmax=b1;
	cmin=c1; cmax=c1;
	dmin=c1; dmax=c1;
	umin=u1; umax=u1;
      }
      if(a1>amax) amax=a1;
      if(a1<amin) amin=a1;
      if(b1>bmax) bmax=b1;
      if(b1<bmin) bmin=b1;
      if(c1>cmax) cmax=c1;
      if(c1<cmin) cmin=c1;
      if(d1>dmax) dmax=d1;
      if(d1<dmin) dmin=d1;
      if(u1>umax) umax=u1;
      if(u1<umin) umin=u1;

      ua = 0.0;
      for(i=0; i<nn; i++)
      {
	ua = ua + Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]] * X[i];
      }
      diff = Math.abs(u1-ua);
      if(diff>maxe) maxe=diff;
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
      k++;
    }
    System.out.println("check_4d k="+k+", amin="+amin+", amax="+amax+
                       ", bmin="+bmin+", bmax="+bmax);
    System.out.println("         cmin="+cmin+", cmax="+cmax+", dmin="+dmin+
                       ", dmax="+dmax+", umin="+umin+", umax="+umax);
    max_err = maxe;
    avg_err = sum/(double)k;
    rms_err = Math.sqrt(sumsq/(double)k);
  } // end check_4d 


  //  n=highest sum of powers, maximum 4 here.
  //  mm=number of independent variables, must be 4 here.
  //  nnn= number of coefficients, returned
  void gen_4d_powers(int n, int mm, int nnn[], // only one value returned
		     int a[], int b[], int c[], int d[])
  {
    int pwrsix[]={0,1,5,15,34,69}; // start of each order
    int pwrs[][]={{0,0,0,0},{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1},
		  {2,0,0,0},{1,1,0,0},{1,0,1,0},{1,0,0,1},{0,2,0,0},
		  {0,1,1,0},{0,1,0,1},{0,0,2,0},{0,0,1,1},{0,0,0,2},
		  {3,0,0,0},{2,1,0,0},{2,0,1,0},{2,0,0,1},{1,2,0,0},
		  {1,1,1,0},{1,1,0,1},{1,0,2,0},{1,0,1,1},{1,0,0,2},
		  {0,3,0,0},{0,2,1,0},{0,2,0,1},{0,1,2,0},{0,1,1,1},
		  {0,0,3,0},{0,0,2,1},{0,0,1,2},{0,0,0,3},{4,0,0,0},
		  {3,1,0,0},{3,0,1,0},{3,0,0,1},{2,2,0,0},{2,1,1,0},
		  {2,1,0,1},{2,0,2,0},{2,0,1,1},{2,0,0,2},{1,3,0,0},
		  {1,2,1,0},{1,2,0,1},{1,1,2,0},{1,1,1,1},{1,1,0,2},
		  {1,0,3,0},{1,0,2,1},{1,0,1,2},{1,0,0,3},{0,4,0,0},
		  {0,3,1,0},{0,3,0,1},{0,2,2,0},{0,2,1,1},{0,2,0,2},
		  {0,1,3,0},{0,1,2,1},{0,1,1,2},{0,1,0,3},{0,0,4,0},
		  {0,0,3,1},{0,0,2,2},{0,0,1,3},{0,0,0,4}};
    if(mm != 4) System.out.println(
		"ERROR, this only good for mm=4 independent variables");
    System.out.println("terms used to find fit ");
    for(int i=0; i<=n; i++) // power requested
    {
      for(int ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
      {
	a[ii]=pwrs[ii][0];
	b[ii]=pwrs[ii][1];
	c[ii]=pwrs[ii][2];
	d[ii]=pwrs[ii][3];
	System.out.println(ii+"   a^"+a[ii]+" * b^"+b[ii]+" * c^"+c[ii]+
			   " * d^"+d[ii]);
      }
      System.out.println(" ");
    }
    nnn[0] = pwrsix[n+1];
  } /* end gen_4d_powers */
    
  void simeq(final double A[][], final double Y[], double X[])
  {
    // solve real linear equations for X where Y = A * X
    // method: Gauss-Jordan elimination using maximum pivot
    // usage:  simeq(A,Y,X);
    //    Translated to java by : Jon Squire , 26 March 2003
    //    First written by Jon Squire December 1959 for IBM 650, translated to
    //    other languages  e.g. Fortran converted to Ada converted to C
    //    then converted to java
    int n=A.length;
    int m=n+1;
    double B[][]=new double[n][m];  // working matrix
    int row[]=new int[n];           // row interchange indicies
    int hold , I_pivot;             // pivot indicies
    double pivot;                   // pivot element value
    double abs_pivot;
    if(A[0].length!=n || Y.length!=n || X.length!=n)
    {
      System.out.println("Error in Matrix.solve, inconsistent array sizes.");
    }
    // build working data structure
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        B[i][j] = A[i][j];
      }
      B[i][n] = Y[i];
    }
    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k;
    }
    //  begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = Math.abs(pivot);
      I_pivot = k;
      for(int i=k; i<n; i++)
      {
        if(Math.abs(B[row[i]][k]) > abs_pivot)
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = Math.abs(pivot);
        }
      }
      // have pivot, interchange row indicies
      hold = row[k];
      row[k] = row[I_pivot];
      row[I_pivot] = hold;
      // check for near singular
      if(abs_pivot < 1.0E-10)
      {
        for(int j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = 0.0;
        }
        System.out.println("redundant row (singular) "+row[k]);
      } // singular, delete row
      else
      {
        // reduce about pivot
        for(int j=k+1; j<n+1; j++)
        {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if( i != k)
          {
            for(int j=k+1; j<n+1; j++)
            {
              B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j];
            }
          }
        }
      }
      //  finished inner reduction
    }
    //  end main reduction loop
    //  build  X  for return, unscrambling rows
    for(int i=0; i<n; i++)
    {
      X[i] = B[row[i]][n];
    }
  } // end simeq 

  public static void main (String[] args)
  {
    new least_square_fit_4d();
  } // end main

} // end class least_square_fit_4d