integrate_fp.h

Go to the documentation of this file.

#ifndef AMPGEN_INTEGRATE_FP
#define AMPGEN_INTEGRATE_FP 1
 
#include <tuple>
#include "AmpGen/simd/utils.h"
#include "AmpGen/MetaUtils.h"
 
namespace AmpGen {
  namespace integrate_fp_constants {
    static const double xl2 = 0.358568582800318073; // lambda_2
    static const double xl4 = 0.948683298050513796; // lambda_4
    static const double xl5 = 0.688247201611685289; // lambda_5
    static const double w2 = 980. / 6561;           // weights/2^n
    static const double w4 = 200. / 19683;
    static const double wp2 = 245. / 486; // error weights/2^n
    static const double wp4 = 25. / 729;
 
    static const double wn1[14] = {-0.193872885230909911, -0.555606360818980835, -0.876695625666819078, -1.15714067977442459, -1.39694152314179743,
                                   -1.59609815576893754,  -1.75461057765584494,  -1.87247878880251983,  -1.94970278920896201, -1.98628257887517146,
                                   -1.98221815780114818,  -1.93750952598689219,  -1.85215668343240347,  -1.72615963013768225};
 
    static const double wn3[14] = {0.0518213686937966768,  0.0314992633236803330,  0.0111771579535639891,  -0.00914494741655235473, -0.0294670527866686986,
                                   -0.0497891581567850424, -0.0701112635269013768, -0.0904333688970177241, -0.110755474267134071,   -0.131077579637250419,
                                   -0.151399685007366752,  -0.171721790377483099,  -0.192043895747599447,  -0.212366001117715794};
 
    static const double wn5[14]
      = {0.871183254585174982e-01, 0.435591627292587508e-01, 0.217795813646293754e-01, 0.108897906823146873e-01, 0.544489534115734364e-02,
         0.272244767057867193e-02, 0.136122383528933596e-02, 0.680611917644667955e-03, 0.340305958822333977e-03, 0.170152979411166995e-03,
         0.850764897055834977e-04, 0.425382448527917472e-04, 0.212691224263958736e-04, 0.106345612131979372e-04};
 
    static const double wpn1[14] = {-1.33196159122085045, -2.29218106995884763, -3.11522633744855959, -3.80109739368998611, -4.34979423868312742,
                                    -4.76131687242798352, -5.03566529492455417, -5.17283950617283939, -5.17283950617283939, -5.03566529492455417,
                                    -4.76131687242798352, -4.34979423868312742, -3.80109739368998611, -3.11522633744855959};
 
