Using ScalarIntegrationModule

In y3_cluster_cpp there is a cubacpp wrapper written by Marc that compiles into a cosmosis module to do high-dimensional numerical integration using the Cuhre algorithm.

It has been known throughout history as ScalarIntegrationModule

Some definitions

• Grid Point: Fixed numbers that you wish to evaluate your integral at. Some examples include radii at which to evaluate a cluster profile or cluster observables.
• Integration Volume: The total volume of integration defined by the lower and upper limits of each integration variable.

Configuration Parameters

In your cosmosis ini file, you will need to include configuration several parameters for your integration module.

The following configuration parameters are required every time

• file: (str) The path to the compiled .so for your module.
• eps_rel: (double) The relative difference between successive integration steps to signify convergence.
• eps_abs: (double) The absolute difference between successive integration steps to signify convergence.
• max_eval: (int) The maximum number of times to evaluate your integrand before giving up on convergence.

• use_cartesian_product: (bool)

  • If True, your integral is evaluated at every combination of grid points and integration volumes for each MCMC step.
  • If False, your integral is evaluated for each element in the zipped sequence of grid points and integration volumes. You need to define an equal number of grid points and integration volumes.

In addition to the parameters above, you will also need to define grid points and integration volumes. Integration volumes are defined by appending _low and _high to the string that defines the variable of integration. For example, given two integration variables x and y, you define the integration volumes by writing

x_low = 0.2 0.4
x_high = 0.4 0.6
y_low = 1. 4.
y_high = 4. 8.

in your cosmosis ini file. To define grid points, you simply use the string that defines the variable. If you want to evaluate your integral at three different values of a grid points radii, you would write

radii = 0.5 0.75

in your cosmosis ini file.

There are two provided functions for reading integration volumes and grid points from the ini files, wall_of_numbers and cartesian_product.

If you read the volumes above using the wall_of_numbers scheme, you will evaluate the integral over two volumes at each MCMC step:

• x\in[0.2, 0.4], y\in[1, 4]
• x\in[0.4, 0.6], y\in[4, 8]

If you read the volumes above using the cartesian_product scheme, you will evaluate the integral over four volumes at each MCMC step:

• x\in[0.2, 0.4], y\in[1, 4]
• x\in[0.4, 0.6], y\in[1, 4]
• x\in[0.2, 0.4], y\in[4, 8]
• x\in[0.4, 0.6], y\in[4, 8]

If you have multiple grid points, the same logic applies when you use these two different functions.

Finally, if you set use_cartesian_product = False and read in both your volumes and grid points with the wall_of_numbers scheme, you will evaluate two integrations in this example per MCMC step:

• x\in[0.2, 0.4], y\in[1, 4] \mathrm{with\ radii} = 0.5
• x\in[0.4, 0.6], y\in[4, 8] \mathrm{with\ radii} = 0.75

And if you set use_cartesian_product = True and read in both your volumes and grid points with the wall_of_numbers scheme, you will evaluate four integrations in this example per MCMC step:

• x\in[0.2, 0.4], y\in[1, 4] \mathrm{with\ radii} = 0.5
• x\in[0.4, 0.6], y\in[4, 8] \mathrm{with\ radii} = 0.5
• x\in[0.2, 0.4], y\in[1, 4] \mathrm{with\ radii} = 0.75
• x\in[0.4, 0.6], y\in[4, 8] \mathrm{with\ radii} = 0.75

Integration Using ScalarIntegrationModule

This wrapper assumes that you have written a class that contains several member variables and member functions with predefined names and call signatures. There is an example class that we can just copy and fill in details where needed to make the process as simple as possible.

1. To begin, copy the contents of src/modules/MattsExampleScalarIntegrand.cc into a new file. This file will contain integration module so chose a descriptive name!

