/* interpolation/test2d.c

 * 

 * Copyright 2012 David Zaslavsky

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

#include <config.h>

#include <math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_ieee_utils.h>

#include <gsl/gsl_test.h>

#include <gsl/gsl_interp.h>

#include <gsl/gsl_interp2d.h>

#include <gsl/gsl_spline2d.h>

/* tests a single evaluator function from the low-level interface */

static int

test_single_low_level(

  double (*evaluator)(const gsl_interp2d *, const double[], const double[],

                      const double[], const double, const double,

                      gsl_interp_accel *, gsl_interp_accel *),

  int (*evaluator_e)(const gsl_interp2d *, const double[], const double[],

                     const double[], const double, const double,

                     gsl_interp_accel *, gsl_interp_accel *, double *),

  const gsl_interp2d * interp, const double xarr[], const double yarr[],

  const double zarr[], const double x, const double y,

  gsl_interp_accel * xa, gsl_interp_accel * ya,

  const double expected_results[], size_t i)

{

  if (expected_results != NULL)

    {

      int status;

      double result = evaluator(interp, xarr, yarr, zarr, x, y, xa, ya);

      gsl_test_rel(result, expected_results[i], 1e-10,

                   "low level %s %d", gsl_interp2d_name(interp), i);

      status = evaluator_e(interp, xarr, yarr, zarr, x, y, xa, ya, &result);

      if (status == GSL_SUCCESS)

        {

          gsl_test_rel(result, expected_results[i], 1e-10,

                       "low level _e %s %d", gsl_interp2d_name(interp), i);

        }

    }

  return 0;

}

/* tests a single evaluator function from the high-level interface */

static int

test_single_high_level(

    double (*evaluator)(const gsl_spline2d *, const double, const double,

                        gsl_interp_accel *, gsl_interp_accel *),

    int (*evaluator_e)(const gsl_spline2d *, const double, const double,

                       gsl_interp_accel *, gsl_interp_accel *, double *),

    const gsl_spline2d * interp, const double x, const double y,

    gsl_interp_accel * xa, gsl_interp_accel * ya,

    const double expected_results[], size_t i)

{

  if (expected_results != NULL)

    {

      int status;

      double result = evaluator(interp, x, y, xa, ya);

      gsl_test_rel(result, expected_results[i], 1e-10,

                   "high level %s %d", gsl_spline2d_name(interp), i);

      status = evaluator_e(interp, x, y, xa, ya, &result);

      if (status == GSL_SUCCESS)

        {

          gsl_test_rel(result, expected_results[i], 1e-10,

                       "high level _e %s %d", gsl_spline2d_name(interp), i);

        }

    }

  return 0;

}

/*

 * Tests that a given interpolation type reproduces the data points

 * it is given, and then tests that it correctly reproduces additional

 * values.

 */

static int

test_interp2d(

const double xarr[], const double yarr[], const double zarr[], /* data */

              size_t xsize, size_t ysize,                /* array sizes */

              const double xval[], const double yval[],  /* test points */

              const double zval[],                       /* expected results */

              const double zxval[], const double zyval[],

              const double zxxval[], const double zyyval[], const double zxyval[],

              size_t test_size,                          /* number of test points */

              const gsl_interp2d_type * T)

