#include <iostream>
#include <iomanip>
using namespace std;

#include "MpIeee.hh"
#include "CMpIeee.hh"
#include "ArithmosIO.hh"

/* interpolation/test.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */


/* Author:  G. Jungman
 */
#include <config.h>
#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_ieee_utils.h>

int
 test_bsearch(void)
{
  MpIeee x_array[5] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ), MpIeee( "4.0" ) };
  size_t index_result;
  int  status=  0;
  int  s;

  /* check an interior point */
  index_result = gsl_interp_bsearch(x_array, 1.5, 0, 4);
  s = (index_result != 1);
  status += s;
  gsl_test (s, "simple bsearch");

  /* check that we get the last interval if x == last value */
  index_result = gsl_interp_bsearch(x_array, 4.0, 0, 4);
  s = (index_result != 3);
  status += s;
  gsl_test (s, "upper endpoint bsearch");

  /* check that we get the first interval if x == first value */
  index_result = gsl_interp_bsearch(x_array, 0.0, 0, 4);
  s = (index_result != 0);
  status += s;
  gsl_test (s, "lower endpoint bsearch");

  /* check that we get correct interior boundary behaviour */
  index_result = gsl_interp_bsearch(x_array, 2.0, 0, 4);
  s = (index_result != 2);
  status += s;
  gsl_test (s, "degenerate bsearch");

  /* check out of bounds above */
  index_result = gsl_interp_bsearch(x_array, 10.0, 0, 4);
  s = (index_result != 3);
  status += s;
  gsl_test (s, "out of bounds bsearch +");

  /* check out of bounds below */
  index_result = gsl_interp_bsearch(x_array, -10.0, 0, 4);
  s = (index_result != 0);
  status += s;
  gsl_test (s, "out of bounds bsearch -");

  return status;
}



typedef MpIeee TEST_FUNC(MpIeee);
typedef struct _xy_table xy_table;

struct _xy_table
  {
    MpIeee * x;
    MpIeee * y;
    size_t n;
  };

xy_table make_xy_table (MpIeee x[], MpIeee y[], size_t n);

xy_table
make_xy_table (MpIeee x[], MpIeee y[], size_t n)
{
  xy_table t;
  t.x = x;
  t.y = y;
  t.n = n;
  return t;
}

static int
 test_interp(
  const xy_table * data_table,
  const gsl_interp_type * T,
  xy_table * test_table,
  xy_table * test_d_table,
  xy_table * test_i_table
  )
{
  int  status=  0;
  size_t i;

  gsl_interp_accel *a = gsl_interp_accel_alloc ();
  gsl_interp *interp = gsl_interp_alloc (T, data_table->n);

  gsl_interp_init (interp, data_table->x, data_table->y, data_table->n);

  for (i = 0; i < test_table->n; i++)
    {
      MpIeee x=  test_table->x[i];
      MpIeee y;
      MpIeee deriv;
      MpIeee integ;
      MpIeee diff_y;MpIeee  diff_deriv;MpIeee  diff_integ;
      gsl_interp_eval_e (interp, data_table->x, data_table->y, x, a, &y);
      gsl_interp_eval_deriv_e (interp, data_table->x, data_table->y, x, a, &deriv);
      gsl_interp_eval_integ_e (interp, data_table->x, data_table->y, 0.0, x, a, &integ);
      diff_y = y - test_table->y[i];
      diff_deriv = deriv - test_d_table->y[i];
      diff_integ = integ - test_i_table->y[i];
      if (fabs (diff_y) > 1.e-10 || fabs(diff_deriv) > 1.0e-10 || fabs(diff_integ) > 1.0e-10) {
        status++;
      }
    }

  gsl_interp_accel_free (a);
  gsl_interp_free (interp);

  return status;
}

