Search code examples
c++functionphysicsintegral

workaround to local functions in c++ for numerical integration


For a physics related calculation, I need to evaluate numerically a four dimensional integral that depends on several parameters that need to be varied; such that they can't be globally defined. I am using the Cuba package which provides the following function:

Cuhre(NDIM, NCOMP, Integrand , USERDATA, NVEC,EPSREL, EPSABS, VERBOSE | LAST,MINEVAL,
 MAXEVAL, KEY,STATEFILE, SPIN,&nregions, &neval, &fail, integral, error, prob);

The argument "Integrand" is a function that is expected to be defined in the following fashion:

int Integrand(const int *ndim, const cubareal xx[],const int *ncomp,
cubareal ff[], void *userdata) ;

My problem is that I would like to integrate a function that depends on continuous parameters that I would like to vary within the program. What I mean by this is that my integrand function depends on several extra parameters, say a,b,c.

I would then like to define something like the Integrand function locally, i.e. in a scope where the parameters a, b and c are fixed, so that it will have the correct form to be passed as an argument to the Cuhre function.

I then thought that I could perhaps put everything into a class such as:

 class integrate{
  private: // parameters

  public:
//constructor
int Integrand(...){
  // calculation involving the parameters 
}

  };

and define the function as a method for that class. This however didn't work because I can't pass non-static object methods to a function, and, making the method static would also force me to make the parameters static which means I would not be able to initialise them within the program.

Is there any workaround to this?


Solution

  • [Edit: as you misunderstood the userdata pointer, as I see from the comments, let me show you how to use it:]

    struct Parameters
    {
        int a, b, c; // TODO: appropriate defaults
    };
    
    int integrand
    (
        const int* ndim, const cubareal xx[],
        const int* ncomp, cubareal ff[],
        void* userdata
    )
    {
        Parameters* p = reinterpret_cast<Parameters*>(userdata);
        /* calculations using p->a, p->b, p->c */
    }
    
    int main(int argc, char* argv[])
    {
        Parameters p;
        /* fill in parameters appropriately, e. g. parsing command line parameters */
        Cuhre(/* other parameters */, &integrand, &p, /* remaining parameters */);
    }
    

    My previous solution remains as an alternative...

    .cpp file:

    namespace
    {
        int g_a = 0; // appropriate default values...
        int g_b = 0;
        int g_c = 0;
    };
    
    void setIntegrationParameters (int a, int b, int c)
    {
        g_a = a;
        g_b = b;
        g_c = c;
    }
    
    int integrand
    (
        const int* ndim, const cubareal xx[],
        const int* ncomp, cubareal ff[],
        void* userdata
    )
    {
        /* calculations using g_a, g_b, g_c */
    }
    

    .h file:

    void setIntegrationParameters (int a, int b, int c);
    int integrand
    (
        const int* ndim, const cubareal xx[],
        const int* ncomp, cubareal ff[],
        void* userdata
    );
    

    And from within your main program, you can e. g. parse command line parameters to calculate a, b, c and set them as parameters for the integration via setIntegrationParameters.

    (You could do exactly the same with C code, of course then using C-Casts and static global variables instead of anonymous namespace...)