{

  gsl_interp_accel *xa = gsl_interp_accel_alloc();

  gsl_interp_accel *ya = gsl_interp_accel_alloc();

  int status = 0;

  size_t xi, yi, zi, i;

  gsl_interp2d * interp = gsl_interp2d_alloc(T, xsize, ysize);

  gsl_spline2d * interp_s = gsl_spline2d_alloc(T, xsize, ysize);

  unsigned int min_size = gsl_interp2d_type_min_size(T);

  gsl_test_int(min_size, T->min_size, "gsl_interp2d_type_min_size on %s", gsl_interp2d_name(interp));

  gsl_interp2d_init(interp, xarr, yarr, zarr, xsize, ysize);

  gsl_spline2d_init(interp_s, xarr, yarr, zarr, xsize, ysize);

  /* First check that the interpolation reproduces the given points */

  for (xi = 0; xi < xsize; xi++)

    {

      double x = xarr[xi];

      for (yi = 0; yi < ysize; yi++)

        {

          double y = yarr[yi];

          zi = gsl_interp2d_idx(interp, xi, yi);

          test_single_low_level(&gsl_interp2d_eval, &gsl_interp2d_eval_e,

                                interp, xarr, yarr, zarr, x, y,

                                xa, ya, zarr, zi);

          test_single_low_level(&gsl_interp2d_eval_extrap,

                                &gsl_interp2d_eval_e_extrap, interp,

                                xarr, yarr, zarr, x, y, xa, ya, zarr, zi);

          test_single_high_level(&gsl_spline2d_eval, &gsl_spline2d_eval_e,

                                 interp_s, x, y, xa, ya, zarr, zi);

        }

    }

  /* Then check additional points provided */

  for (i = 0; i < test_size; i++)

    {

      double x = xval[i];

      double y = yval[i];

      test_single_low_level(&gsl_interp2d_eval,         &gsl_interp2d_eval_e,          interp, xarr, yarr, zarr, x, y, xa, ya, zval, i);

      test_single_low_level(&gsl_interp2d_eval_deriv_x, &gsl_interp2d_eval_deriv_x_e,  interp, xarr, yarr, zarr, x, y, xa, ya, zxval, i);

      test_single_low_level(&gsl_interp2d_eval_deriv_y, &gsl_interp2d_eval_deriv_y_e,  interp, xarr, yarr, zarr, x, y, xa, ya, zyval, i);

      test_single_low_level(&gsl_interp2d_eval_deriv_xx,&gsl_interp2d_eval_deriv_xx_e, interp, xarr, yarr, zarr, x, y, xa, ya, zxxval, i);

      test_single_low_level(&gsl_interp2d_eval_deriv_yy,&gsl_interp2d_eval_deriv_yy_e, interp, xarr, yarr, zarr, x, y, xa, ya, zyyval, i);

      test_single_low_level(&gsl_interp2d_eval_deriv_xy,&gsl_interp2d_eval_deriv_xy_e, interp, xarr, yarr, zarr, x, y, xa, ya, zxyval, i);

      test_single_high_level(&gsl_spline2d_eval,         &gsl_spline2d_eval_e,          interp_s, x, y, xa, ya, zval, i);

      test_single_high_level(&gsl_spline2d_eval_deriv_x, &gsl_spline2d_eval_deriv_x_e,  interp_s, x, y, xa, ya, zxval, i);

      test_single_high_level(&gsl_spline2d_eval_deriv_y, &gsl_spline2d_eval_deriv_y_e,  interp_s, x, y, xa, ya, zyval, i);

      test_single_high_level(&gsl_spline2d_eval_deriv_xx,&gsl_spline2d_eval_deriv_xx_e, interp_s, x, y, xa, ya, zxxval, i);

      test_single_high_level(&gsl_spline2d_eval_deriv_yy,&gsl_spline2d_eval_deriv_yy_e, interp_s, x, y, xa, ya, zyyval, i);

      test_single_high_level(&gsl_spline2d_eval_deriv_xy,&gsl_spline2d_eval_deriv_xy_e, interp_s, x, y, xa, ya, zxyval, i);

      test_single_low_level(&gsl_interp2d_eval_extrap, &gsl_interp2d_eval_e_extrap, interp, xarr, yarr, zarr, x, y, xa, ya, zval, i);

    }

  gsl_interp_accel_free(xa);

  gsl_interp_accel_free(ya);

  gsl_interp2d_free(interp);

  gsl_spline2d_free(interp_s);

  return status;

}

/*

 * Tests bilinear interpolation using a symmetric function, f(x,y)==f(y,x),

 * and diagonal interpolation points (x,y) where x==y. If these tests don't

 * pass, something is seriously broken.

 */

static int

test_bilinear_symmetric()

{

  int status;

  double xarr[] = {0.0, 1.0, 2.0, 3.0};

  double yarr[] = {0.0, 1.0, 2.0, 3.0};

  double zarr[] = {1.0, 1.1, 1.2, 1.3,

                   1.1, 1.2, 1.3, 1.4,

                   1.2, 1.3, 1.4, 1.5,

                   1.3, 1.4, 1.5, 1.6};

  double xval[] = {0.0, 0.5, 1.0, 1.5, 2.5, 3.0};

  double yval[] = {0.0, 0.5, 1.0, 1.5, 2.5, 3.0};

  double zval[] = {1.0, 1.1, 1.2, 1.3, 1.5, 1.6};

  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);

  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);

  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,

                         NULL, NULL, NULL, NULL, NULL, test_size,

                         gsl_interp2d_bilinear);

  gsl_test(status, "bilinear interpolation with symmetric values");

  return status;

}

/*

 * Tests bilinear interpolation with an asymmetric function, f(x,y)!=f(y,x),

 * and off-diagonal interpolation points (x,y) where x and y may or may not

 * be equal.

 */

static int

test_bilinear_asymmetric_z()

{

  int status;

  double xarr[] = {0.0, 1.0, 2.0, 3.0};

  double yarr[] = {0.0, 1.0, 2.0, 3.0};

  double zarr[] = {1.0, 1.1, 1.2, 1.4,

                   1.3, 1.4, 1.5, 1.7,

                   1.5, 1.6, 1.7, 1.9,

                   1.6, 1.9, 2.2, 2.3};

  double xval[] = { 0.0, 0.5, 1.0, 1.5,  2.5, 3.0,

                    1.3954, 1.6476, 0.824957,

                    2.41108,  2.98619, 1.36485 };

  double yval[] = {0.0, 0.5, 1.0, 1.5,  2.5, 3.0,

                   0.265371, 2.13849, 1.62114,

                   1.22198, 0.724681, 0.0596087 };

  /* results computed using Mathematica 9.0.1.0 */

  double zval[] = {1.0, 1.2, 1.4, 1.55, 2.025, 2.3,

                   1.2191513, 1.7242442248, 1.5067237,

                   1.626612, 1.6146423, 1.15436761};

  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);

  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);

  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,

                         NULL, NULL, NULL, NULL, NULL, test_size,

                         gsl_interp2d_bilinear);

  gsl_test(status, "bilinear interpolation with asymmetric z values");

  return status;

}

