/* test_nuderiv2dg.c   test  nuderiv2dg  with max 28 points */
/*                     order 1  minimum points 3       */
/*                     order 2  minimum points 6       */
/*                     order 3  minimum points 10      */
/*                     order 4  minimum points 15      */
/*                     order 5  minimum points 21      */
/*                     order 6  minimum points 28      */
/*                     going beyond minimum points does not help much  */
/*                     max order may be less if there are dependencies */

#include <stdio.h>
#include <stdlib.h>
#include "nuderiv2dg.h"
#include "inverse.h"

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

/* 4th order test case */
static double f(double x, double y)
{
  return x*x*x*x + 2.0*x*x*y*y + 3.0*y*y*y*y + 
       4.0*x*x*x + 5.0*x*x*y   + 6.0*y*y*y   +
       7.0*x*x   + 8.0*x*y     + 9.0*y*y     +
      10.0*x     +11.0*y       +12.0;
}

static double fx(double x, double y)
{
  return 4.0*x*x*x + 4.0*x*y*y +
        12.0*x*x   +10.0*x*y   +
        14.0*x     + 8.0*y     +
        10.0;
}

static double fy(double x, double y)
{
  return 4.0*x*x*y + 12.0*y*y*y +
         5.0*x*x   + 18.0*y*y   +
         8.0*x     + 18.0*y     +
        11.0;
}

static double fxx(double x, double y)
{
  return 12.0*x*x + 4.0*y*y +
         24.0*x   +10.0*y   +
         14.0;
}

static double fxy(double x, double y)
{
  return 8.0*x*y +
        10.0*x   +
        8.0;
}

static double fyy(double x, double y)
{
  return 4.0*x*x + 36.0*y*y +
        36.0*y              +
        18.0;
}

static double fxxx(double x, double y)
{
  return 24.0*x +
         24.0;
}

static double fxxy(double x, double y)
{
  return 8.0*y +
        10.0;
}

static double fxyy(double x, double y)
{
  return 8.0*x;
}

static double fyyy(double x, double y)
{
  return 72.0*y +
         36.0;
}

static double fxxxx(double x, double y)
{
  return 24.0;
}

static double fxxxy(double x, double y)
{
  return 0.0;
}

static double fxxyy(double x, double y)
{
  return 8.0;
}

static double fxyyy(double x, double y)
{
  return 0.0;
}

static double fyyyy(double x, double y)
{
  return 72.0;
}


/* 2nd order test case */
static double f2(double x, double y)
{
  return 1.0*x*x + 2.0*x*y + 3.0*y*y +
         4.0*x   + 5.0*y   + 6.0;
}

static double f2x(double x, double y)
{
  return  2.0*x + 2.0*y + 4.0;
}

static double f2y(double x, double y)
{
  return 2.0*x + 6.0*y + 5.0;
}

static double f2xx(double x, double y)
{
  return 2.0;
}

static double f2xy(double x, double y)
{
  return 2.0;
}

static double f2yy(double x, double y)
{
  return 6.0;
}

/* max order may be less due to dependencies */
static int xpow[] = {0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0,
                     5, 4, 3, 2, 1, 0, 6, 5, 4, 3, 2, 1, 0,
                     7, 6, 5, 4, 3, 2, 1, 0};
static int ypow[] = {0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4,
                     0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6,
                     0, 1, 2, 3, 4, 5, 6, 7};
static double xpwr[8], ypwr[8];


