在 C++ 中使用 GSL Monte Carlo 积分进行多维函数

Question

我正在尝试在我生成的 C++ 代码中使用 GSL Monte Carlo 集成。 这个想法是有一个布朗运动函数 (brownian)，它用于另一个函数 (g) 以执行 4-dim 数值积分。

double brownian(const double &x,const double &x0,const double & sigma,const double & t) {
   double a1 = (1/sqrt(2.0 * M_PI * sigma * t));
   double b1 = exp(-((x - x0) * (x - x0))/(2.0 * sigma * t));
   double res = a1 * b1;
  return res;
}

double g(double *k, size_t dim, void *p[]){

 const double& xA0 = *(double *)p[0];
 const double& xB0 = *(double *)p[1];
 const double& yA0 = *(double *)p[2];
 const double& yB0 = *(double *)p[3];
 const double& sigma = *(double *)p[4];
 const double& t1 = *(double *)p[5];


  double temp_pbxA = brownian(k[0], xA0, sigma,t1);
  double temp_pbxB = brownian(k[1], xB0, sigma, t1);
  double temp_pbyA = brownian(k[2], yA0, sigma,t1);
  double temp_pbyB = brownian(k[3], yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);
}

double integrate_test(const double xl, const double xu, const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){

  double res, err;

 const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G;
  G.f = &g;
  G.dim = 4;
  G.params = { xA0, xB0, yA0, yB0, sigma, t1};

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);


    gsl_monte_plain_state *s = gsl_monte_plain_alloc (4);
    gsl_monte_plain_integrate (&G, xl, xu, 4, calls, r, s,
                               &res, &err);
    gsl_monte_plain_free (s);

return res;
    }

但是，当我尝试编译代码时，出现以下错误：

error: assigning to 'double (*)(double *, size_t, void *)' (aka 'double (*)(double *, unsigned long, void *)') from incompatible type 'double (*)(double *, size_t, void **)' (aka 'double (*)(double *, unsigned long, void **)'): type mismatch at 3rd parameter ('void *' vs 'void **')
  G.f = &g;

我不明白为什么它使用void *p[]，因为它是void**。 有什么建议吗？

Answer 1

这是您的代码与可编译的代码之间的差异。这并不意味着我的代码是正确的：你只会得到你所要求的——因为你没有提供一个完整的最低限度的工作示例，任何试图帮助你的人都不得不猜测你的问题到底是什么。

13c13
< double g(double *k, size_t dim, void *p){
---
> double g(double *k, size_t dim, void *p[]){
15,20c15,20
<  const double& xA0 = ((double *)p)[0];
<  const double& xB0 = ((double *)p)[1];
<  const double& yA0 = ((double *)p)[2];
<  const double& yB0 = ((double *)p)[3];
<  const double& sigma = ((double *)p)[4];
<  const double& t1 = ((double *)p)[5];
---
>  const double& xA0 = *(double *)p[0];
>  const double& xB0 = *(double *)p[1];
>  const double& yA0 = *(double *)p[2];
>  const double& yB0 = *(double *)p[3];
>  const double& sigma = *(double *)p[4];
>  const double& t1 = *(double *)p[5];
31c31
< double integrate_test(const double xl[], const double xu[], const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){
---
> double integrate_test(const double xl, const double xu, const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){
41,42c41
<   double params[] = { xA0, xB0, yA0, yB0, sigma, t1};
<   G.params = params;
---
>   G.params = { xA0, xB0, yA0, yB0, sigma, t1};

如您所见，您的编译器至少抱怨了两个问题。

1. gsl_monte_function 期望组件 f 具有签名 double (* f) (double * x, size_t dim, void * params)，请参阅 https://www.gnu.org/software/gsl/doc/html/montecarlo.html 。我必须强调，通常 GSL 中的参数是通过 void 指针传递给必须事先定义的结构的，但在这里你很幸运，所有参数都只是doubles，所以为了节省额外的代码，你可以使用旧的，普通 C 数组

double params[] = { xA0, xB0, yA0, yB0, sigma, t1};

并将其地址作为G.params 传递。然后你将使用像这样的 C 语法在 g 中检索它们

const double& xA0 = ((double *)p)[0];

（将 void 指针转换为 double*，然后将其视为 C 样式数组的名称）。

1. gsl_monte_plain_integrate 期望它的第二个和第三个参数是 C 风格的数组。见https://www.gnu.org/software/gsl/doc/html/montecarlo.html。这就是为什么您需要在integrate_test 中更改xl 和xu 的类型。

GSL 很难使用。如有疑问，请随时询问。