static int

test_bicubic()

{

  int status;

  double xarr[] = {0.0, 1.0, 2.0, 3.0};

  double yarr[] = {0.0, 1.0, 2.0, 3.0};

  double zarr[] = {1.0, 1.1, 1.2, 1.3,

                   1.1, 1.2, 1.3, 1.4,

                   1.2, 1.3, 1.4, 1.5,

                   1.3, 1.4, 1.5, 1.6};

  double xval[] = {1.0, 1.5, 2.0};

  double yval[] = {1.0, 1.5, 2.0};

  double zval[] = {1.2, 1.3, 1.4};

  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);

  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);

  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,

                         NULL, NULL, NULL, NULL, NULL, test_size,

                         gsl_interp2d_bicubic);

  gsl_test(status, "bicubic interpolation on linear function");

  return status;

}

static int

test_bicubic_nonlinear()

{

  int status;

  double xarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};

  double yarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};

  /* least common multiple of x and y */

  double zarr[] = { 1,  2,  3,  4,  5,  6,  7,  8,

                    2,  2,  6,  4, 10,  6, 14,  8,

                    3,  6,  3, 12, 15,  6, 21, 24,

                    4,  4, 12,  4, 20, 12, 28,  8,

                    5, 10, 15, 20,  5, 30, 35, 40,

                    6,  6,  6, 12, 30,  6, 42, 24,

                    7, 14, 21, 28, 35, 42,  7, 56,

                    8,  8, 24,  8, 40, 24, 56,  8};

  double xval[] = {1.4, 2.3, 4.7, 3.3, 7.5, 6.6, 5.1};

  double yval[] = {1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3};

  /* results computed using GSL 1D cubic interpolation twice */

  double zval[] = { 1.4, 3.11183531264736, 8.27114315792559, 5.03218982537718,

                    22.13230634702637, 23.63206834997871, 17.28553080971182 };

  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);

  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);

  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,

                         NULL, NULL, NULL, NULL, NULL, test_size,

                         gsl_interp2d_bicubic);

  gsl_test(status, "bicubic interpolation on nonlinear symmetric function");

  return status;

}

/* This function contributed by Andrew W. Steiner <awsteiner@gmail.com> */

static int

test_bicubic_nonlinear_nonsq()

{

  int status;

  double xarr[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};

  double yarr[] = {1.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0};

  double zarr[] = { 1,  2,  3,  4,  5,  6,  7,  8, 9, 10,

                    2,  2,  6,  4, 10,  6, 14,  8, 11, 12,

                    3,  6,  3, 12, 15,  6, 21, 24, 13, 14,

                    4,  4, 12,  4, 20, 12, 28,  8, 15, 16,

                    5, 10, 15, 20,  5, 30, 35, 40, 17, 18,

                    6,  6,  6, 12, 30,  6, 42, 24, 19, 20,

                    7, 14, 21, 28, 35, 42,  7, 56, 21, 22,

                    8,  8, 24,  8, 40, 24, 56,  8, 23, 24};

  double xval[] = {1.4, 2.3, 9.7, 3.3, 9.5, 6.6, 5.1};

  double yval[] = {1.0, 1.8, 1.9, 2.5, 2.7, 4.1, 3.3};

  /* results computed using GSL 1D cubic interpolation twice */

  double zval[] = { 1.4, 2.46782030941187003, 10.7717721621846465,

                    4.80725067958096375, 11.6747032398627297,

                    11.2619968682970111, 9.00168877916872567};

  size_t xsize = sizeof(xarr) / sizeof(xarr[0]);

  size_t ysize = sizeof(yarr) / sizeof(yarr[0]);

  size_t test_size = sizeof(xval) / sizeof(xval[0]);

  status = test_interp2d(xarr, yarr, zarr, xsize, ysize, xval, yval, zval,

                         NULL, NULL, NULL, NULL, NULL, test_size,

                         gsl_interp2d_bicubic);

  gsl_test(status, "bicubic interpolation on nonlinear symmetric function");

  return status;

}

/* runs all the tests */

int

test_interp2d_main(void)

{

  int status = 0;

  status += test_bilinear_symmetric();

  status += test_bilinear_asymmetric_z();

  status += test_bicubic();

  status += test_bicubic_nonlinear();

  status += test_bicubic_nonlinear_nonsq();

  return status;

}