// test_nuderiv2dg.java   test  nuderiv2dg  with 28 points 
//          order 1  minimum independent points 3   Ux ... on first degree
//          order 2  minimum independent points 6   Uxx ... on second degree
//          order 3  minimum independent points 10  Uxxx ... on third degree
//          order 4  minimum independent points 15  Uxxxx ... on fourth degree
//          order 5  minimum independent points 21  Uxxxxx ... on fifth degree
//          going beyond minimum points does not help much 
class test_nuderiv2dg
{
  test_nuderiv2dg()
  {
    int debug = 0;
    int m;
    int n = 21; // max 5th order this test 
    int ordpt[] = {0, 3, 6, 10, 15, 21, 28};
    double x[] = {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[] = {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};
    double U[] = new double[28];
    double est, actual, dx, dy;
    double cx[] = new double[28];
    double cy[] = new double[28];
    double cxx[] = new double[28];
    double cxy[] = new double[28];
    double cyy[] = new double[28];
    double cxxx[] = new double[28];
    double cxxy[] = new double[28];
    double cxyy[] = new double[28];
    double cyyy[] = new double[28];
    double cxxxx[] = new double[28];
    double cxxxy[] = new double[28];
    double cxxyy[] = new double[28];
    double cxyyy[] = new double[28];
    double cyyyy[] = new double[28];
    int ip[] = new int[28];
    double error, maxerror;

    nuderiv2dg N = new nuderiv2dg();

    System.out.println("test_nuderiv2dg.java running");
    for(int i=0; i<n; i++)
    {
      U[i] = f(x[i],y[i]);
      dx = fx(x[i],y[i]);
      dy = fy(x[i],y[i]);
      System.out.println("point "+i+"  x="+x[i]+", y="+y[i]+
                         ", U="+U[i]+", dx="+dx+", dy="+dy);
    } 

    // f(x,y) is 4th order, thus 4 is minimum order for this test
    for(int order=4; order<=5; order++)
    {
      n = ordpt[order];
      maxerror = 0.0;
      for(int point=0; point<n; point++) // Ux 
      {
        m = N.nu2dx(order, n, point, x, y, cx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cx[i]*U[ip[i]];
        actual = fx(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Ux 

      for(int point=0; point<n; point++) // Uy 
      {
        m = N.nu2dy(order, n, point, x, y, cy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cy[i]*U[ip[i]];
        actual = fy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uy 

      if(n<6) continue;
      for(int point=0; point<n; point++) // Uxx 
      {
        m = N.nu2dxx(order, n, point, x, y, cxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxx[i]*U[ip[i]];
        actual = fxx(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxx 

      for(int point=0; point<n; point++) // Uxy 
      {
        m = N.nu2dxy(order, n, point, x, y, cxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxy[i]*U[ip[i]];
        actual = fxy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxy 

      for(int point=0; point<n; point++) // Uyy 
      {
        m = N.nu2dyy(order, n, point, x, y, cyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyy[i]*U[ip[i]];
        actual = fyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyy 

      if(n<10) continue;
      for(int point=0; point<n; point++) // Uxxx 
      {
        m = N.nu2dxxx(order, n, point, x, y, cxxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxx[i]*U[ip[i]];
        actual = fxxx(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxx 

      for(int point=0; point<n; point++) // Uxxy 
      {
        m = N.nu2dxxy(order, n, point, x, y, cxxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxy[i]*U[ip[i]];
        actual = fxxy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxy 

      for(int point=0; point<n; point++) // Uxyy 
      {
        m = N.nu2dxyy(order, n, point, x, y, cxyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxyy[i]*U[ip[i]];
        actual = fxyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxyy 

      for(int point=0; point<n; point++) // Uyyy 
      {
        m = N.nu2dyyy(order, n, point, x, y, cyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyyy[i]*U[ip[i]];
        actual = fyyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyyy 

      if(n<15) continue;
      for(int point=0; point<n; point++) // Uxxxx 
      {
        m = N.nu2dxxxx(order, n, point, x, y, cxxxx, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxxx[i]*U[ip[i]];
        actual = fxxxx(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxxx  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxxx 

      for(int point=0; point<n; point++) // Uxxxy 
      {
        m = N.nu2dxxxy(order, n, point, x, y, cxxxy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxxy[i]*U[ip[i]];
        actual = fxxxy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxxy %2d,  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxxy 

      for(int point=0; point<n; point++) // Uxxyy 
      {
        m = N.nu2dxxyy(order, n, point, x, y, cxxyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxxyy[i]*U[ip[i]];
        actual = fxxyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxxyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxxyy 

      for(int point=0; point<n; point++) // Uxyyy 
      {
        m = N.nu2dxyyy(order, n, point, x, y, cxyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cxyyy[i]*U[ip[i]];
        actual = fxyyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cxyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uxyyy 

      for(int point=0; point<n; point++) // Uyyyy 
      {
        m = N.nu2dyyyy(order, n, point, x, y, cyyyy, ip);
        System.out.println("order="+order+", n="+n+", point="+point);
        est = 0.0;
        for(int i=0; i<m; i++) est += cyyyy[i]*U[ip[i]];
        actual = fyyyy(x[point],y[point]);
        error = actual-est;
        maxerror = Math.max(maxerror,Math.abs(error));
        System.out.println("cyyyy  est="+est+", actual="+actual+
                           ", error="+ error);
      } // end Uyyyy 

      System.out.println("maxerror for order="+order+", n="+n+
                         ", is "+maxerror);
    } // end order loop 
    System.out.println("test_nuderiv2dg.java finished");
  }

  // test functions
  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;
  }

  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;
  }

  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;
  }

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

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

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

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

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

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

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

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

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

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

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

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

  public static void main (String[] args)
  {
    new test_nuderiv2dg();
  }
} // end  class test_nuderiv2dg