#!c++
class ExampleScalarIntegrand {

2. This is simply the name of your class, chose a descriptive class name! Everywhere you see ExampleScalarIntegrand below, replace with your class name.

#!c++
public:
  // Define the data-type describing a grid point; this should be an
  // instance of std::array<double, N> with N set to the number
  // of different parameters being varied in the grid.
  // The alias we define must be grid_point_t.
  using grid_point_t = std::array<double, 1>; // we only vary radius.

3. This is where you specify the number of grid dimensions that you have. All you need to change here is the number in std::array<double, N> from 1 to the number of grid dimensions that you have.

#!c++
private:
  // We define the type alias volume_t to be the right dimensionality
  // of integration volume for our integrand. If we were to change the
  // number of arguments required by the function call operator (below),
  // we would need to also modify this type alias to keep consistent.
  using volume_t = cubacpp::IntegrationVolume<2>;

4. This is where you specify the number of dimensions you wish to integrate over. Change cubacpp::IntegrationVolume<N> from 2 to the number of dimensions you wish to integrate over.

#!c++
  // State obtained from configuration. These things should be set in the
  // constructor.
  // <none in this example>

5. This is where you declare anything that you need to compute the integrand that doesn't depend on any model parameters (or anything else that varies between MCMC steps). These may be read in from configuration files or could be types that do not depend on model parameters.

#!c++
  // State obtained from each sample.
  // If there were a type X that did not have a default constructor,
  // we would use std::optional<X> as our data member.
  double sigma_8_;
  std::optional<HMF_t> hmf_;

6. In this section, we declare anything that depends on model parameters varied between MCMC steps. Some things do not make sense to be declared without a cosmology (or other model parameters) and thus will not have a default constructor. These should be declared using std::optional<X>.

#!c++
  // State set for current 'bin' to be integrated.
  double radius_;

7. Next up are the grid points. In this example we have only 1 grid point, radius. If you have several grid points, they should all be declared here.

#!c++
  public:
  // Initialize my integrand object from the parameters read
  // from the relevant block in the CosmoSIS ini file.
  explicit ExampleScalarIntegrand(cosmosis::DataBlock& config);

#!c++
ExampleScalarIntegrand::ExampleScalarIntegrand(DataBlock&)
  : sigma_8_(), hmf_(), radius_()
{}

8. This is the class constructor and is run only once during the setup stage of a cosmosis run. This should read in any configuration information and initialize all of the member variables.

#!c++
  // Set any data members from values read from the current sample.
  // Do not attempt to copy the sample!.
  void set_sample(cosmosis::DataBlock& sample);

#!c++
void
ExampleScalarIntegrand::set_sample(DataBlock& sample)
{
  sigma_8_ = sample.view<double>("cosmological_parameters", "sigma_8");
  // If we had a data member of type std::optional<X>, we would set the
  // value using std::optional::emplace(...) here. emplace takes a set
  // of arguments that it passes to the constructor of X.
  hmf_.emplace(sample);
}

9. This function updates the members that depend on the model parameters with the current MCMC step. This function should update each member declared in step 6. Note the difference between members declared using their type and using std::optional.

#!c++
  // Set the data for the current bin.
  void set_grid_point(grid_point_t const& pt);

#!c++
void
ExampleScalarIntegrand::set_grid_point(grid_point_t const& grid_point)
{
  radius_ = grid_point[0];
}

10. This function updates each member that corresponds to the grid points that were declared in step 7. Note that the vector grid_point has a specific ordering set below in step 14 and you should ensure that these orderings are consistent.

#!c++
  // The function to be integrated. All arguments to this function must be of
  // type double, and there must be at least two of them (because our
  // integration routine does not work for functions of one variable). The
  // function is const because calling it does not change the state of the
  // object.
  double operator()(double x, double y) const;

#!c++
// This math is totally non-physical and stupid, but it uses all the values
// provided.
double
ExampleScalarIntegrand::operator()(double x, double y) const
{
  // For any data members of type std::optional<X>, we have to use operator*
  // to access the X object (as if we were dereferencing a pointer).
  auto const delta = y - sigma_8_;
  return (x / radius_) + (delta * delta);
}

11. This function overloads the () operator and should contain your integrand. This function is passed each of your integration variables in the order specified in step 13 and you should ensure that the ordering between these two steps is consistent. Should return the value of your integrand.

#!c++
  // module_label() is a non-member (static) function that returns the label for
  // this module. The name this returns
  // is the name that must be used in the 'ini file' for configuring the module
  // made with this class.
  // We return char const* rather than std::string to avoid some needless memory
  // allocations.
  static char const* module_label();

#!c++
char const*
ExampleScalarIntegrand::module_label()
{
  return "example_scalar_integrand";
}

12. This is the name of the module. It is used to save data to the datablock and to obtain configuration information from the ini file. You should chose a descriptive name that will not collide with other modules in your cosmosis pipeline.

#!c++
  // The following non-member (static) function creates a vector of integration
  // volumes (the type alias defined above) based on the parameters read from
  // the configuration block for the module.
  static std::vector<volume_t> make_integration_volumes(
    cosmosis::DataBlock& cfg);

#!c++
std::vector<ExampleScalarIntegrand::volume_t>
ExampleScalarIntegrand::make_integration_volumes(cosmosis::DataBlock& cfg)
{
  return y3_cluster::make_integration_volumes_wall_of_numbers(
    cfg, ExampleScalarIntegrand::module_label(), "x", "y");
}

13. This function returns the integration volumes to use. There are two functions defined in src/utils/make_integration_volumes.hh that should cover the vast majority of use cases. You can write your own function if you need to by simply copying what is done in these functions. The two functions, make_integration_volumes_cartesian_product and make_integration_volumes_wall_of_numbers, are detailed in the previous sections above. Both functions take at least 2 parameters: cfg (the datablock configuration information), the module label, and then a string for each integration dimension with the same name as the relevant configuration parameters (see above). The order that these dimensions are specified are the order that they are passed into your integrand.

#!c++
  // The following non-member (static) function creates a vector of grid points
  // on which the integration results are to be evaluated, based on parameters
  // read from the configuration block for the module.
  static std::vector<grid_point_t> make_grid_points(cosmosis::DataBlock& cfg);
};