static int
 test_linear(void)
{
  int  s;

  MpIeee data_x[4] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ) };
  MpIeee data_y[4] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ) };
  MpIeee test_x[6] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "1.5" ), MpIeee( "2.5" ), MpIeee( "3.0" ) };
  MpIeee test_y[6] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "1.5" ), MpIeee( "2.5" ), MpIeee( "3.0" ) };
  MpIeee test_dy[6] =  { MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ) };
  MpIeee test_iy[6] =  { MpIeee( "0.0" ), MpIeee( "0.125" ), MpIeee( "0.5" ), MpIeee( "9.0" )/MpIeee( "8.0" ), MpIeee( "25.0" )/MpIeee( "8.0" ), MpIeee( "9.0" )/MpIeee( "2.0" ) };

  xy_table data_table = make_xy_table(data_x, data_y, 4);
  xy_table test_table = make_xy_table(test_x, test_y, 6);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 6);

  s = test_interp (&data_table, gsl_interp_linear, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "linear interpolation");
  return s;
}

static int
 test_polynomial(void)
{
  int  s;

  MpIeee data_x[4] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ) };
  MpIeee data_y[4] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ) };
  MpIeee test_x[6] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "1.5" ), MpIeee( "2.5" ), MpIeee( "3.0" ) };
  MpIeee test_y[6] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "1.5" ), MpIeee( "2.5" ), MpIeee( "3.0" ) };
  MpIeee test_dy[6] =  { MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ) };
  MpIeee test_iy[6] =  { MpIeee( "0.0" ), MpIeee( "0.125" ), MpIeee( "0.5" ), MpIeee( "9.0" )/MpIeee( "8.0" ), MpIeee( "25.0" )/MpIeee( "8.0" ), MpIeee( "9.0" )/MpIeee( "2.0" ) };

  xy_table data_table = make_xy_table(data_x, data_y, 4);
  xy_table test_table = make_xy_table(test_x, test_y, 6);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 6);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 6);

  s = test_interp (&data_table, gsl_interp_polynomial, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "polynomial interpolation");
  return s;
}


static int
 test_cspline(void)
{
  int  s;

  MpIeee data_x[3] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee data_y[3] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee test_x[4] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee test_y[4] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee test_dy[4] =  { MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ) };
  MpIeee test_iy[4] =  { MpIeee( "0.0" ), MpIeee( "0.125" ), MpIeee( "0.5" ), MpIeee( "2.0" ) };

  xy_table data_table = make_xy_table(data_x, data_y, 3);
  xy_table test_table = make_xy_table(test_x, test_y, 4);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 4);

  s = test_interp (&data_table, gsl_interp_cspline, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "cspline interpolation");
  return s;
}


static int
 test_akima(void)
{
  int  s;

  MpIeee data_x[5] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ), MpIeee( "4.0" ) };
  MpIeee data_y[5] =  { MpIeee( "0.0" ), MpIeee( "1.0" ), MpIeee( "2.0" ), MpIeee( "3.0" ), MpIeee( "4.0" ) };
  MpIeee test_x[4] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee test_y[4] =  { MpIeee( "0.0" ), MpIeee( "0.5" ), MpIeee( "1.0" ), MpIeee( "2.0" ) };
  MpIeee test_dy[4] =  { MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ), MpIeee( "1.0" ) };
  MpIeee test_iy[4] =  { MpIeee( "0.0" ), MpIeee( "0.125" ), MpIeee( "0.5" ), MpIeee( "2.0" ) };

  xy_table data_table = make_xy_table(data_x, data_y, 5);
  xy_table test_table = make_xy_table(test_x, test_y, 4);
  xy_table test_d_table = make_xy_table(test_x, test_dy, 4);
  xy_table test_i_table = make_xy_table(test_x, test_iy, 4);

  s = test_interp (&data_table, gsl_interp_akima, &test_table, &test_d_table, &test_i_table);
  gsl_test (s, "akima interpolation");
  return s;
}


int 
main (int argc, char **argv)
{
  int  status=  0;

  gsl_ieee_env_setup ();

  argc = 0;    /* prevent warnings about unused parameters */
  argv = 0;

  status += test_bsearch();
  status += test_linear();
  status += test_polynomial();
  status += test_cspline();
  status += test_akima();

  exit (gsl_test_summary());
}
