integrate.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2015-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#ifndef SRC_INCLUDE_INTEGRATE_H_
#define SRC_INCLUDE_INTEGRATE_H_

#include <cuba.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>

#include <algorithm>
#include <iostream>
#include <memory>
#include <sstream>
#include <string>
#include <tuple>
#include <utility>

#include "cxx14compat.h"
#include "fpenvironment.h"
#include "random.h"

namespace smash {

struct GslWorkspaceDeleter {
  constexpr GslWorkspaceDeleter() = default;

  void operator()(gsl_integration_cquad_workspace *ptr) const {
    if (ptr == nullptr) {
      return;
    }
    gsl_integration_cquad_workspace_free(ptr);
  }
};

class Result : public std::pair<double, double> {
  using Base = std::pair<double, double>;

 public:
  using Base::pair;

  operator double() const { return Base::first; }

  double value() const { return Base::first; }

  double error() const { return Base::second; }

  void check_error(const std::string &integration_name,
                   double relative_tolerance = 5e-4,
                   double absolute_tolerance = 1e-9) const {
    const double allowed_error =
        std::max(absolute_tolerance, value() * relative_tolerance);
    if (error() > allowed_error) {
      std::stringstream error_msg;
      error_msg << integration_name << " resulted in I = " << value() << " ± "
                << error()
                << ", but the required precision is either absolute error < "
                << absolute_tolerance << " or relative error < "
                << relative_tolerance << std::endl;
      throw std::runtime_error(error_msg.str());
    }
  }
};

class Integrator {
 public:
  explicit Integrator(size_t workspace_size = 1000)
      : workspace_(gsl_integration_cquad_workspace_alloc(workspace_size)) {}

  template <typename F>
  Result operator()(double a, double b, F &&fun) {
    Result result;
    const gsl_function gslfun{
        /* important! The lambda cannot use captures, otherwise the
         * conversion to a function pointer type is impossible. */
        [](double x, void *type_erased) -> double {
          auto &&f = *static_cast<F *>(type_erased);
          return f(x);
        },
        &fun};
    // We disable float traps when calling GSL code we cannot control.
    DisableFloatTraps guard;
    const int error_code = gsl_integration_cquad(
        &gslfun, a, b, accuracy_absolute_, accuracy_relative_, workspace_.get(),
        &result.first, &result.second,
        nullptr /* Don't store the number of evaluations */);
    if (error_code) {
      std::stringstream err;
      err << "GSL 1D deterministic integration: " << gsl_strerror(error_code);
      throw std::runtime_error(err.str());
    }
    result.check_error("GSL 1D deterministic integration", accuracy_relative_,
                       accuracy_absolute_);
    return result;
  }

 private:
  std::unique_ptr<gsl_integration_cquad_workspace, GslWorkspaceDeleter>
      workspace_;

  const double accuracy_absolute_ = 1.0e-9;

  const double accuracy_relative_ = 5.0e-4;
};

class Integrator1dMonte {
 public:
  explicit Integrator1dMonte(size_t num_calls = 1E6)
      : state_(gsl_monte_plain_alloc(1)),
        rng_(gsl_rng_alloc(gsl_rng_mt19937)),
        number_of_calls_(num_calls) {
    gsl_monte_plain_init(state_);
    // initialize the GSL RNG with a random seed
    const uint32_t seed = random::uniform_int(0ul, ULONG_MAX);
    gsl_rng_set(rng_, seed);
  }

  ~Integrator1dMonte() {
    gsl_monte_plain_free(state_);
    gsl_rng_free(rng_);
  }

  template <typename F>
  Result operator()(double min, double max, F &&fun) {
    Result result = {0, 0};

    const double lower[1] = {min};
    const double upper[1] = {max};

    if (max <= min)
      return result;

    const gsl_monte_function monte_fun{
        // trick: pass integrand function as 'params'
        [](double *x, size_t /*dim*/, void *params) -> double {
          auto &&f = *static_cast<F *>(params);
          return f(x[0]);
        },
        1, &fun};

    const int error_code =
        gsl_monte_plain_integrate(&monte_fun, lower, upper, 1, number_of_calls_,
                                  rng_, state_, &result.first, &result.second);
    if (error_code) {
      std::stringstream err;
      err << "GSL 1D Monte-Carlo integration: " << gsl_strerror(error_code);
      throw std::runtime_error(err.str());
    }

    result.check_error("GSL 1D Monte-Carlo integration");

    return result;
  }