#!c++
std::vector<ExampleScalarIntegrand::grid_point_t>
ExampleScalarIntegrand::make_grid_points(cosmosis::DataBlock& cfg)
{
  return y3_cluster::make_grid_points_cartesian_product(
    cfg, ExampleScalarIntegrand::module_label(), "radii");
}

14. This function reads the grid points from the configuration file. There are two functions defined in src/utils/make_grid_points.hh that should cover the vast majority of use cases. You can write your own function if you need to by simply copying what is done in these functions. The two functions, make_grid_points_cartesian_product and make_grid_points_wall_of_numbers, are detailed in the previous sections above. Both functions take at least 2 parameters: cfg (the datablock configuration information, the module label, and then a string for each grid point that can be found in the configuration file. The order that these grid points are specified are the order that they appear in the vector `grid_points` in step 10 above.

#!c++
DEFINE_COSMOSIS_SCALAR_INTEGRATION_MODULE(ExampleScalarIntegrand)

15. Finally, replace the name of your class in this line at the end of the file. This line of code is what actually converts your class into the integration module.

16. The last step is to add your module to src/modules/CMakeLists.txt

#!bash
add_library(ExampleScalarIntegrand MODULE MattsExampleScalarIntegrand.cc)
target_link_libraries(ExampleScalarIntegrand PRIVATE cosmosis utils models
                      ${GSL_LIBRARIES} ${CUBA_LIBRARIES})
target_include_directories(ExampleScalarIntegrand PRIVATE ${CMAKE_SOURCE_DIR}/src
                           ${CUBACPP_DIR} ${GSL_INCLUDE_DIRS} ${EIGEN3_INCLUDE_DIR})