Trying to integrate an interpolated function with C GSL

Question

I am new to the C language and I'm having some trouble with type declaration and pointers.

I'm trying to interpolate a dataset then integrate it, both numerical methods are done using the C GSL library. The code is adapted from GSL's manual examples on monte carlo integration and 2d interpolation. Before using a dataset, I'm doing a test using a known and very simple 2D function.

The 2D interpolation part is working well, and I can obtain values from the interpolated function, interpFun.

But I cannot use this interpolated function as input to the integration algorithm.

#include <stdio.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>

int
main()
{
  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
  const size_t N = 100;             /* number of points to interpolate */
  double xa[] = { 0.0, 1.0 }; /* define unit square */
  double ya[] = { 0.0, 1.0 };
  double xl[] = { 0., 0.};
  double xu[] = { 1., 1.}, xp[] = {1.0, 0.0}, zz[] = {.2,.6};
  const size_t nx = sizeof(xa) / sizeof(double); /* x grid points */
  const size_t ny = sizeof(ya) / sizeof(double); /* y grid points */
  double *za = malloc(nx * ny * sizeof(double));
  double fun( double *x, size_t dim, void *params);
  gsl_spline2d *spline = gsl_spline2d_alloc(T, nx, ny);
  gsl_interp_accel *xacc = gsl_interp_accel_alloc();
  gsl_interp_accel *yacc = gsl_interp_accel_alloc();
  size_t i, j;
  const gsl_rng_type *S; 
  gsl_rng *r; 
  double res, err;
  struct fparams { gsl_spline2d *spline; gsl_interp_accel *xacc; gsl_interp_accel *yacc; };

  /* set z grid values */
  gsl_spline2d_set(spline, za, 0, 0, fun( xl, 2, 0) );
  gsl_spline2d_set(spline, za, 0, 1, fun( xa, 2, 0) );
  gsl_spline2d_set(spline, za, 1, 1, fun( xu, 2, 0) );
  gsl_spline2d_set(spline, za, 1, 0, fun( xp, 2, 0) );

  /* initialize interpolation */
  gsl_spline2d_init(spline, xa, ya, za, nx, ny);

  struct fparams params1 = { spline, xacc, yacc };

  double interpFun( double *x, size_t dim, struct fparams p ) {
  (void)(dim);
  struct fparams *par = &p;
  return gsl_spline2d_eval(par->spline, x[0], x[1], par->xacc, par->yacc);
   };

  gsl_monte_function G = { &interpFun, 2, &params1 };
  size_t calls = 500000;
  gsl_rng_env_setup ();
  S = gsl_rng_default;
  r = gsl_rng_alloc (S);
  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (2);
    gsl_monte_vegas_integrate (&G, xl, xu, 2, 10000, r, s, &res, &err);
    printf ("result = %.6f\n", res);
    printf ("sigma = %.6f\n", err);
    printf ("converging...\n");
    do
      {
       gsl_monte_vegas_integrate (&G, xl, xu, 2, calls/5, r, s, &res, &err);
       printf ("result = %.6f sigma = %.6f "
               "chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s));
     }
    while (fabs (gsl_monte_vegas_chisq(s)-1.0) > 0.5);
    printf ("result = %.6f\n", res);
    printf ("sigma = %.6f\n", err);
    gsl_monte_vegas_free (s);
  }

  gsl_rng_free (r);
  gsl_spline2d_free(spline);
  gsl_interp_accel_free(xacc);
  gsl_interp_accel_free(yacc);
  free(za);

  return 0;
}

double fun( double *x, size_t dim, void *p ) {
 (void)(dim);
 (void)(p);
 return 3*x[0]+5*x[1];
}

This code returns the error message:

inter_test.c: In function ‘main’:
inter_test.c:55:3: warning: initialization from incompatible pointer type [enabled by default]
   gsl_monte_function G = { &interpFun, 2, &params1 };
   ^
inter_test.c:55:3: warning: (near initialization for ‘G.f’) [enabled by default]

If I try integrating fun instead of interpFun, there is no error message, but both are the same type.

Answer 1

The interpFun() function as defined in your program has type double (*)(double *, size_t, struct fparams), but according to the GSL docs you need the type to be double (*)(double *, size_t, void *). Here is a fix:

double interpFun( double *x, size_t dim, void *p ) {
    (void)(dim);
    struct fparams *par = (struct fparams *) p;
    return gsl_spline2d_eval(par->spline, x[0], x[1], par->xacc, par->yacc);
}

Additionally, C does not allow nested function definitions, though there is a GCC extension for this.