Попытка интегрировать интерполированную функцию с C GSL

Question

Я новичок в языке C, и у меня возникают проблемы с объявлением типа и указателями.

Я пытаюсь интерполировать набор данных, а затем интегрировать его, и численные методы выполняются с использованием библиотеки C GSL. Код адаптирован из примеров руководства GSL по интеграции monte carlo и интерполяции 2d. Прежде чем использовать набор данных, я делаю тест, используя известную и очень простую 2D-функцию.

2D-интерполяционная часть работает хорошо, и я могу получить значения от интерполированной функции, interpFun.

Но я не могу использовать эту интерполированную функцию в качестве входного сигнала для алгоритма интеграции.

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

Этот код возвращает сообщение об ошибке:

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] 

Если я пытаюсь интегрировать fun вместо interpFun, нет никакого сообщения об ошибке, но и к тому же типу.

Answer 1

interpFun() функция, как это определено в программе имеет тип double (*)(double *, size_t, struct fparams), но по the GSL docs вам нужен тип, чтобы быть double (*)(double *, size_t, void *). Вот исправление:

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

Кроме того, C не допускает определения вложенных функций, хотя есть GCC extension for this.