    static const double wpn3[14] = {0.0445816186556927292, -0.0240054869684499309, -0.0925925925925925875, -0.161179698216735251, -0.229766803840877915,
                                    -0.298353909465020564, -0.366941015089163228,  -0.435528120713305891,  -0.504115226337448555, -0.572702331961591218,
                                    -0.641289437585733882, -0.709876543209876532,  -0.778463648834019195,  -0.847050754458161859};
  }
#if INSTRUCTION_SET == INSTRUCTION_SET_AVX2d
  template <unsigned dim, unsigned width, typename fcn>
  std::tuple<double, double, unsigned> integrate_fp(const fcn &F, const std::array<double, dim> &ctr, const std::array<double, dim> &wth) {
    using namespace integrate_fp_constants;
    auto eval = [&F](auto z, const auto &p, const unsigned &ind) {
      z[ind] += p;
      return F(z);
    };
    auto eval_2 = [&F](auto z, const auto &p, const unsigned &ind, const auto &p2, const unsigned int &ind2) {
      z[ind] += p;
      z[ind2] += p2;
      return F(z);
    };
    unsigned idvaxn = 0;
    std::array<real_v, dim> z;
    double rgnvol = get_power<2, dim>::value;
    for(unsigned j = 0; j < dim; j++) {
      z[j] = ctr[j];
      rgnvol *= wth[j]; // region volume
    }
    double sum1{0}, sum2{0}, sum3{0}, sum4{0}, sum5{0};
    sum1 = F(z).at(0);
    double difmax = 0;
    const real_v xl(-xl4, -xl2, +xl2, +xl4);
    const real_v dxp(-xl4, +xl4, -xl4, +xl4);
    const real_v dxp2(-xl4, -xl4, +xl4, +xl4);
 
    // loop over coordinates
    for(unsigned j = 0; j != dim; j++) {
      auto fb = eval(z, xl * wth[j], j).to_array();
      auto f2 = fb[1] + fb[2];
      auto f3 = fb[0] + fb[3];
      sum2 += f2;
      sum3 += f3;
      double dif = std::abs(7 * f2 - f3 - 12 * sum1);
      idvaxn = dif >= difmax ? j : idvaxn;
      difmax = dif >= difmax ? dif : difmax;
      for(unsigned k = j + 1; k < dim; ++k) sum4 += utils::sum_elements(eval_2(z, dxp * wth[j], j, dxp2 * wth[k], k));
    }
    for(int j = 0; j < get_power<2, dim>::value; j += 4) {
      // auto loop_mask =
      for(int k = 0; k != dim; ++k)
        z[k] = real_v(ctr[k] + (0x1 & ((j + 0) >> k) ? +1 : -1) * xl5 * wth[k], ctr[k] + (0x1 & ((j + 1) >> k) ? +1 : -1) * xl5 * wth[k],
                      ctr[k] + (0x1 & ((j + 2) >> k) ? +1 : -1) * xl5 * wth[k], ctr[k] + (0x1 & ((j + 3) >> k) ? +1 : -1) * xl5 * wth[k]);
      sum5 += utils::sum_elements(F(z));
    }
 
    auto rgncmp = rgnvol * (wpn1[dim - 2] * sum1 + wp2 * sum2 + wpn3[dim - 2] * sum3 + wp4 * sum4);
    auto rgnval = rgnvol * (wn1[dim - 2] * sum1 + w2 * sum2 + wn3[dim - 2] * sum3 + w4 * sum4 + wn5[dim - 2] * sum5);
    auto rgnerr = std::abs(rgnval - rgncmp); // compares estim error with expected error
    return {rgnval, rgnerr, idvaxn};
  }
#endif
  template <unsigned dim, typename fcn>
  std::tuple<double, double, unsigned> integrate_fp_scalar(const fcn &F, const std::array<double, dim> &ctr, const std::array<double, dim> &wth) {
    using namespace integrate_fp_constants;
    auto eval = [&F](auto &z, const auto &p, const unsigned &ind) {
      auto zt = z[ind];
      z[ind] += p;
      auto v = F(z);
      z[ind] = zt;
      return v;
    };
    auto eval_2 = [&F](auto &z, const auto &p, const unsigned &ind, const auto &p2, const unsigned int &ind2) {
      auto zt = z[ind];
      auto zt2 = z[ind2];
      z[ind] += p;
      z[ind2] += p2;
      auto v = F(z);
      z[ind] = zt;
      z[ind2] = zt2;
      return v;
    };
 
    unsigned idvaxn = 0;
    std::array<double, dim> z = ctr;
    double rgnvol = get_power<2, dim>::value;
    for(unsigned j = 0; j < dim; j++) rgnvol *= wth[j]; // region volume
    double sum1 = F(z);
    double sum2 = 0;
    double sum3 = 0;
    double sum4 = 0;
    double sum5 = 0;
 
    double difmax = 0;
 
    // loop over coordinates
    for(unsigned j = 0; j != dim; j++) {
      auto f2 = eval(z, -xl2 * wth[j], j) + eval(z, +xl2 * wth[j], j);
      auto f3 = eval(z, -xl4 * wth[j], j) + eval(z, +xl4 * wth[j], j);
 
      sum2 += f2; // sum func eval with different weights separately
      sum3 += f3;
      double dif = std::abs(7 * f2 - f3 - 12 * sum1);
      idvaxn = dif >= difmax ? j : idvaxn;
      difmax = dif >= difmax ? dif : difmax;
      for(unsigned k = j + 1; k < dim; ++k) {
        sum4 += eval_2(z, xl4 * wth[j], j, xl4 * wth[k], k);
        sum4 += eval_2(z, xl4 * wth[j], j, -xl4 * wth[k], k);
        sum4 += eval_2(z, -xl4 * wth[j], j, xl4 * wth[k], k);
        sum4 += eval_2(z, -xl4 * wth[j], j, -xl4 * wth[k], k);
      }
    }
    for(int j = 0; j != get_power<2, dim>::value; ++j) {
      for(int k = 0; k != dim; ++k) z[k] = ctr[k] + (0x1 & (j >> k) ? +1 : -1) * xl5 * wth[k];
      sum5 += F(z);
    }
 
    auto rgncmp = rgnvol * (wpn1[dim - 2] * sum1 + wp2 * sum2 + wpn3[dim - 2] * sum3 + wp4 * sum4);
    auto rgnval = rgnvol * (wn1[dim - 2] * sum1 + w2 * sum2 + wn3[dim - 2] * sum3 + w4 * sum4 + wn5[dim - 2] * sum5);
    auto rgnerr = std::abs(rgnval - rgncmp); // compares estim error with expected error
    return {rgnval, rgnerr, idvaxn};
  }
 
}
 
#endif