 private:
  gsl_monte_plain_state *state_;

  gsl_rng *rng_;

  const std::size_t number_of_calls_;
};

class Integrator2d {
 public:
  explicit Integrator2d(size_t num_calls = 1E6)
      : state_(gsl_monte_plain_alloc(2)),
        rng_(gsl_rng_alloc(gsl_rng_mt19937)),
        number_of_calls_(num_calls) {
    gsl_monte_plain_init(state_);
    // initialize the GSL RNG with a random seed
    const uint32_t seed = random::uniform_int(0ul, ULONG_MAX);
    gsl_rng_set(rng_, seed);
  }

  ~Integrator2d() {
    gsl_monte_plain_free(state_);
    gsl_rng_free(rng_);
  }

  template <typename F>
  Result operator()(double min1, double max1, double min2, double max2,
                    F &&fun) {
    Result result = {0, 0};

    const double lower[2] = {min1, min2};
    const double upper[2] = {max1, max2};

    if (max1 <= min1 || max2 <= min2)
      return result;

    const gsl_monte_function monte_fun{
        // trick: pass integrand function as 'params'
        [](double *x, size_t /*dim*/, void *params) -> double {
          auto &&f = *static_cast<F *>(params);
          return f(x[0], x[1]);
        },
        2, &fun};

    gsl_monte_plain_integrate(&monte_fun, lower, upper, 2, number_of_calls_,
                              rng_, state_, &result.first, &result.second);

    return result;
  }

 private:
  gsl_monte_plain_state *state_;

  gsl_rng *rng_;

  const std::size_t number_of_calls_;
};

template <typename F>
struct Integrand2d {
  double min1;
  double diff1;
  double min2;
  double diff2;
  F f;
};

class Integrator2dCuhre {
 public:
  explicit Integrator2dCuhre(int num_calls = 1e6, double epsrel = 5e-4,
                             double epsabs = 1e-9)
      : maxeval_(num_calls), epsrel_(epsrel), epsabs_(epsabs) {}

  template <typename F>
  Result operator()(double min1, double max1, double min2, double max2, F fun) {
    Result result = {0., 0.};

    if (max1 < min1 || max2 < min2) {
      bool tolerable = (max1 - min1 > -1.e-16) && (max2 - min2 > -1.e-16);
      if (tolerable) {
        return result;
      }
      std::stringstream err;
      err << "Integrator2dCuhre got wrong integration limits: [" << min1 << ", "
          << max1 << "], [" << min2 << ", " << max2 << "]";
      throw std::invalid_argument(err.str());
    }

    Integrand2d<F> f_with_limits = {min1, max1 - min1, min2, max2 - min2, fun};

    const integrand_t cuhre_fun{[](const int * /* ndim */, const cubareal xx[],
                                   const int * /* ncomp */, cubareal ff[],
                                   void *userdata) -> int {
      auto i = static_cast<Integrand2d<F> *>(userdata);
      /* We have to transform the integrand to the unit cube.
       * This is what Cuba expects. */
      ff[0] = (i->f)(i->min1 + i->diff1 * xx[0], i->min2 + i->diff2 * xx[1]) *
              i->diff1 * i->diff2;
      return 0;
    }};

    const int ndim = 2;
    const int ncomp = 1;
    void *userdata = &f_with_limits;
    const int nvec = 1;
    const int flags = 0;  // Use the defaults.
    const int mineval = 0;
    const int maxeval = maxeval_;
    const int key = -1;  // Use the default.
    const char *statefile = nullptr;
    void *spin = nullptr;

    Cuhre(ndim, ncomp, cuhre_fun, userdata, nvec, epsrel_, epsabs_, flags,
          mineval, maxeval, key, statefile, spin, &nregions_, &neval_, &fail_,
          &result.first, &result.second, &prob_);

    if (fail_) {
      std::stringstream err;
      err << "After " << neval_ << " evaluations "
          << "Cuhre integration from Cuba reports error code " << fail_;
      throw std::runtime_error(err.str());
    }
    result.check_error("Cuba integration ", epsrel_, epsabs_);

    return result;
  }

 private:
  int maxeval_;
  double epsrel_;
  double epsabs_;
  int nregions_;
  int neval_;
  int fail_;
  double prob_;
};

}  // namespace smash

#endif  // SRC_INCLUDE_INTEGRATE_H_