Интегральная оптимизация Монте-Карло - PullRequest
Новая
       13

Интегральная оптимизация Монте-Карло

0 голосов
/ 01 ноября 2019

Я использую метод Монте-Карло, реализованный в библиотеке gsl. Мне нужно вычислить много повторений этого интеграла, меняя параметр в подынтегральном выражении. Поэтому мне нужно сделать мою подпрограмму быстрой. Кажется, что самая трудоемкая часть - это оценка подынтегральной функции в случайных точках. Как я могу сделать оценку быстрее в моем конкретном случае? Вот минимальный пример: 

#include <gsl/gsl_rng.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>

double q=0.0;
double mu=0.001;
double eta=0.1;
double kF=1.0;
double Kcut=10;
long int Nmax=10000000;
int Nwu=1000000;
double w=1;


struct my_f_params { double y;};

double 
g (double *k, size_t dim, void *p)
{
  double A;
  struct my_f_params * fp = (struct my_f_params *)p;

  double PQ=q*q+k[1]*k[1]-2*q*k[1]*cos(k[3])+mu;
  double QK=k[0]*k[0]+k[1]*k[1]-2*k[0]*k[1]* (cos(k[2])*cos(k[3])+cos(k[4])*sin(k[2])*sin(k[3]))+mu;
  double KPQ=q*q+k[0]*k[0]+k[1]*k[1]+2*k[0]*cos(k[2])*(q-k[1]*cos(k[3]))+2*k[1]*       (q*cos(k[3])+k[0]*cos(k[4])*sin(k[2])*sin(k[3]));
  double denFreq=fp->y-0.5*(k[0]*k[0]+k[1]*k[1]+KPQ);
  double vol=k[0]*k[0]*k[1]*k[1]*sin(k[2])*sin(k[3]);
  if (sqrt(KPQ) < kF) {
    A = vol*denFreq*(1/QK-1/PQ)/(QK*(pow(denFreq,2)+eta*eta));
  }
  else {
    A = 0;
  }

  return A; 
}


int
main (void)
{
  double res, err;
  double xl[5] = {0, kF, 0, 0, 0};
  double xu[5] = {kF, Kcut, M_PI, M_PI, 2*M_PI};
  const gsl_rng_type *T;
  gsl_rng *r;
  gsl_monte_function G;
  size_t calls = Nmax;

  gsl_rng_env_setup ();
  struct my_f_params params;
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  params.y=w;
  G.f=&g;
  G.dim=5;
  G.params=&params;

  {
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (5);
    gsl_monte_vegas_integrate (&G, xl, xu, 5, Nwu, r, s,&res, &err);
    do
      {
        gsl_monte_vegas_integrate (&G, xl, xu, 5, calls/5, r, s,&res, &err);
      }
    while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);
    gsl_monte_vegas_free (s);
  }

  printf ("%.6f   %.6f   %.6f\n", w,res,err);
  gsl_rng_free (r);



  return 0;
}

1 Ответ

1 голос
/ 03 ноября 2019

Расширяя комментарий Боба, вы можете использовать sincos для вычисления sin и cos того же аргумента (k[2] и k[3]) и определить kF_sqr для инициализации в главноми используйте это в функции g, чтобы избежать вызова sqrt. С этими оптимизациями быстрый и грязный тест на моей машине показал ускорение на 5% по сравнению с вашим кодом.

...