/* test_lsfit.c  2 dimension, n = 2, 3, 4, 5 fit */
#include <stdio.h>
#include <stdlib.h>
#include "lsfit.h"

#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))
#undef max
#define max(x,y) ((x)>(y)?(x):(y))

int main(int argc, char *argv[])
{
  int m; /* dimension */
  int n; /* order, max power */
  double x, y, z, t, u, v, w;
  double xg[6] = {0.1,  0.5,  1.2,  1.9,  2.7,  3.8};
  double yg[6] = {0.2,  0.6,  1.3,  2.1,  2.9,  4.1};
  double zg[6] = {0.25, 0.55, 1.22, 2.05, 2.6,  3.9};
  double tg[6] = {0.15, 0.45, 1.32, 1.95, 2.8,  3.7};
  double ug[6] = {0.08, 0.43, 1.42, 1.85, 2.95, 3.5};
  double vg[6] = {0.07, 0.41, 1.52, 1.65, 2.85, 3.6};
  double wg[6] = {0.09, 0.40, 1.40, 2.70, 3.30, 3.4};
  int i, j, k, l, nx, ny, nz, nt, nu, nv, nw, last;
  int ii, iii, iiii, iiiii, iiiiii, iiiiiii;
  double u1, v1, w1, err, diff, maxerr;

  nx = 6;
  ny = 6;
  nz = 6;
  nt = 6;
  nu = 6;
  nv = 6;
  nw = 6;
  m = 2; /* 2D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      u = x*x+2.0*x*y+3.0*y*y+4.0*x+5.0*y+6.0;
      if(i==nx-1 && j==ny-1) last = 1;
      fit2_data(x, y, u, last);
    }
  }
  /* check */
  printf("2D nx=%d, ny=%d \n", nx, ny);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      u1 = eval2_poly(x, y);
      u = x*x+2.0*x*y+3.0*y*y+4.0*x+5.0*y+6.0;
      err = abs(u1-u);
      maxerr = max(maxerr, err);
      printf("x=%f, y=%f, u=%f, eval=%f, err=%15.7e\n",
             x, y, u, u1, err);
    }
  }
  printf("2D maxerr = %g \n\n", maxerr);
  fit_free();

  m = 3; /* 3D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      for(k=0; k<nz; k++)
      {
        z = zg[k];
        u = x*x+2.0*y*y+3.0*z*z+4.0*x*y+5.0*x*z+6.0*y*z+
            7.0*x+8.0*y+9.0*z+10.0;
        if(i==nx-1 && j==ny-1 && k==nz-1) last = 1;
        fit3_data(x, y, z, u, last);
      }
    }
  }
  /* check */
  printf("3D nx=%d, ny=%d, nz=%d \n", nx, ny, nz);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      for(k=0; k<nz; k++)
      {
        z = zg[k];
        u1 = eval3_poly(x, y, z);
        u = x*x+2.0*y*y+3.0*z*z+4.0*x*y+5.0*x*z+6.0*y*z+
            7.0*x+8.0*y+9.0*z+10.0;
        err = abs(u1-u);
        maxerr = max(maxerr, err);
        /* printf("x=%f, y=%f, z=%f, u=%f, eval=%f, err=%15.7e\n",
	   x, y, z, u, u1, err); */
      }
    }
  }
  printf("3D maxerr = %g \n\n", maxerr);
  fit_free();

  m = 4; /* 4D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      for(k=0; k<nz; k++)
      {
        z = zg[k];
        for(l=0; l<nt; l++)
        {
          t = tg[l];
          u = 1.0+2.0*x+3.0*y+4.0*z+5.0*t +
              6.0*x*x+7.0*x*y+8.0*y*y+9.0*x*z +
              10.0*y*z+11.0*z*z+12.0*x*t+13.0*y*t +
              14.0*z*t+15.0*t*t;
          if(i==nx-1 && j==ny-1 && k==nz-1 && l==nt-1) last = 1;
          fit4_data(x, y, z, t, u, last);
	}
      }
    }
  }
  /* check */
  printf("4D nx=%d, ny=%d, nz=%d, nt=%d \n", nx, ny, nz, nt);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(j=0; j<ny; j++)
    {
      y = yg[j];
      for(k=0; k<nz; k++)
      {
        z = zg[k];
        for(l=0; l<nt; l++)
        {
          t = tg[l];
          u1 = eval4_poly(x, y, z, t);
          u = 1.0+2.0*x+3.0*y+4.0*z+5.0*t +
              6.0*x*x+7.0*x*y+8.0*y*y+9.0*x*z +
              10.0*y*z+11.0*z*z+12.0*x*t+13.0*y*t +
              14.0*z*t+15.0*t*t;
          err = abs(u1-u);
          maxerr = max(maxerr, err);
          /* printf("x=%f, y=%f, z=%f, t=%f, u=%f, eval=%f, err=%15.7e\n",
	     x, y, z, t, u, u1, err); */
	}
      }
    }
  }
  printf("4D maxerr = %g \n\n", maxerr);
  fit_free();

  m = 5; /* 5D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            w = 1.0+2.0*x+3.0*y+4.0*z+5.0*t + 6.0*u +
                7.0*x*x+8.0*x*y+9.0*y*y+10.0*x*z + 11.0*x*u +
                12.0*y*u+13.0*z*z+14.0*x*t+15.0*y*t +
                16.0*z*t+17.0*t*t+18.0*u*u;
            if(i==nx-1 && ii==ny-1 && iii==nz-1 && iiii==nt-1 && iiiii==nu-1)
               last = 1;
            fit5_data(x, y, z, t, u, w, last);
	  }
	}
      }
    }
  }
  /* check */
  printf("5D nx=%d, ny=%d, nz=%d, nt=%d, nu=%d \n",
         nx, ny, nz, nt, nu);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            w1 = eval5_poly(x, y, z, t, u);
            w = 1.0+2.0*x+3.0*y+4.0*z+5.0*t + 6.0*u +
                7.0*x*x+8.0*x*y+9.0*y*y+10.0*x*z + 11.0*x*u +
                12.0*y*u+13.0*z*z+14.0*x*t+15.0*y*t +
                16.0*z*t+17.0*t*t+18.0*u*u;
            err = abs(w1-w);
            maxerr = max(maxerr, err);
            /* printf("x=%f, y=%f, z=%f, t=%f, u=%f, w=%f, eval=%f, err=%15.7e\n",
	       x, y, z, t, u, w, w1, err); */ 
	  }
	}
      }
    }
  }
  printf("5D maxerr = %g \n\n", maxerr);
  fit_free();

  m = 6; /* 6D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=0; iiiiii<nv; iiiiii++)
            {
              v = vg[iiiii];
              w = 1.0+2.0*x+3.0*y+4.0*z+5.0*t + 6.0*u + 7.0*v +
                  8.0*x*x+9.0*x*y+10.0*y*y+11.0*x*z+12.0*x*u+13.0*x*v+
                  14.0*y*u+15.0*z*z+16.0*x*t+17.0*y*t +
                  18.0*z*t+19.0*t*t+20.0*u*u+21.0*v*v;
              if(i==nx-1 && ii==ny-1 && iii==nz-1 && iiii==nt-1 && 
                 iiiii==nu-1 && iiiiii==nv-1)
                last = 1;
              fit6_data(x, y, z, t, u, v, w, last);
	    }
	  }
	}
      }
    }
  }
  /* check */
  printf("6D nx=%d, ny=%d, nz=%d, nt=%d, nu=%d, nv=%d \n",
         nx, ny, nz, nt, nu, nv);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=0; iiiiii<nv; iiiiii++)
            {
              v = vg[iiiii];
              w1 = eval6_poly(x, y, z, t, u, v);
              w = 1.0+2.0*x+3.0*y+4.0*z+5.0*t + 6.0*u + 7.0*v +
                  8.0*x*x+9.0*x*y+10.0*y*y+11.0*x*z+12.0*x*u+13.0*x*v+
                  14.0*y*u+15.0*z*z+16.0*x*t+17.0*y*t +
                  18.0*z*t+19.0*t*t+20.0*u*u+21.0*v*v;
              err = abs(w1-w);
              maxerr = max(maxerr, err);
              /* printf("x=%f, y=%f, z=%f, t=%f, u=%f, w=%f, eval=%f, err=%15.7e\n",
	         x, y, z, t, u, w, w1, err); */ 
	    }
	  }
	}
      }
    }
  }
  printf("6D maxerr = %g \n\n", maxerr);
  fit_free();

  // 7D
  m = 7; /* 7D */
  n = 2; /* max power */
  fit_init(n,m);
  last = 0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=0; iiiiii<nv; iiiiii++)
            {
              v = vg[iiiiii];
              for(iiiiiii=0; iiiiiii<nw; iiiiiii++)
              {
                w = wg[iiiiiii];
                v1 = 1.0+2.0*x+3.0*y+4.0*z+5.0*t+6.0*u+7.0*v+8.0*w +
                     9.0*x*x+10.0*x*y+11.0*y*y+12.0*x*z+13.0*x*u+14.0*x*v+
                    15.0*y*u+16.0*z*z+17.0*x*t+18.0*y*t +
                    19.0*z*t+20.0*t*t+21.0*u*u+22.0*v*v+23.0*w*w+
		    24.0*x*w+25*y*w+26.0*z*w+27.0*u*w+28.0*v*w;
                if(i==nx-1 && ii==ny-1 && iii==nz-1 && iiii==nt-1 && 
                   iiiii==nu-1 && iiiiii==nv-1 && iiiiiii==nw-1)
                  last = 1;
                fit7_data(x, y, z, t, u, v, w, v1, last);
	      }
	    }
	  }
	}
      }
    }
  }
  /* check 7 */
  printf("7D nx=%d, ny=%d, nz=%d, nt=%d, nu=%d, nv=%d, nw=%d \n",
         nx, ny, nz, nt, nu, nv, nw);
  maxerr = 0.0;
  for(i=0; i<nx; i++)
  {
    x = xg[i];
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      for(iii=0; iii<nz; iii++)
      {
        z = zg[iii];
        for(iiii=0; iiii<nt; iiii++)
        {
          t = tg[iiii];
          for(iiiii=0; iiiii<nu; iiiii++)
          {
            u = ug[iiiii];
            for(iiiiii=0; iiiiii<nv; iiiiii++)
            {
              v = vg[iiiiii];
              for(iiiiiii=0; iiiiiii<nw; iiiiiii++)
              {
                w = wg[iiiiiii];
		w1 = eval7_poly(x, y, z, t, u, v, w);
                v1 = 1.0+2.0*x+3.0*y+4.0*z+5.0*t+6.0*u+7.0*v+8.0*w+
                     9.0*x*x+10.0*x*y+11.0*y*y+12.0*x*z+13.0*x*u+14.0*x*v+
                    15.0*y*u+16.0*z*z+17.0*x*t+18.0*y*t +
                    19.0*z*t+20.0*t*t+21.0*u*u+22.0*v*v+23.0*w*w+
		    24.0*x*w+25*y*w+26.0*z*w+27.0*u*w+28.0*v*w;

                err = abs(w1-v1);
                maxerr = max(maxerr, err);
                /* printf("x=%f, y=%f, z=%f, t=%f, u=%f, w=%f, eval=%f, err=%15.7e\n",
	           x, y, z, t, u, w, w1, err); */
	      } 
	    }
	  }
	}
      }
    }
  }
  printf("7D maxerr = %g \n\n", maxerr);
  fit_free();


  printf("test_lsfit.c finished\n");
  return 0;
} /* end test_lsfit.c */