int main(int argc, char *argv[])
{
  int i, j;
  int debug = 0;
  int n = 28;
  int order = 4;
  double x[28] = {0.05, 0.14, 0.27, 0.42, 0.57, 0.71, 0.86,
                  0.04, 0.18, 0.31, 0.43, 0.63, 0.82, 0.95,
                  0.07, 0.16, 0.33, 0.41, 0.59, 0.82, 0.99,
                  0.11, 0.21, 0.35, 0.45, 0.65, 0.75, 0.97};
  double y[28] = {0.06, 0.05, 0.07, 0.10, 0.13, 0.15, 0.17,
                  0.30, 0.33, 0.36, 0.39, 0.31, 0.32, 0.37,
                  0.65, 0.77, 0.71, 0.67, 0.67, 0.75, 0.76,
                  0.92, 0.97, 0.93, 0.98, 0.87, 0.91, 0.99};
  int n2 = 9; /* will get cut down to n2i */
  int order2 = 2;
  double x2[9] = {1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0};
  double y2[9] = {1.0, 2.0, 3.0, 1.0, 2.0, 3.0, 1.0, 2.0, 3.0};
  double xa[36], ya[36]; /* for variable */
  int ip[36];
  int ixy[36], ni, nj, n2i, nd;
  double U[36], xg[36], yg[36]; /* for additional singularity testing */
  double est, actual, dx, dy;
  double cx[36], cy[36], cxx[36], cxy[36], cyy[36], cxxx[36], cxxy[36];
  double cxyy[36], cyyy[36], cxxxx[36], cxxxy[36], cxxyy[36], cxyyy[36];
  double cyyyy[36];
  double error, maxerror;
  int point;
  int stat = 0;

  printf("test_nuderiv2dg.c running \n");
  for(i=0; i<n; i++)
  {
    U[i] = f(x[i],y[i]);
    dx = fx(x[i],y[i]);
    dy = fy(x[i],y[i]);
    printf("point %2d  x=%5.3f, y=%5.3f, U=%7.3f, dx=%f, dy=%f \n",
           i, x[i], y[i], U[i], dx, dy);
  } 

  /* just one test case */
  {
    maxerror = 0.0;
    if(n<15) printf("n<15, thus expect error to grow, f(x,y) is 4th order\n");
    for(point=0; point<n; point++) /* Ux */
    {
      stat = nu2dx(order, n, point, x, y, cx, ip);
      printf("order=%1d, n=%2d, stat=%2d, point=%2d \n",
             order, n, stat, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cx[i]*U[ip[i]];
      actual = fx(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cx %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Ux */

    for(point=0; point<n; point++) /* Uy */
    {
      stat = nu2dy(order, n, point, x, y, cy, ip);
      printf("order=%1d, n=%2d, stat=%2d, point=%2d \n",
             order, n, stat, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cy[i]*U[ip[i]];
      actual = fy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uy */

    for(point=0; point<n; point++) /* Uxx */
    {
      stat = nu2dxx(order, n, point, x, y, cxx, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxx[i]*U[ip[i]];
      actual = fxx(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxx %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxx */

    for(point=0; point<n; point++) /* Uxy */
    {
      stat = nu2dxy(order, n, point, x, y, cxy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxy[i]*U[ip[i]];
      actual = fxy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxy */

    for(point=0; point<n; point++) /* Uyy */
    {
      stat = nu2dyy(order, n, point, x, y, cyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cyy[i]*U[ip[i]];
      actual = fyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uyy */

    for(point=0; point<n; point++) /* Uxxx */
    {
      stat = nu2dxxx(order, n, point, x, y, cxxx, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxxx[i]*U[ip[i]];
      actual = fxxx(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxxx %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxxx */

    for(point=0; point<n; point++) /* Uxxy */
    {
      stat = nu2dxxy(order, n, point, x, y, cxxy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxxy[i]*U[ip[i]];
      actual = fxxy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxxy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxxy */

    for(point=0; point<n; point++) /* Uxyy */
    {
      stat = nu2dxyy(order, n, point, x, y, cxyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxyy[i]*U[ip[i]];
      actual = fxyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxyy */

    for(point=0; point<n; point++) /* Uyyy */
    {
      stat = nu2dyyy(order, n, point, x, y, cyyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cyyy[i]*U[ip[i]];
      actual = fyyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cyyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uyyy */

    for(point=0; point<n; point++) /* Uxxxx */
    {
      stat = nu2dxxxx(order, n, point, x, y, cxxxx, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxxxx[i]*U[ip[i]];
      actual = fxxxx(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxxxx %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxxxx */

    for(point=0; point<n; point++) /* Uxxxy */
    {
      stat = nu2dxxxy(order, n, point, x, y, cxxxy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxxxy[i]*U[ip[i]];
      actual = fxxxy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxxxy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxxxy */

    for(point=0; point<n; point++) /* Uxxyy */
    {
      stat = nu2dxxyy(order, n, point, x, y, cxxyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxxyy[i]*U[ip[i]];
      actual = fxxyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxxyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxxyy */

    for(point=0; point<n; point++) /* Uxyyy */
    {
      stat = nu2dxyyy(order, n, point, x, y, cxyyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cxyyy[i]*U[ip[i]];
      actual = fxyyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cxyyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uxyyy */

    for(point=0; point<n; point++) /* Uyyyy */
    {
      stat = nu2dyyyy(order, n, point, x, y, cyyyy, ip);
      printf("order=%1d, n=%2d, point=%2d \n", order, n, point);
      est = 0.0;
      for(i=0; i<stat; i++) est += cyyyy[i]*U[ip[i]];
      actual = fyyyy(x[point],y[point]);
      error = actual-est;
      maxerror = max(maxerror,abs(error));
      printf("cyyyy %2d, est=%g, actual=%g, error=%g \n",
             point, est, actual, error);
    } /* end Uyyyy */

    printf("\nmaxerror for order=%d, n=%d, is %g \n\n", order, n, maxerror);
  } /* end while */

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