// fit_data_file.java  least square fit of data file 
//                     java -cp . fit_data_file  vars  order  file_name
//                     vars  1,2,3 or 4  number of independent variables
//                     order < 100 for one variable
//                     order <  33 for two variables
//                     order <  10 for three variables
//                     order <   5 for four variables

import java.io.*;

public class fit_data_file
{
  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];
  boolean first = true; // uses by data_set, fit, then check 
  String data_file_name;
  FileReader handle;
  BufferedReader file_in;
  String input_line = new String("");
  int found=0;
  int line=0; // input line count
  double tol = 1.0E-4; // minimum value of coefficient

  fit_data_file(int vars, int order, String file_name)
  {
    data_file_name=file_name;
    int n=order; // max depends on number of variables
    int mm=vars; // max 4 dimensions
    int nnn[] = new int[1]; // number of terms computed
  
    System.out.println("fit_data_file.java");
    System.out.println("fit to "+n+" degree polynomial in "+mm+" variables");
    if(mm==1)
    {
      gen_1d_powers(n, mm, nnn, Apwr);
      double A[][] = new double[nnn[0]][nnn[0]];
      double C[] = new double[nnn[0]];
      double Y[] = new double[nnn[0]];
      fit_1d(n, mm, nnn, A, Y, C);
      check_1d(n, mm, nnn, C);
    }
    else if(mm==2)
    {
      gen_2d_powers(n, mm, nnn, Apwr, Bpwr);
      double A[][] = new double[nnn[0]][nnn[0]];
      double C[] = new double[nnn[0]];
      double Y[] = new double[nnn[0]];
      fit_2d(n, mm, nnn, A, Y, C);
      check_2d(n, mm, nnn, C);
    }
    else if(mm==3)
    {
      gen_3d_powers(n, mm, nnn, Apwr, Bpwr, Cpwr);
      double A[][] = new double[nnn[0]][nnn[0]];
      double C[] = new double[nnn[0]];
      double Y[] = new double[nnn[0]];
      fit_3d(n, mm, nnn, A, Y, C);
      check_3d(n, mm, nnn, C);
    }
    else if(mm==4)
    {
      if(n<1 || n>4)
      {
        System.out.println("order must be 1,2,3 or 4  to fit 4 variables");
	System.exit(1);
      }
      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);
    }
    else
    {
      System.out.println("number of variables  must be 1,2,3 or 4  to fit");
      System.exit(1);
    }
  } // end fit_data_file

  void fit_1d(int n, int mm, int nnn[], double A[][], double Y[], double C[])
  {
    int i, j, k, nn;
    double Av[] = new double[100]; // at least n
    double term_i, term_j;
    double au[] = new double[2];
    boolean debug = false;

    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_set1d(au)==1)
    {
      if(debug) System.out.println("have data a="+au[0]);
      Av[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*au[0]; // a
      }
      for(i=0; i<nn; i++)
      {
	term_i = Av[Apwr[i]];
	for(j=0; j<nn; j++)
	{
	  term_j = Av[Apwr[j]];
	  A[i][j] = A[i][j] + term_i * term_j;
	}
	Y[i] = Y[i] + au[1] * term_i;
      }
    }
    if(debug)
    {
      for(i=0; i<nn; i++)
      {
	for(j=0; j<nn; j++)
	{
	  System.out.println("A["+i+"]["+j+"]="+A[i][j]);
	}
        System.out.println("Y["+i+"]="+Y[i]);
      }
    } // end debug
    simeq(A, Y, C);
    System.out.println("Polynomial Coefficients");
    for(i=0; i<nn; i++)
    {
      if(Math.abs(C[i])>tol)
      {
	System.out.println(C[i]+" * a^"+Apwr[i]);
      }
    }
  } // end fit_1d 

  void fit_2d(int n, int mm, int nnn[], double A[][], double Y[], double C[])
  {
    int i, j, k, nn;
    double Av[] = new double[33]; // at least n
    double Bv[] = new double[33]; // powers of variables
    double term_i, term_j;
    double abu[] = new double[3];
    boolean debug = false;

    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_set2d(abu)==1)
    {
      if(debug) System.out.println("have data a="+abu[0]);
      Av[0] = 1.0;
      Bv[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*abu[0]; // a
        Bv[i] = Bv[i-1]*abu[1]; // b
      }
      for(i=0; i<nn; i++)
      {
	term_i = Av[Apwr[i]] * Bv[Bpwr[i]];
	for(j=0; j<nn; j++)
	{
	  term_j = Av[Apwr[j]] * Bv[Bpwr[j]];
	  A[i][j] = A[i][j] + term_i * term_j;
	}
	Y[i] = Y[i] + abu[2] * term_i;
      }
    }
    simeq(A, Y, C);
    System.out.println("Polynomial Coefficients");
    for(i=0; i<nn; i++)
    {
      if(Math.abs(C[i])>tol)
      {
	System.out.println(C[i]+" * a^"+Apwr[i]+" * b^"+Bpwr[i]);
      }
    }
  } // end fit_2d 

  void fit_3d(int n, int mm, int nnn[], double A[][], double Y[], double C[])
  {
    int i, j, k, nn;
    double Av[] = new double[20]; // at least n
    double Bv[] = new double[20]; // powers of variables
    double Cv[] = new double[20];
    double term_i, term_j;
    double abcu[] = new double[4];
    boolean debug = false;

    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_set3d(abcu)==1)
    {
      if(debug) System.out.println("have data a="+abcu[0]);
      Av[0] = 1.0;
      Bv[0] = 1.0;
      Cv[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*abcu[0]; // a
        Bv[i] = Bv[i-1]*abcu[1]; // b
	Cv[i] = Cv[i-1]*abcu[2]; // c
      }
      for(i=0; i<nn; i++)
      {
	term_i = Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]];
	for(j=0; j<nn; j++)
	{
	  term_j = Av[Apwr[j]] * Bv[Bpwr[j]] * Cv[Cpwr[j]];
	  A[i][j] = A[i][j] + term_i * term_j;
	}
	Y[i] = Y[i] + abcu[3] * term_i;
      }
    }
    simeq(A, Y, C);
    System.out.println("Polynomial Coefficients");
    for(i=0; i<nn; i++)
    {
      if(Math.abs(C[i])>tol)
      {
	System.out.println(C[i]+" * a^"+Apwr[i]+" * b^"+Bpwr[i]+
                                " * c^"+Cpwr[i]);
      }
    }
  } // end fit_3d 

  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 term_i, term_j;
    double abcdu[] = new double[5];
    boolean debug = false;

    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)==1)
    {
      if(debug) System.out.println("have data a="+abcdu[0]);
      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);
    System.out.println("Polynomial Coefficients");
    for(i=0; i<nn; i++)
    {
      if(Math.abs(C[i])>tol)
      {
	System.out.println(C[i]+" * a^"+Apwr[i]+" * b^"+Bpwr[i]+
                                " * c^"+Cpwr[i]+" * d^"+Dpwr[i]);
      }
    }
  } // end fit_4d 


  void gen_1d_powers(int n, int mm, int nnn[], // only one value returned
		     int a[])
  {
    int i; /* need more, or generalize, for more variables */

    if(mm != 1) System.out.println(
                "ERROR, this only good for 1 independent variable");
    for(i=0; i<=n; i++) a[i]=i;
    nnn[0] = n+1;
    System.out.println("trying "+nnn[0]+" terms");
  } // end gen_1d_powers

  void gen_2d_powers(int n, int mm, int nnn[], // only one value returned
		     int a[], int b[])
  {
    int i,j; /* need more, or generalize, for more variables */
    int ii = 0; /* pointer to next available a[ii], b[ii] */  
    int nn=1;  /* power being generated */
    boolean debug = true;

    if(mm != 2) System.out.println(
                "ERROR, this only good for 2 independent variables");
    a[0]=0;
    b[0]=0;
    if(debug) System.out.println("terms used to find fit");
    if(debug) System.out.println("a^0 b^0 \n");
  
    i=nn;
    j=0;
  
    while(nn<=n) /* n is highest sum of powers */
    {
      ii++;
      a[ii] = i;
      b[ii] = j;
      System.out.println(ii+"   a^"+a[ii]+" * b^"+b[ii]);
      if(i==0 && j==nn) /* increment nn, set i,j */
      {
        nn++;
        i = nn;
        j = 0;
        if(debug) System.out.println("");
      }
      else
      {
        i--;
        j++;
      }
    }
    nnn[0] = ii+1;
    System.out.println("trying "+nnn[0]+" terms");
  } // end gen_2d_powers

  void gen_3d_powers(int n, int mm, int nnn[], // only one value returned
		     int a[], int b[], int c[])
  {
    int i,j,k; /* need more, or generalize, for more variables */
    int ii = 0; /* pointer to next available a[ii], b[ii], c[ii] */  
    int nn=1;  /* power being generated */
    boolean debug = true;

    if(mm != 3) System.out.println(
                "ERROR, this only good for 3 independent variables");
    a[0]=0;
    b[0]=0;
    c[0]=0;
    if(debug) System.out.println("terms used to find fit");
    if(debug) System.out.println("a^0 b^0 c^0 \n");
  
    i=nn;
    j=0;
    k=0;
  
    while(nn<=n) /* n is highest sum of powers */
    {
      ii++;
      a[ii] = i;
      b[ii] = j;
      c[ii] = k;
      System.out.println(ii+"   a^"+a[ii]+" * b^"+b[ii]+" * c^"+c[ii]);
      if(i==0 && j==0 && k==nn) /* increment nn, set i,j,k */
      {
        nn++;
        i = nn;
        j = 0;
        k = 0;
        if(debug) System.out.println("");
      }
      else if(i==nn)
      {
        i--;
        j=1; /* k should be zero */
      }
      else if(j==nn)
      {
        j--; /* i should be zero */
        k=1;
      }
      else if(j>0)
      {
        j--;
        k++;
      }
      else if(j==0)
      {
        i--;
        j=nn-i;
        k=0;
      }
    }
    nnn[0] = ii+1;
    System.out.println("trying "+nnn[0]+" terms");
  } // end gen_3d_powers

  void gen_4d_powers(int n, int mm, int nnn[], // only one value returned
		     int a[], int b[], int c[], int d[])
  {
    //  n=highest sum of powers, maximum 4 here.
    //  mm=number of independent variables, maximum 4 for 4 variables.
    //  nnn= number of coefficients, returned
    int pwrs44ix[]={0,1,5,15,34,66}; // start of each order
    int pwrs4[][]={{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,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=pwrs44ix[i]; ii<pwrs44ix[i+1]; ii++)
      {
	a[ii]=pwrs4[ii][0];
	b[ii]=pwrs4[ii][1];
	c[ii]=pwrs4[ii][2];
	d[ii]=pwrs4[ii][3];
	System.out.println(ii+"   a^"+a[ii]+" * b^"+b[ii]+" * c^"+c[ii]+
			      " * d^"+d[ii]);
      }
      System.out.println(" ");
    }
    nnn[0] = pwrs44ix[n+1];
    System.out.println("trying "+nnn[0]+" terms");
  } /* end gen_4d_powers */


  int data_set1d(double au[])
  {
    boolean debug = false;
    double a=0.0, u=0.0;
    int index, last, len; // for parsing a, u
    try
    {
      if(first)
      {
        handle  = new FileReader(data_file_name);
        file_in = new BufferedReader(handle);
        input_line = file_in.readLine();
        line = 0;
        found = 1;
        first = false;
      }
      if(input_line != null)
      {
	if(debug) System.out.println(input_line);
        // process input
	index = 0; // start of u
	if(debug) System.out.println("start u index="+index);
	last = input_line.indexOf(' ',index+1); // end of u
	if(debug) System.out.println("last of u ="+last);
	u = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("u="+u);
	index = input_line.indexOf(' ',last); // start of a
	if(debug) System.out.println("start a index="+index);
	last = input_line.length();
	if(debug) System.out.println("last of a ="+last);
	a = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("a="+a);
        input_line = file_in.readLine();
        line++;
      }
      else
      {
        file_in.close();
        System.out.println(line+" lines read");
        first = true;
        found = 0;
        return found;
      }
    }
    catch(Exception exception) // catches all exceptions
    {
      // do nothing, keep going
      System.out.println("file read exception");
    }
    au[0] = a;
    au[1] = u;
    return found;
  } // end data_set1d

  int data_set2d(double abu[])
  {
    boolean debug = false;
    double a=0.0, b=0.0, u=0.0;
    int index, last, len; // for parsing a, b, u
    try
    {
      if(first)
      {
        handle  = new FileReader(data_file_name);
        file_in = new BufferedReader(handle);
        input_line = file_in.readLine();
        line = 0;
        found = 1;
        first = false;
      }
      if(input_line != null)
      {
	if(debug) System.out.println(input_line);
        // process input
	index = 0; // start of u
	if(debug) System.out.println("start u index="+index);
	last = input_line.indexOf(' ',index+1); // end of u
	if(debug) System.out.println("last of u ="+last);
	u = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("u="+u);
	index = input_line.indexOf(' ',last); // start of a
	if(debug) System.out.println("start a index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of a ="+last);
	a = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("a="+a);
	index = input_line.indexOf(' ',last); // start of b
	if(debug) System.out.println("start b index="+index);
	last = input_line.length();
	if(debug) System.out.println("last of b ="+last);
	b = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("b="+b);
        input_line = file_in.readLine();
        line++;
      }
      else
      {
        file_in.close();
        System.out.println(line+" lines read");
        first = true;
        found = 0;
        return found;
      }
    }
    catch(Exception exception) // catches all exceptions
    {
      // do nothing, keep going
      System.out.println("file read exception");
    }
    abu[0] = a;
    abu[1] = b;
    abu[2] = u;
    return found;
  } // end data_set2d

  int data_set3d(double abcu[])
  {
    boolean debug = false;
    double a=0.0, b=0.0, c=0.0, u=0.0;
    int index, last, len; // for parsing a, b, c, u
    try
    {
      if(first)
      {
        handle  = new FileReader(data_file_name);
        file_in = new BufferedReader(handle);
        input_line = file_in.readLine();
        line = 0;
        found = 1;
        first = false;
      }
      if(input_line != null)
      {
	if(debug) System.out.println(input_line);
        // process input
	index = 0; // start of u
	if(debug) System.out.println("start u index="+index);
	last = input_line.indexOf(' ',index+1); // end of u
	if(debug) System.out.println("last of u ="+last);
	u = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("u="+u);
	index = input_line.indexOf(' ',last); // start of a
	if(debug) System.out.println("start a index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of a ="+last);
	a = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("a="+a);
	index = input_line.indexOf(' ',last); // start of b
	if(debug) System.out.println("start b index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of b ="+last);
	b = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("b="+b);
	index = input_line.indexOf(' ',last); // start of c
	if(debug) System.out.println("start c index="+index);
	last = input_line.length();
	if(debug) System.out.println("last of c ="+last);
	c = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("c="+c);
        input_line = file_in.readLine();
        line++;
      }
      else
      {
        file_in.close();
        System.out.println(line+" lines read");
        first = true;
        found = 0;
        return found;
      }
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("file read exception");
      try{file_in.close();}catch(Exception exdumb){}
      System.out.println(line+" lines read");
      first = true;
      found = 0;
      return found;
    }
    abcu[0] = a;
    abcu[1] = b;
    abcu[2] = c;
    abcu[3] = u;
    return found;
  } // end data_set3d

  int data_set4d(double abcdu[])
  {
    boolean debug = false;
    double a=0.0, b=0.0, c=0.0, d=0.0, u=0.0;
    int index, last, len; // for parsing a, b, c, d, u
    try
    {
      if(first)
      {
        handle  = new FileReader(data_file_name);
        file_in = new BufferedReader(handle);
        input_line = file_in.readLine();
        line = 0;
        found = 1;
        first = false;
      }
      if(input_line != null)
      {
	if(debug) System.out.println(input_line);
        // process input
	index = 0; // start of u
	if(debug) System.out.println("start u index="+index);
	last = input_line.indexOf(' ',index+1); // end of u
	if(debug) System.out.println("last of u ="+last);
	u = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("u="+u);
	index = input_line.indexOf(' ',last); // start of a
	if(debug) System.out.println("start a index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of a ="+last);
	a = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("a="+a);
	index = input_line.indexOf(' ',last); // start of b
	if(debug) System.out.println("start b index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of b ="+last);
	b = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("b="+b);
	index = input_line.indexOf(' ',last); // start of c
	if(debug) System.out.println("start c index="+index);
	last = input_line.indexOf(' ',index+1);
	if(debug) System.out.println("last of c ="+last);
	c = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("c="+c);
	index = input_line.indexOf(' ',last); // start of d
	if(debug) System.out.println("start d index="+index);
	last = input_line.length();
	if(debug) System.out.println("last of d ="+last);
	d = Double.parseDouble(input_line.substring(index+1,last));
	if(debug) System.out.println("d="+d);
        input_line = file_in.readLine();
        line++;
      }
      else
      {
        file_in.close();
        System.out.println(line+" lines read");
        first = true;
        found = 0;
        return found;
      }
    }
    catch(Exception exception) // catches all exceptions
    {
      // do nothing, keep going
      System.out.println("file read exception");
    }
    abcdu[0] = a;
    abcdu[1] = b;
    abcdu[2] = c;
    abcdu[3] = d;
    abcdu[4] = u;
    return found;
  } // end data_set_4d


  void check_1d(int n, int m, int nnn[], double C[])
  {
    double a1, u1, ua, diff;
    double sumsq = 0.0;
    double rms_err;
    double avg_err;
    double max_err;
    double sum = 0.0;
    double maxe = 0.0;
    double maxet = 0.0;
    double amin=0.0, amax=0.0, umin=0.0, umax=0.0;
    double Av[] = new double[200]; /* at least n */
    double au[] = new double[2];
    int i, k;
    int nn = nnn[0];
    boolean debug = false;

    k = 0;
    while(data_set1d(au)==1)
    {
      a1 = au[0];
      u1 = au[1];
      Av[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*a1;
      }
      if(k==0)
      {
	amin=a1; amax=a1;
	umin=u1; umax=u1;
      }
      if(a1>amax) amax=a1;
      if(a1<amin) amin=a1;
      if(u1>umax) umax=u1;
      if(u1<umin) umin=u1;

      ua = 0.0;
      for(i=0; i<nn; i++)
      {
	if(Math.abs(C[i])>tol)
	{
	  ua = ua + Av[Apwr[i]] * C[i];
	}
      }
      diff = Math.abs(u1-ua);
      if(debug) System.out.println("sample="+k+", u_fit="+ua+", u_data="+u1+
                                   ", diff="+(u1-ua));
      if(diff>maxe) maxe=diff;
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
      if(debug) // compare to all terms
      {
        ua = 0.0;
        for(i=0; i<nn; i++)
        {
	  ua = ua + Av[Apwr[i]] * C[i];
        }
        diff = Math.abs(u1-ua);
        if(debug) System.out.println("all   ="+k+", u_fit="+ua+", u_data="+u1+
                                   ", diff="+(u1-ua));
        if(diff>maxet) maxet=diff;
      } // end debug
      k++;
    }
    System.out.println("check_1d pts="+k+", amin="+amin+", amax="+amax+
                       ", umin="+umin+", umax="+umax);
    max_err = maxe;
    avg_err = sum/(double)k;
    rms_err = Math.sqrt(sumsq/(double)k);
    System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                       ", max_err="+max_err);
    if(debug)System.out.println("max_err with all ="+maxet);
    System.out.println(" ");
  } // end check_1d

  void check_2d(int n, int m, int nnn[], double C[])
  {
    double a1, b1, u1, ua, diff;
    double sumsq = 0.0;
    double rms_err;
    double avg_err;
    double max_err;
    double sum = 0.0;
    double maxe = 0.0;
    double amin=0.0, amax=0.0, bmin=0.0, bmax=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 abu[] = new double[3];
    int i, k;
    int nn = nnn[0];
    boolean debug = false;

    k = 0;
    while(data_set2d(abu)==1)
    {
      a1 = abu[0];
      b1 = abu[1];
      u1 = abu[2];
      Av[0] = 1.0;
      Bv[0] = 1.0;
      for(i=1; i<=n; i++)
      {
	Av[i] = Av[i-1]*a1;
	Bv[i] = Bv[i-1]*b1;
      }
      if(k==0)
      {
	amin=a1; amax=a1;
	bmin=b1; bmax=b1;
	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(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]] * C[i];
      }
      diff = Math.abs(u1-ua);
      if(debug) System.out.println("sample="+k+", u_fit="+ua+", u_data="+u1+
                                   ", diff="+(u1-ua));
      if(diff>maxe) maxe=diff;
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
      k++;
    }
    System.out.println("check_2d pts="+k+", amin="+amin+", amax="+amax+
                       ", bmin="+bmin+", bmax="+bmax);
    System.out.println("                    umin="+umin+", umax="+umax);
    max_err = maxe;
    avg_err = sum/(double)k;
    rms_err = Math.sqrt(sumsq/(double)k);
    System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                       ", max_err="+max_err);
    System.out.println(" ");
  } // end check_2d


  void check_3d(int n, int m, int nnn[], double C[])
  {
    double a1, b1, c1, u1, ua, diff;
    double sumsq = 0.0;
    double sum = 0.0;
    double rms_err;
    double avg_err;
    double max_err;
    double maxe = 0.0;
    double amin=0.0, amax=0.0, bmin=0.0, bmax=0.0, cmin=0.0;
    double cmax=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];
    int i, k;
    int nn = nnn[0];
    double abcu[] = new double[4];
  
    k = 0;
    while(data_set3d(abcu)==1)
    {
      a1=abcu[0];
      b1=abcu[1];
      c1=abcu[2];
      u1=abcu[3];
      Av[0] = 1.0;
      Bv[0] = 1.0;
      Cv[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;
      }

      if(k==0)
      {
	amin=a1; amax=a1;
	bmin=b1; bmax=b1;
	cmin=c1; cmax=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(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]] * C[i];
      }
      diff = Math.abs(u1-ua);
      if(diff>maxe) maxe=diff;
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
      k++;
    }
    System.out.println("check_3d k="+k+", amin="+amin+", amax="+amax+
                       ", bmin="+bmin+", bmax="+bmax);
    System.out.println("         cmin="+cmin+", cmax="+cmax+
                       ", umin="+umin+", umax="+umax);
    max_err = maxe;
    avg_err = sum/(double)k;
    rms_err = Math.sqrt(sumsq/(double)k);
    System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                       ", max_err="+max_err);
    System.out.println(" ");
  } // end check_3d

  void check_4d(int n, int mm, int nnn[], double C[])
  {
    double a1, b1, c1, d1, u1, ua, diff;
    double sumsq = 0.0;
    double sum = 0.0;
    double rms_err;
    double avg_err;
    double max_err;
    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];
  
    k = 0;
    while(data_set4d(abcdu)==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]] * C[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);
    System.out.println("rms_err="+rms_err+", avg_err="+avg_err+
                       ", max_err="+max_err);
    System.out.println(" ");
  } // end check_4d 

  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, value set to zero
      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)
  {
    if(args.length<3)
    {
      System.out.println("fit_data_file running");
      System.out.println("provide on command line:");
      System.out.println("number-of-variables 1,2,3 or 4");
      System.out.println("maximum order of polynomial");
      System.out.println("file name");
      System.exit(1);
    }
    int vars = Integer.parseInt(args[0]);
    int order = Integer.parseInt(args[1]);
    new fit_data_file(vars, order, args[2]);
  } // end main

} // end class fit_data_file