// test_lsfit2.java

import java.io.*;
import java.text.*;

public class test_lsfit2
{
  int nx = 21;
  int ny = 23;
  int dof = (nx-2)*(ny-2);  // Degrees of Freedom                        
  double hx, hy;
  double xg[] = new double[nx]; // X grid, does not need to be uniform
  double yg[] = new double[ny];
  double us[] = new double[dof];
  double ug[][] = new double[nx][ny]; 
  double uA[][] = new double[nx][ny]; 
  double xmin = -1.0;
  double xmax =  1.0;
  double ymin = -1.0;
  double ymax =  1.0;
    double x, y, u, U, err, maxerr;
  boolean last = false;
  DecimalFormat f4 = new DecimalFormat("0.000");

  public test_lsfit2()
  {
    System.out.println("test_lsfit2.java");
    System.out.println("fitting");
    System.out.println("6D expect 1.0 to 28.0 coefficients");
    System.out.println(" ");
    // initialize hx, hy, xg, yg  
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    for(int i=0; i<nx; i++) xg[i] = xmin+hx*(double)i;
    for(int ii=0; ii<ny; ii++) yg[ii] = ymin+hy*(double)ii;

    for(int k=2; k<=6; k++)
    {
      for(int i=0; i<nx; i++)
      {
        x = xg[i];
        for(int ii=0; ii<ny; ii++)
        {
	  y = yg[ii];
          if(k==2) U = u2(x,y);
          if(k==3) U = u3(x,y);
          if(k==4) U = u4(x,y);
          if(k==5) U = u5(x,y);
          if(k==6) U = u6(x,y);
	  uA[i][ii] = U;
        } // ii
      } // i
      System.out.println("ug initialized\n");

      System.out.println("test_lsfit2.java on data 2,"+k);
      lsfit LSd = new lsfit(2, k, dof); // dimension, power, points
      LSd.debug = false; // to turn off debug printing in lsfit.java
      last = false;
      for(int i=0; i<nx; i++)
      {
        x = xg[i];
        for(int ii=0; ii<ny; ii++)
        {
	  y = yg[ii];
	  u = uA[i][ii];
          if(i == nx-1 && ii== ny-1) last = true;
          LSd.fit2_data(x, y, u, last);
        } // ii
      } // i
      for(int i=0; i<LSd.nnn; i++)
        System.out.println("C["+i+"]="+LSd.C[i]+"  x^"+LSd.Apwr[i]+
                                                " y^"+LSd.Bpwr[i]);
      // check
      maxerr = 0.0;
      for(int i=0; i<nx; i++)
      {
        x = xg[i];
        for(int ii=0; ii<ny; ii++)
        {
	  y = yg[ii];
	  u = LSd.eval2_poly(x, y);
          err = Math.abs(u-uA[i][ii]);
	  maxerr = Math.max(maxerr,err);
        } // ii
      } // i
      System.out.println("lsfit error on data, maxerr="+maxerr);
    } // k
    System.out.println("end test_lsfit2.java on data");
    System.out.println(" ");

    System.out.println("test_lsfit2.java on bpoly");
    System.out.println("nx="+nx+", ny="+ny+" 4D");
    lsfit LS = new lsfit(2, 4, dof); // dimension, power, points
    LS.debug = false; // to turn off debug printing in lsfit.java
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
        uA[i][ii] = u4(x,y);
        if(i==0 || i==nx-1 || ii==0 || ii==ny-1)
	  ug[i][ii] = uA[i][ii];
        else 
          ug[i][ii] = 0.0;
      } // ii
    } // i
    LS.fit2_bpoly(nx, ny, xg, yg, ug);
    for(int i=0; i<LS.nnn; i++)
      System.out.println("C["+i+"]="+LS.C[i]+"  x^"+LS.Apwr[i]+
                                              " y^"+LS.Bpwr[i]);
    // check boundary
    maxerr = 0.0;
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
        if(i==0 || i==nx-1 || ii==0 || ii==ny-1)
	{
	  u = LS.eval2_poly(x, y);
          err = Math.abs(u-uA[i][ii]);
	  maxerr = Math.max(maxerr,err);
          //System.out.println("want="+f4.format(uA[i][ii])+
          //                   "  fit="+f4.format(u)+
	  //		     "  err="+(u-uA[i][ii]));
	} // end if
      } // ii
    } // i
    System.out.println("test_lsfit2.java on boundary bpoly, maxerr="+maxerr);
    // check inside, using same coefficients
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
	y = yg[ii];
	u = LS.eval2_poly(x, y);
        err = Math.abs(u-uA[i][ii]);
	maxerr = Math.max(maxerr,err);
        //System.out.println("want="+f4.format(uA[i][ii])+
        //                   "  fit="+f4.format(u)+
	//		     "  err="+(u-uA[i][ii]));
      } // ii
    } // i
    System.out.println("test_lsfit2.java on inside bpoly, maxerr="+maxerr);
    System.out.println("end test_lsfit2.java on bpoly");

    System.out.println("finished test2_lsfit.java");
  } // end test_lsfit2

  double u2(double x, double y) // all possible 6 terms
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y;
  }

  double u3(double x, double y) // all possible 10 terms
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y +
           7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y + 10.0*y*y*y;
  }

  double u4(double x, double y) // all possible 15 terms
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y +
           7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y + 10.0*y*y*y +
          11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 
	  14.0*x*y*y*y + 15.0*y*y*y*y;
  }

  double u5(double x, double y) // all possible 21 terms
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y +
           7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y + 10.0*y*y*y +
          11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 
          14.0*x*y*y*y + 15.0*y*y*y*y +
          16.0*x*x*x*x*x + 17.0*x*x*x*x*y + 18.0*x*x*x*y*y + 
	  19.0*x*x*y*y*y + 20.0*x*y*y*y*y + 21.0*y*y*y*y*y;
  }

  double u6(double x, double y) // all possible 28 terms
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*x*x + 5.0*x*y + 6.0*y*y +
           7.0*x*x*x + 8.0*x*x*y + 9.0*x*y*y + 10.0*y*y*y +
          11.0*x*x*x*x + 12.0*x*x*x*y + 13.0*x*x*y*y + 
          14.0*x*y*y*y + 15.0*y*y*y*y +
          16.0*x*x*x*x*x + 17.0*x*x*x*x*y + 18.0*x*x*x*y*y + 
          19.0*x*x*y*y*y + 20.0*x*y*y*y*y + 21.0*y*y*y*y*y +
          22.0*x*x*x*x*x*x + 23.0*x*x*x*x*x*y + 24.0*x*x*x*x*y*y + 
          25.0*x*x*x*y*y*y + 26.0*x*x*y*y*y*y + 27.0*x*y*y*y*y*y +
	  28.0*y*y*y*y*y*y;
  }

  public static void main (String[] args)
  {
    new test_lsfit2();
  } // end main of test_lsfit2.java
} // end class test_lsfit2