// test_lsfit4.java

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

public class test_lsfit4
{
  int nx = 5;
  int ny = 5;
  int nz = 5;
  int nt = 5;
  int dof = (nx-2)*(ny-2)*(nz-2)*(nt-2);  // Degrees of Freedom
  double hx, hy, hz, ht;
  double xg[] = new double[nx]; // X grid, does not need to be uniform
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double tg[] = new double[nt];
  double ug[][][][] = new double[nx][ny][nz][nt]; 
  double uA[][][][] = new double[nx][ny][nz][nt]; 
  double xmin = -1.0;
  double xmax =  1.0;
  double ymin = -1.0;
  double ymax =  1.0;
  double zmin = -1.0;
  double zmax =  1.0;
  double tmin = -1.0;
  double tmax =  1.0;
  double x, y, z, t, U, Ue;
  double err, maxerr;
  int ierr = 0;
  boolean last = false;
  DecimalFormat f4 = new DecimalFormat("0.000");

  public test_lsfit4()
  {
    System.out.println("test_lsfit4.java");
    System.out.println(" ");
    // initialize hx, hy, hy, ht, xg, yg, zg, tg
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    hz = (zmax-zmin)/(double)(nz-1);
    ht = (tmax-tmin)/(double)(nt-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 iii=0; iii<nz; iii++) zg[iii] = zmin+hz*(double)iii;
    for(int iiii=0; iiii<nt; iiii++) tg[iiii] = tmin+ht*(double)iiii;
    System.out.println("nx="+nx+", ny="+ny+", nz="+nz+", nt="+nt);
    System.out.println("hx="+hx+", hy="+hy+", hz="+hz+", ht="+ht);


    System.out.println("test_lsfit4.java (4,1) data fit");
    lsfit LS1 = new lsfit(4, 1, dof); // dimension, power, points
    LS1.debug = false;
    last = false;
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
            U = u1(x,y,z,t);
	    if(i==nx-1 && ii==ny-1 && iii==nz-1 && iiii==nt-1) last = true;
	    LS1.fit4_data(x, y, z, t, U, last);
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java (4,1) data fit coefficients");
    for(int i=0; i<LS1.nnn; i++)
	System.out.println("C["+i+"]="+LS1.C[i]+"  x^"+LS1.Apwr[i]+
                           " y^"+LS1.Bpwr[i]+" z^"+LS1.Cpwr[i]+
                           " t^"+LS1.Dpwr[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];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
            Ue = u1(x, y, z, t);
	    U = LS1.eval4_poly(x, y, z, t);
            err = Math.abs(U-Ue);
            if(err>maxerr) ierr = 1000*i+100*ii+10*iii+iiii;
            maxerr = Math.max(maxerr,err);
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java (4,1) data fit maxerr="+maxerr+
                       " at "+ierr);
    System.out.println(" ");

    System.out.println("test_lsfit4.java (4,2) data fit");
    lsfit LS2 = new lsfit(4, 2, dof); // dimension, power, points
    LS2.debug = false;
    last = false;
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
            U = u2(x,y,z,t);
	    if(i==nx-1 && ii==ny-1 && iii==nz-1 && iiii==nt-1) last = true;
	    LS2.fit4_data(x, y, z, t, U, last);
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java (4,2) data fit coefficients");
    for(int i=0; i<LS2.nnn; i++)
	System.out.println("C["+i+"]="+LS2.C[i]+"  x^"+LS2.Apwr[i]+
                           " y^"+LS2.Bpwr[i]+" z^"+LS2.Cpwr[i]+
                           " t^"+LS2.Dpwr[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];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
            Ue = u2(x, y, z, t);
	    U = LS2.eval4_poly(x, y, z, t);
            err = Math.abs(U-Ue);
            maxerr = Math.max(maxerr,err);
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java (4,2) data fit maxerr="+maxerr);
    System.out.println(" ");

    // set up test data
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
	    uA[i][ii][iii][iiii] = u2(x,y,z,t);
	    if(i==0 || ii==0 || i==nx-1 || ii==ny-1 || iii==0 || iii==nz-1 ||
               iiii==0 || iiii==nt-1)
	    {
		ug[i][ii][iii][iiii] = uA[i][ii][iii][iiii];
	    }
	    else // not boundary
	    {
	      ug[i][ii][iii][iiii] = 0.0;
	    }
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("Ua, ug initialized");

    // boundary fit
    lsfit LSb = new lsfit(4, 2, dof); // dimension, power, points
    LSb.debug = false;
    LSb.fit4_bpoly(nx, ny, nz, nt, xg, yg, zg, tg, ug);
    System.out.println("test_lsfit4.java (4,2) data fit coefficients");
    for(int i=0; i<LSb.nnn; i++)
	System.out.println("C["+i+"]="+LSb.C[i]+"  x^"+LSb.Apwr[i]+
                           " y^"+LSb.Bpwr[i]+" z^"+LSb.Cpwr[i]+
                           " t^"+LSb.Dpwr[i]);

    // check boundary
    for(int i=0; i<nx; i++)
    {
      x = xg[i];
      for(int ii=0; ii<ny; ii++)
      {
	y = yg[ii];
	for(int iii=0; iii<nz; iii++)
	{
	  z = zg[iii];
	  for(int iiii=0; iiii<nt; iiii++)
	  {
	    t = tg[iiii];
	    if(i==0 || ii==0 || i==nx-1 || ii==ny-1 || iii==0 || iii==nz-1 ||
               iiii==0 || iiii==nt-1)
	    {
	      U = LSb.eval4_poly(x, y, z, t);
              err = Math.abs(U-uA[i][ii][iii][iiii]);
              maxerr = Math.max(maxerr,err);
	    }
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java, on boundary maxerr="+maxerr);

    // inside boundary fit
    // check boundary
    for(int i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(int ii=1; ii<ny-1; ii++)
      {
	y = yg[ii];
	for(int iii=1; iii<nz-1; iii++)
	{
	  z = zg[iii];
	  for(int iiii=1; iiii<nt-1; iiii++)
	  {
	    t = tg[iiii];
	    U = LSb.eval4_poly(x, y, z, t);
            err = Math.abs(U-uA[i][ii][iii][iiii]);
            maxerr = Math.max(maxerr,err);
	  } // iiii
	} // iii
      } // ii
    } // i
    System.out.println("test_lsfit4.java inside boundary, maxerr="+maxerr);

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

  double u1(double x, double y, double z, double t)
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*z + 5.0*t;
  }

  double u2(double x, double y, double z, double t)
  {
    return 1.0 + 2.0*x + 3.0*y + 4.0*z + 5.0*t +
           6.0*x*x + 7.0*x*y + 8.0*x*z + 9.0*x*t +
          10.0*y*y + 11.0*y*z + 12.0*y*t + 13.0*z*z +
          14.0*z*t + 15.0*t*t;
  }

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