MISERIntegrator.cpp

Go to the documentation of this file.

/*
 *  CepGen: a central exclusive processes event generator
 *  Copyright (C) 2013-2024  Laurent Forthomme
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#include <gsl/gsl_monte_miser.h>
 
#include "CepGen/Core/Exception.h"
#include "CepGen/Integration/GSLIntegrator.h"
#include "CepGen/Modules/IntegratorFactory.h"
 
namespace cepgen {
  class MISERIntegrator final : public GSLIntegrator {
  public:
    explicit MISERIntegrator(const ParametersList& params)
        : GSLIntegrator(params), ncvg_(steer<int>("numFunctionCalls")) {}
 
    static ParametersDescription description() {
      auto desc = GSLIntegrator::description();
      desc.setDescription("MISER adaptive importance sampling integrator");
      desc.add<int>("numFunctionCalls", 50'000)
          .setDescription("Number of function calls per phase space point evaluation");
      desc.add<double>("estimateFraction", 0.1);
      desc.add<int>("minCalls", 16 * 10);
      desc.add<int>("minCallsPerBisection", 32 * 16 * 10);
      desc.add<double>("alpha", 2.);
      desc.add<double>("dither", 0.1);
      return desc;
    }
 
    Value integrate(Integrand& integrand) override {
      setIntegrand(integrand);
      std::unique_ptr<gsl_monte_miser_state, decltype(&gsl_monte_miser_free)> miser_state(
          gsl_monte_miser_alloc(function_->dim), gsl_monte_miser_free);
      miser_state->verbose = verbosity_;
      gsl_monte_miser_params_get(miser_state.get(), &miser_params_);
      miser_params_.estimate_frac = steer<double>("estimateFraction");
      miser_params_.min_calls = steer<int>("minCalls");
      miser_params_.min_calls_per_bisection = steer<int>("minCallsPerBisection");
      miser_params_.alpha = steer<double>("alpha");
      miser_params_.dither = steer<double>("dither");
      gsl_monte_miser_params_set(miser_state.get(), &miser_params_);
 
      CG_DEBUG("Integrator:build") << "MISER parameters:\n\t"
                                   << "Number of calls: " << miser_params_.min_calls << ", "
                                   << "per bisection: " << miser_params_.min_calls_per_bisection << ",\n\t"
                                   << "Estimate fraction: " << miser_params_.estimate_frac << ",\n\t"
                                   << "α-value: " << miser_params_.alpha << ",\n\t"
                                   << "Dither: " << miser_params_.dither << ".";
 
      // launch the full integration
      double result, abserr;
      int res = gsl_monte_miser_integrate(function_.get(),
                                          &xlow_[0],
                                          &xhigh_[0],
                                          function_->dim,
                                          ncvg_,
                                          rnd_gen_->engine<gsl_rng>(),
                                          miser_state.get(),
                                          &result,
                                          &abserr);
 
      if (res != GSL_SUCCESS)
        throw CG_FATAL("Integrator:integrate") << "Error while performing the integration!\n\t"
                                               << "GSL error: " << gsl_strerror(res) << ".";
 
      return Value{result, abserr};
    }
 
  private:
    const int ncvg_;
    gsl_monte_miser_params miser_params_{};
  };
 
}  // namespace cepgen
REGISTER_INTEGRATOR("MISER", MISERIntegrator);