From 6a3e8ea978ad98ec22418abb0611c8738f333efc Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 14 Apr 2023 18:17:28 +0200 Subject: [PATCH 01/12] Port polygamma to CFFI --- cern_polygamma.c | 117 +++++++++++++++++++++++++++++++ src/ekore/cffilib.py | 22 ++++++ src/ekore/harmonics/polygamma.py | 14 +++- 3 files changed, 152 insertions(+), 1 deletion(-) create mode 100644 cern_polygamma.c create mode 100644 src/ekore/cffilib.py diff --git a/cern_polygamma.c b/cern_polygamma.c new file mode 100644 index 000000000..156601630 --- /dev/null +++ b/cern_polygamma.c @@ -0,0 +1,117 @@ +// gcc -shared -o cernlibc.so -fPIC cern_polygamma.c + +#include +#include + +double getdouble(double* ptr) { + return *ptr; +} + +int cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { + const double _Complex Z = re_in + im_in * _Complex_I; + const double DELTA = 5e-13; + const double R1 = 1; + const double HF = R1/2; + const double C1 = pow(M_PI,2); + const double C2 = 2*pow(M_PI,3); + const double C3 = 2*pow(M_PI,4); + const double C4 = 8*pow(M_PI,5); + + // SGN is originally indexed 0:4 -> no shift + const double SGN[] = {-1,+1,-1,+1,-1}; + // FCT is originally indexed -1:4 -> shift +1 + const double FCT[] = {0,1,1,2,6,24}; + + // C is originally indexed 1:6 x 0:4 -> swap indices and shift new last -1 + const double C[5][6] = { + { + 8.33333333333333333e-2, + -8.33333333333333333e-3, + 3.96825396825396825e-3, + -4.16666666666666667e-3, + 7.57575757575757576e-3, + -2.10927960927960928e-2 + }, + { + 1.66666666666666667e-1, + -3.33333333333333333e-2, + 2.38095238095238095e-2, + -3.33333333333333333e-2, + 7.57575757575757576e-2, + -2.53113553113553114e-1 + }, + { + 5.00000000000000000e-1, + -1.66666666666666667e-1, + 1.66666666666666667e-1, + -3.00000000000000000e-1, + 8.33333333333333333e-1, + -3.29047619047619048e+0 + }, + { + 2.00000000000000000e+0, + -1.00000000000000000e+0, + 1.33333333333333333e+0, + -3.00000000000000000e+0, + 1.00000000000000000e+1, + -4.60666666666666667e+1 + }, + {10., -7., 12., -33., 130., -691.} + }; + double _Complex U=Z; + double X=creal(U); + double A=cabs(X); + if (K < 0 || K > 4) + return -1; + // raise NotImplementedError("Order K has to be in [0:4]") + const int A_as_int = A; + if (cabs(cimag(U)) < DELTA && cabs(X+A_as_int) < DELTA) + return -2; + // raise ValueError("Argument Z equals non-positive integer") + const unsigned int K1=K+1; + if (X < 0) + U=-U; + double _Complex V=U; + double _Complex H=0; + if (A < 15) { + H=1/cpow(V,K1); + for (int i = 1; i < 14 - A_as_int + 1 ; ++i) { + V=V+1; + H=H+1/cpow(V,K1); + } + V=V+1; + } + double _Complex R=1/cpow(V,2); + double _Complex P=R*C[K][6-1]; + for (int i = 5; i>1-1; i--) + P=R*(C[K][i-1]+P); + H=SGN[K]*(FCT[K+1]*H+(V*(FCT[K-1+1]+P)+HF*FCT[K+1])/cpow(V,K1)); + if (0 == K) + H=H+clog(V); + if (X < 0){ + V=M_PI*U; + X=creal(V); + const double Y=cimag(V); + const double A=sin(X); + const double B=cos(X); + const double T=tanh(Y); + P=(B-A*T*_Complex_I)/(A+B*T*_Complex_I); + if (0 == K) + H=H+1/U+M_PI*P; + else if (1 == K) + H=-H+1/cpow(U,2)+C1*(cpow(P,2)+1); + else if (2 == K) + H=H+2/cpow(U,3)+C2*P*(cpow(P,2)+1); + else if (3 == K) { + R=cpow(P,2); + H=-H+6/cpow(U,4)+C3*((3*R+4)*R+1); + } else if (4 == K) { + R=cpow(P,2); + H=H+24/cpow(U,5)+C4*P*((3*R+5)*R+2); + } + } + //return creal(H); + *re_out = creal(H); + *im_out = cimag(H); + return 0; +} diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py new file mode 100644 index 000000000..ef98d65a6 --- /dev/null +++ b/src/ekore/cffilib.py @@ -0,0 +1,22 @@ +"""CFFI interface.""" +from cffi import FFI + +ffi = FFI() +cernlibc = ffi.dlopen("./cernlibc.so") + +ffi.cdef( + """ +int cern_polygamma(const double re_in, const double im_in, const unsigned int K, void* re_out, void* im_out); +double getdouble(void* ptr); +""" +) + +# we need to "activate" the actual function first +c_cern_polygamma = cernlibc.cern_polygamma +c_getdouble = cernlibc.getdouble + +# allocate the pointers and get their addresses +re_y = ffi.new("double*") +im_y = ffi.new("double*") +re_y_address = int(ffi.cast("uintptr_t", re_y)) +im_y_address = int(ffi.cast("uintptr_t", im_y)) diff --git a/src/ekore/harmonics/polygamma.py b/src/ekore/harmonics/polygamma.py index a2c7e442f..a03620da2 100644 --- a/src/ekore/harmonics/polygamma.py +++ b/src/ekore/harmonics/polygamma.py @@ -6,9 +6,21 @@ import numba as nb import numpy as np +from ..cffilib import c_cern_polygamma, c_getdouble, im_y_address, re_y_address + + +@nb.njit() +def cern_polygamma(z, k): + """Compute polygamma from C.""" + # bypass complex and pointers + res = c_cern_polygamma(np.real(z), np.imag(z), k, re_y_address, im_y_address) + re_y = c_getdouble(re_y_address) + im_y = c_getdouble(im_y_address) + return re_y + im_y * 1j + @nb.njit(cache=True) -def cern_polygamma(Z, K): # pylint: disable=all +def _cern_polygamma(Z, K): # pylint: disable=all r"""Compute the polygamma functions :math:`\psi_k(z)`. Reimplementation of ``WPSIPG`` (C317) in `CERNlib `_ From 80483348d7a3d8be6aadec80373a3151358fdaf9 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 10 Jul 2023 17:06:36 +0200 Subject: [PATCH 02/12] Move to c++ --- cern_polygamma.c => cern_polygamma.hpp | 72 +++++++++++--------------- cinterface.cpp | 22 ++++++++ src/ekore/cffilib.py | 10 ++-- 3 files changed, 57 insertions(+), 47 deletions(-) rename cern_polygamma.c => cern_polygamma.hpp (60%) create mode 100644 cinterface.cpp diff --git a/cern_polygamma.c b/cern_polygamma.hpp similarity index 60% rename from cern_polygamma.c rename to cern_polygamma.hpp index 156601630..d47083113 100644 --- a/cern_polygamma.c +++ b/cern_polygamma.hpp @@ -1,14 +1,7 @@ -// gcc -shared -o cernlibc.so -fPIC cern_polygamma.c +#include +#include -#include -#include - -double getdouble(double* ptr) { - return *ptr; -} - -int cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { - const double _Complex Z = re_in + im_in * _Complex_I; +std::complex cern_polygamma(const std::complex Z, const unsigned int K) { const double DELTA = 5e-13; const double R1 = 1; const double HF = R1/2; @@ -58,60 +51,55 @@ int cern_polygamma(const double re_in, const double im_in, const unsigned int K, }, {10., -7., 12., -33., 130., -691.} }; - double _Complex U=Z; - double X=creal(U); - double A=cabs(X); + std::complex U=Z; + double X=U.real(); + double A=abs(X); if (K < 0 || K > 4) - return -1; - // raise NotImplementedError("Order K has to be in [0:4]") + throw std::invalid_argument("Order K has to be in [0:4]"); const int A_as_int = A; - if (cabs(cimag(U)) < DELTA && cabs(X+A_as_int) < DELTA) - return -2; - // raise ValueError("Argument Z equals non-positive integer") + if (abs(U.imag()) < DELTA && abs(X+A_as_int) < DELTA) + throw std::invalid_argument("Argument Z equals non-positive integer"); const unsigned int K1=K+1; if (X < 0) U=-U; - double _Complex V=U; - double _Complex H=0; + std::complex V=U; + std::complex H=0; if (A < 15) { - H=1/cpow(V,K1); + H=1./pow(V,K1); for (int i = 1; i < 14 - A_as_int + 1 ; ++i) { - V=V+1; - H=H+1/cpow(V,K1); + V=V+1.; + H=H+1./pow(V,K1); } - V=V+1; + V=V+1.; } - double _Complex R=1/cpow(V,2); - double _Complex P=R*C[K][6-1]; + std::complex R=1./pow(V,2); + std::complex P=R*C[K][6-1]; for (int i = 5; i>1-1; i--) P=R*(C[K][i-1]+P); - H=SGN[K]*(FCT[K+1]*H+(V*(FCT[K-1+1]+P)+HF*FCT[K+1])/cpow(V,K1)); + H=SGN[K]*(FCT[K+1]*H+(V*(FCT[K-1+1]+P)+HF*FCT[K+1])/pow(V,K1)); if (0 == K) - H=H+clog(V); + H=H+log(V); if (X < 0){ V=M_PI*U; - X=creal(V); - const double Y=cimag(V); + X=V.real(); + const double Y=V.imag(); const double A=sin(X); const double B=cos(X); const double T=tanh(Y); - P=(B-A*T*_Complex_I)/(A+B*T*_Complex_I); + P=std::complex(B,-A*T)/std::complex(A,+B*T); if (0 == K) - H=H+1/U+M_PI*P; + H=H+1./U+M_PI*P; else if (1 == K) - H=-H+1/cpow(U,2)+C1*(cpow(P,2)+1); + H=-H+1./pow(U,2)+C1*(pow(P,2)+1.); else if (2 == K) - H=H+2/cpow(U,3)+C2*P*(cpow(P,2)+1); + H=H+2./pow(U,3)+C2*P*(pow(P,2)+1.); else if (3 == K) { - R=cpow(P,2); - H=-H+6/cpow(U,4)+C3*((3*R+4)*R+1); + R=pow(P,2); + H=-H+6./pow(U,4)+C3*((3.*R+4.)*R+1.); } else if (4 == K) { - R=cpow(P,2); - H=H+24/cpow(U,5)+C4*P*((3*R+5)*R+2); + R=pow(P,2); + H=H+24./pow(U,5)+C4*P*((3.*R+5.)*R+2.); } } - //return creal(H); - *re_out = creal(H); - *im_out = cimag(H); - return 0; + return H; } diff --git a/cinterface.cpp b/cinterface.cpp new file mode 100644 index 000000000..d30a2c5ce --- /dev/null +++ b/cinterface.cpp @@ -0,0 +1,22 @@ +// g++ -shared -o cernlibc2.so -fPIC cinterface.cpp + +#include +#include +#include +#include "cern_polygamma.hpp" + +extern "C" double c_getdouble(double* ptr) { + return *ptr; +} + +extern "C" int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { + const std::complex Z (re_in, im_in); + try { + const std::complex H = cern_polygamma(Z, K); + *re_out = H.real(); + *im_out = H.imag(); + } catch (std::invalid_argument &e) { + return -1; + } + return 0; +} diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index ef98d65a6..451b1362c 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -2,18 +2,18 @@ from cffi import FFI ffi = FFI() -cernlibc = ffi.dlopen("./cernlibc.so") +cernlibc = ffi.dlopen("./cernlibc2.so") ffi.cdef( """ -int cern_polygamma(const double re_in, const double im_in, const unsigned int K, void* re_out, void* im_out); -double getdouble(void* ptr); +int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, void* re_out, void* im_out); +double c_getdouble(void* ptr); """ ) # we need to "activate" the actual function first -c_cern_polygamma = cernlibc.cern_polygamma -c_getdouble = cernlibc.getdouble +c_cern_polygamma = cernlibc.c_cern_polygamma +c_getdouble = cernlibc.c_getdouble # allocate the pointers and get their addresses re_y = ffi.new("double*") From 6146a7c129f725ba38befda6d04154ed88983113 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 10 Jul 2023 18:52:31 +0200 Subject: [PATCH 03/12] Implement ad/u/s/as1.gamma_ns --- cinterface.cpp | 22 ---------- .../unpolarized/space_like/__init__.py | 27 +++++++++++++ src/ekore/cffilib.py | 9 ++++- .../unpolarized/space_like/as1.hpp | 21 ++++++++++ .../unpolarized/space_like/init.hpp | 28 +++++++++++++ src/ekorepp/cinterface.cpp | 40 +++++++++++++++++++ src/ekorepp/constants.hpp | 27 +++++++++++++ .../ekorepp/harmonics/polygamma.hpp | 25 ++++++++++++ src/ekorepp/harmonics/w1.hpp | 22 ++++++++++ 9 files changed, 198 insertions(+), 23 deletions(-) delete mode 100644 cinterface.cpp create mode 100644 src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp create mode 100644 src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp create mode 100644 src/ekorepp/cinterface.cpp create mode 100644 src/ekorepp/constants.hpp rename cern_polygamma.hpp => src/ekorepp/harmonics/polygamma.hpp (78%) create mode 100644 src/ekorepp/harmonics/w1.hpp diff --git a/cinterface.cpp b/cinterface.cpp deleted file mode 100644 index d30a2c5ce..000000000 --- a/cinterface.cpp +++ /dev/null @@ -1,22 +0,0 @@ -// g++ -shared -o cernlibc2.so -fPIC cinterface.cpp - -#include -#include -#include -#include "cern_polygamma.hpp" - -extern "C" double c_getdouble(double* ptr) { - return *ptr; -} - -extern "C" int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { - const std::complex Z (re_in, im_in); - try { - const std::complex H = cern_polygamma(Z, K); - *re_out = H.real(); - *im_out = H.imag(); - } catch (std::invalid_argument &e) { - return -1; - } - return 0; -} diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 6a3c6cf06..72335de0d 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -20,10 +20,37 @@ from eko import constants +from ....cffilib import ( + c_ad_us_gamma_ns, + c_getdouble, + im_double_3_address, + re_double_3_address, +) from ....harmonics import cache as c from . import aem1, aem2, as1, as1aem1, as2, as3, as4 +@nb.njit() +def c_gamma_ns(order, mode, n, nf): + """Compute gamma_ns from C.""" + # bypass complex and pointers + res = c_ad_us_gamma_ns( + order[0], + mode, + np.real(n), + np.imag(n), + nf, + re_double_3_address, + im_double_3_address, + ) + gamma_ns = np.zeros(order[0], np.complex_) + for j in range(order[0]): + gamma_ns[j] = ( + c_getdouble(re_double_3_address) + c_getdouble(im_double_3_address) * 1j + ) + return gamma_ns + + @nb.njit(cache=True) def gamma_ns(order, mode, n, nf): r"""Compute the tower of the non-singlet anomalous dimensions. diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index 451b1362c..2d0626221 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -2,21 +2,28 @@ from cffi import FFI ffi = FFI() -cernlibc = ffi.dlopen("./cernlibc2.so") +cernlibc = ffi.dlopen("./ekorepp.so") ffi.cdef( """ int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, void* re_out, void* im_out); double c_getdouble(void* ptr); +int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); """ ) # we need to "activate" the actual function first c_cern_polygamma = cernlibc.c_cern_polygamma c_getdouble = cernlibc.c_getdouble +c_ad_us_gamma_ns = cernlibc.c_ad_us_gamma_ns # allocate the pointers and get their addresses re_y = ffi.new("double*") im_y = ffi.new("double*") re_y_address = int(ffi.cast("uintptr_t", re_y)) im_y_address = int(ffi.cast("uintptr_t", im_y)) + +re_double_3 = ffi.new("double[3]") +im_double_3 = ffi.new("double[3]") +re_double_3_address = int(ffi.cast("uintptr_t", re_double_3)) +im_double_3_address = int(ffi.cast("uintptr_t", im_double_3)) diff --git a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp new file mode 100644 index 000000000..47f5b33ab --- /dev/null +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp @@ -0,0 +1,21 @@ +#include "../../../constants.hpp" +#include "../../../harmonics/w1.hpp" + + +namespace ekorepp { +namespace anomalous_dimensions { +namespace unpolarized { +namespace spacelike { +namespace as1 { + +std::complex gamma_ns(std::complex N, const unsigned int _nf) { + const std::complex S1 = harmonics::S1(N); + const std::complex gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0)); + return CF * gamma; +} + +} // namepsace as1 +} // namespace spacelike +} // namespace unpolarized +} // namespace anomalous_dimensions +} // namespace ekorepp diff --git a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp new file mode 100644 index 000000000..6beeae806 --- /dev/null +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp @@ -0,0 +1,28 @@ +#include +#include "./as1.hpp" +#include "../../../harmonics/polygamma.hpp" + + +namespace ekorepp { +namespace anomalous_dimensions { +namespace unpolarized { +namespace spacelike { + +/** + * @brief Compute the tower of the non-singlet anomalous dimensions. + * @param order_qcd QCD perturbative orders + * @param mode sector identifier + * @param n Mellin variable + * @param nf number of active flavors + * @return non-singlet anomalous dimensions + */ +std::vector> gamma_ns(const unsigned int order_qcd, const unsigned int mode, std::complex n, const unsigned int nf) { + std::vector> gamma_ns(order_qcd); + gamma_ns[0] = as1::gamma_ns(n, nf); + return gamma_ns; +} + +} // namespace spacelike +} // namespace unpolarized +} // namespace anomalous_dimensions +} // namespace ekorepp diff --git a/src/ekorepp/cinterface.cpp b/src/ekorepp/cinterface.cpp new file mode 100644 index 000000000..9a243ce67 --- /dev/null +++ b/src/ekorepp/cinterface.cpp @@ -0,0 +1,40 @@ +// g++ -shared -o ekorepp.so -fPIC src/ekorepp/cinterface.cpp + +#include +#include +#include +#include +#include "./harmonics/polygamma.hpp" +#include "./anomalous_dimensions/unpolarized/space_like/init.hpp" + +extern "C" double c_getdouble(double* ptr) { + return *ptr; +} + +extern "C" int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { + const std::complex Z (re_in, im_in); + try { + const std::complex H = ekorepp::harmonics::cern_polygamma(Z, K); + *re_out = H.real(); + *im_out = H.imag(); + } catch (std::invalid_argument &e) { + return -1; + } + return 0; +} + + +extern "C" int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, + const double re_in, const double im_in, const unsigned int nf, double* re_out, double* im_out) { + const std::complex n (re_in, im_in); + try { + const std::vector> res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, n, nf); + for (unsigned int j = 0; j < res.size(); ++j) { + re_out[j] = res[j].real(); + im_out[j] = res[j].imag(); + } + } catch (std::invalid_argument &e) { + return -1; + } + return 0; +} diff --git a/src/ekorepp/constants.hpp b/src/ekorepp/constants.hpp new file mode 100644 index 000000000..f1d2b6f00 --- /dev/null +++ b/src/ekorepp/constants.hpp @@ -0,0 +1,27 @@ +namespace ekorepp { + +/** + * @brief The number of colors. + * Defaults to :math:`N_C = 3`. + */ +const unsigned int NC = 3; + +/** + * @brief The normalization of fundamental generators. + * Defaults to :math:`T_R = 1/2`. + */ +const double TR = 1.0 / 2.0; + +/** + * @brief Second Casimir constant in the adjoint representation. + * Defaults to :math:`N_C = 3`. + */ +const double CA = double(NC); + +/** + * @brief Second Casimir constant in the fundamental representation. + * Defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`. + */ +const double CF = double(double(NC * NC - 1) / (2.0 * NC)); + +} diff --git a/cern_polygamma.hpp b/src/ekorepp/harmonics/polygamma.hpp similarity index 78% rename from cern_polygamma.hpp rename to src/ekorepp/harmonics/polygamma.hpp index d47083113..7deca2af7 100644 --- a/cern_polygamma.hpp +++ b/src/ekorepp/harmonics/polygamma.hpp @@ -1,6 +1,19 @@ +#ifndef EKOREPP_HARMONICS_POLYGAMMA_HPP_ +#define EKOREPP_HARMONICS_POLYGAMMA_HPP_ + #include #include +namespace ekorepp { +namespace harmonics { + +/** + * @brief Compute the polygamma functions :math:`\psi_k(z)`. + * Reimplementation of ``WPSIPG`` (C317) in `CERNlib `_ :cite:`KOLBIG1972221`. + * @param Z argument of polygamma function + * @param K order of polygamma function + * @return k-th polygamma function :math:`\psi_k(z)` + */ std::complex cern_polygamma(const std::complex Z, const unsigned int K) { const double DELTA = 5e-13; const double R1 = 1; @@ -103,3 +116,15 @@ std::complex cern_polygamma(const std::complex Z, const unsigned } return H; } + +std::complex recursive_harmonic_sum(const std::complex base_value, const std::complex n, const unsigned int iterations, const unsigned int weight) { + std::complex fact = 0.0; + for (unsigned int i = 1; i <= iterations + 1; ++i) + fact += pow(1.0 / (n + double(i)), weight); + return base_value + fact; +} + +} // namespace harmonics +} // namespace ekorepp + +#endif // EKOREPP_HARMONICS_POLYGAMMA_HPP_ diff --git a/src/ekorepp/harmonics/w1.hpp b/src/ekorepp/harmonics/w1.hpp new file mode 100644 index 000000000..59c6f00f7 --- /dev/null +++ b/src/ekorepp/harmonics/w1.hpp @@ -0,0 +1,22 @@ +#include +#include "./polygamma.hpp" + +namespace ekorepp { +namespace harmonics { + +/** + * @brief Compute the harmonic sum :math:`S_1(N)`. + * .. math:: + * S_1(N) = \sum\limits_{j=1}^N \frac 1 j = \psi_0(N+1)+\gamma_E + * + * with :math:`\psi_0(N)` the digamma function and :math:`\gamma_E` the + * Euler-Mascheroni constant. + * @param N Mellin moment + * @return (simple) Harmonic sum :math:`S_1(N)` + */ +std::complex S1(const std::complex N) { + return cern_polygamma(N + 1.0, 0) + 0.5772156649015328606065120; +} + +} // namespace harmonics +} // namespace ekorepp From 0a435431f7caf7ad4003e637cc3be6e1f9fbecbb Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 11 Jul 2023 13:03:47 +0200 Subject: [PATCH 04/12] Add ekorepp/ad/u/s/gamma_singlet --- .../unpolarized/space_like/__init__.py | 37 ++++++++-- src/ekore/cffilib.py | 20 +++--- src/ekore/harmonics/polygamma.py | 14 +--- .../unpolarized/space_like/as1.hpp | 68 +++++++++++++++++-- .../unpolarized/space_like/init.hpp | 29 +++++--- src/ekorepp/cinterface.cpp | 42 +++++++----- src/ekorepp/harmonics/cache.hpp | 63 +++++++++++++++++ src/ekorepp/harmonics/polygamma.hpp | 19 +++--- src/ekorepp/harmonics/w1.hpp | 11 ++- src/ekorepp/types.hpp | 16 +++++ 10 files changed, 251 insertions(+), 68 deletions(-) create mode 100644 src/ekorepp/harmonics/cache.hpp create mode 100644 src/ekorepp/types.hpp diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 72335de0d..57545495c 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -22,9 +22,12 @@ from ....cffilib import ( c_ad_us_gamma_ns, + c_ad_us_gamma_singlet, c_getdouble, - im_double_3_address, - re_double_3_address, + im_double_4_address, + im_double_16_address, + re_double_4_address, + re_double_16_address, ) from ....harmonics import cache as c from . import aem1, aem2, as1, as1aem1, as2, as3, as4 @@ -40,17 +43,41 @@ def c_gamma_ns(order, mode, n, nf): np.real(n), np.imag(n), nf, - re_double_3_address, - im_double_3_address, + re_double_4_address, + im_double_4_address, ) gamma_ns = np.zeros(order[0], np.complex_) for j in range(order[0]): gamma_ns[j] = ( - c_getdouble(re_double_3_address) + c_getdouble(im_double_3_address) * 1j + c_getdouble(re_double_4_address + j * 8) + + c_getdouble(im_double_4_address + j * 8) * 1j ) return gamma_ns +@nb.njit() +def c_gamma_singlet(order, n, nf): + """Compute gamma_singlet from C.""" + # bypass complex and pointers + res = c_ad_us_gamma_singlet( + order[0], + np.real(n), + np.imag(n), + nf, + re_double_16_address, + im_double_16_address, + ) + size = order[0] * 4 + gamma_s = np.zeros(size, np.complex_) + for j in range(size): + gamma_s[j] = ( + c_getdouble(re_double_16_address + j * 8) + + c_getdouble(im_double_16_address + j * 8) * 1j + ) + + return gamma_s.reshape(order[0], 2, 2) + + @nb.njit(cache=True) def gamma_ns(order, mode, n, nf): r"""Compute the tower of the non-singlet anomalous dimensions. diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index 2d0626221..c6a031e35 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -6,24 +6,24 @@ ffi.cdef( """ -int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, void* re_out, void* im_out); double c_getdouble(void* ptr); int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); +int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); """ ) # we need to "activate" the actual function first -c_cern_polygamma = cernlibc.c_cern_polygamma c_getdouble = cernlibc.c_getdouble c_ad_us_gamma_ns = cernlibc.c_ad_us_gamma_ns +c_ad_us_gamma_singlet = cernlibc.c_ad_us_gamma_singlet # allocate the pointers and get their addresses -re_y = ffi.new("double*") -im_y = ffi.new("double*") -re_y_address = int(ffi.cast("uintptr_t", re_y)) -im_y_address = int(ffi.cast("uintptr_t", im_y)) +re_double_4 = ffi.new("double[4]") +im_double_4 = ffi.new("double[4]") +re_double_4_address = int(ffi.cast("uintptr_t", re_double_4)) +im_double_4_address = int(ffi.cast("uintptr_t", im_double_4)) -re_double_3 = ffi.new("double[3]") -im_double_3 = ffi.new("double[3]") -re_double_3_address = int(ffi.cast("uintptr_t", re_double_3)) -im_double_3_address = int(ffi.cast("uintptr_t", im_double_3)) +re_double_16 = ffi.new("double[16]") +im_double_16 = ffi.new("double[16]") +re_double_16_address = int(ffi.cast("uintptr_t", re_double_16)) +im_double_16_address = int(ffi.cast("uintptr_t", im_double_16)) diff --git a/src/ekore/harmonics/polygamma.py b/src/ekore/harmonics/polygamma.py index a03620da2..a2c7e442f 100644 --- a/src/ekore/harmonics/polygamma.py +++ b/src/ekore/harmonics/polygamma.py @@ -6,21 +6,9 @@ import numba as nb import numpy as np -from ..cffilib import c_cern_polygamma, c_getdouble, im_y_address, re_y_address - - -@nb.njit() -def cern_polygamma(z, k): - """Compute polygamma from C.""" - # bypass complex and pointers - res = c_cern_polygamma(np.real(z), np.imag(z), k, re_y_address, im_y_address) - re_y = c_getdouble(re_y_address) - im_y = c_getdouble(im_y_address) - return re_y + im_y * 1j - @nb.njit(cache=True) -def _cern_polygamma(Z, K): # pylint: disable=all +def cern_polygamma(Z, K): # pylint: disable=all r"""Compute the polygamma functions :math:`\psi_k(z)`. Reimplementation of ``WPSIPG`` (C317) in `CERNlib `_ diff --git a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp index 47f5b33ab..e85eeadf4 100644 --- a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp @@ -1,5 +1,6 @@ #include "../../../constants.hpp" -#include "../../../harmonics/w1.hpp" +#include "../../../types.hpp" +#include "../../../harmonics/cache.hpp" namespace ekorepp { @@ -8,12 +9,71 @@ namespace unpolarized { namespace spacelike { namespace as1 { -std::complex gamma_ns(std::complex N, const unsigned int _nf) { - const std::complex S1 = harmonics::S1(N); - const std::complex gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0)); +/** + * @brief Compute the leading-order non-singlet anomalous dimension. + * Implements Eq. (3.4) of :cite:`Moch:2004pa`. + * @param c Mellin cache + * @param _nf Number of light flavors + * @return Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` + */ +cmplx gamma_ns(harmonics::Cache& c, const unsigned int _nf) { + const cmplx N = c.N(); + const cmplx S1 = c.get(harmonics::K::S1); + const cmplx gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0)); return CF * gamma; } +/** + * @brief Compute the leading-order quark-gluon anomalous dimension. + * Implements Eq. (3.5) of :cite:`Vogt:2004mw`. + * @param c Mellin cache + * @param nf Number of light flavors + * @return Leading-order quark-gluon anomalous dimension :math:`\\gamma_{qg}^{(0)}(N)` + */ +cmplx gamma_qg(harmonics::Cache& c, const unsigned int nf) { + const cmplx N = c.N(); + const cmplx gamma = -(pow(N, 2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0)); + return 2.0 * TR * 2.0 * nf * gamma; +} + +/** + * @brief Compute the leading-order gluon-quark anomalous dimension. + * Implements Eq. (3.5) of :cite:`Vogt:2004mw`. + * @param c Mellin cache + * @param _nf Number of light flavors + * @return Leading-order gluon-quark anomalous dimension :math:`\\gamma_{gq}^{(0)}(N)` + */ +cmplx gamma_gq(harmonics::Cache& c, const unsigned int _nf) { + const cmplx N = c.N(); + const cmplx gamma = -(pow(N, 2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0)); + return 2.0 * CF * gamma; +} + +/** + * @brief Compute the leading-order gluon-gluon anomalous dimension. + * Implements Eq. (3.5) of :cite:`Vogt:2004mw`. + * @param c Mellin cache + * @param nf Number of light flavors + * @return Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)` + */ +cmplx gamma_gg(harmonics::Cache& c, const unsigned int nf) { + const cmplx N = c.N(); + const cmplx S1 = c.get(harmonics::K::S1); + const cmplx gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0); + return CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * nf; +} + +/** + * @brief Compute the leading-order singlet anomalous dimension matrix. + * @param c Mellin cache + * @param nf Number of light flavors + * @return Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(0)}(N)` + */ +mtrx gamma_singlet(harmonics::Cache& c, const unsigned int nf) { + const cmplx gamma_qq = gamma_ns(c, nf); + return {{gamma_qq, gamma_qg(c, nf)}, {gamma_gq(c, nf), gamma_gg(c, nf)}}; +} + } // namepsace as1 } // namespace spacelike } // namespace unpolarized diff --git a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp index 6beeae806..77134747b 100644 --- a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp @@ -1,6 +1,6 @@ -#include +#include "../../../types.hpp" +#include "../../../harmonics/cache.hpp" #include "./as1.hpp" -#include "../../../harmonics/polygamma.hpp" namespace ekorepp { @@ -10,18 +10,31 @@ namespace spacelike { /** * @brief Compute the tower of the non-singlet anomalous dimensions. - * @param order_qcd QCD perturbative orders + * @param order_qcd QCD perturbative order * @param mode sector identifier - * @param n Mellin variable + * @param c Mellin cache * @param nf number of active flavors - * @return non-singlet anomalous dimensions + * @return Unpolarized, space-like, non-singlet anomalous dimensions */ -std::vector> gamma_ns(const unsigned int order_qcd, const unsigned int mode, std::complex n, const unsigned int nf) { - std::vector> gamma_ns(order_qcd); - gamma_ns[0] = as1::gamma_ns(n, nf); +vctr gamma_ns(const unsigned int order_qcd, const unsigned int mode, harmonics::Cache& c, const unsigned int nf) { + vctr gamma_ns(order_qcd); + gamma_ns[0] = as1::gamma_ns(c, nf); return gamma_ns; } +/** + * @brief Compute the tower of the singlet anomalous dimension matrices. + * @param order_qcd QCD perturbative order + * @param c Mellin cache + * @param nf Number of light flavors + * @return Unpolarized, space-like, singlet anomalous dimension matrices + */ +tnsr gamma_singlet(const unsigned int order_qcd, harmonics::Cache& c, const unsigned int nf) { + tnsr gamma_S(order_qcd); + gamma_S[0] = as1::gamma_singlet(c, nf); + return gamma_S; +} + } // namespace spacelike } // namespace unpolarized } // namespace anomalous_dimensions diff --git a/src/ekorepp/cinterface.cpp b/src/ekorepp/cinterface.cpp index 9a243ce67..c328430d6 100644 --- a/src/ekorepp/cinterface.cpp +++ b/src/ekorepp/cinterface.cpp @@ -1,37 +1,45 @@ // g++ -shared -o ekorepp.so -fPIC src/ekorepp/cinterface.cpp -#include #include -#include -#include -#include "./harmonics/polygamma.hpp" +#include "./types.hpp" +#include "./harmonics/cache.hpp" #include "./anomalous_dimensions/unpolarized/space_like/init.hpp" extern "C" double c_getdouble(double* ptr) { return *ptr; } -extern "C" int c_cern_polygamma(const double re_in, const double im_in, const unsigned int K, double* re_out, double* im_out) { - const std::complex Z (re_in, im_in); +extern "C" int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, + const double re_in, const double im_in, const unsigned int nf, double* re_out, double* im_out) { + const ekorepp::cmplx n (re_in, im_in); try { - const std::complex H = ekorepp::harmonics::cern_polygamma(Z, K); - *re_out = H.real(); - *im_out = H.imag(); + ekorepp::harmonics::Cache c(n); + const ekorepp::vctr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, c, nf); + for (unsigned int j = 0; j < res.size(); ++j) { + re_out[j] = res[j].real(); + im_out[j] = res[j].imag(); + } } catch (std::invalid_argument &e) { return -1; } return 0; } - -extern "C" int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, - const double re_in, const double im_in, const unsigned int nf, double* re_out, double* im_out) { - const std::complex n (re_in, im_in); +extern "C" int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, + double* re_out, double* im_out) { + const ekorepp::cmplx n (re_in, im_in); try { - const std::vector> res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, n, nf); - for (unsigned int j = 0; j < res.size(); ++j) { - re_out[j] = res[j].real(); - im_out[j] = res[j].imag(); + ekorepp::harmonics::Cache c(n); + const ekorepp::tnsr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_singlet(order_qcd, c, nf); + unsigned int j = 0; + for (unsigned int k = 0; k < res.size(); ++k) { + for (unsigned int l = 0; l < res[k].size(); ++l) { + for (unsigned int m = 0; m < res[k][l].size(); ++m) { + re_out[j] = res[k][l][m].real(); + im_out[j] = res[k][l][m].imag(); + ++j; + } + } } } catch (std::invalid_argument &e) { return -1; diff --git a/src/ekorepp/harmonics/cache.hpp b/src/ekorepp/harmonics/cache.hpp new file mode 100644 index 000000000..469e5c568 --- /dev/null +++ b/src/ekorepp/harmonics/cache.hpp @@ -0,0 +1,63 @@ +#ifndef EKOREPP_HARMONICS_CACHE_HPP_ +#define EKOREPP_HARMONICS_CACHE_HPP_ + +#include +#include "../types.hpp" +#include "./w1.hpp" + +namespace ekorepp { +namespace harmonics { + +/** @brief available elements */ +enum K { + S1 +}; + +/** @brief Cache type */ +typedef std::unordered_map mp; + +/** @brief A simple caching system */ +class Cache { + /** @brief Mellin base value */ + cmplx n; + + /** @brief cache list */ + mp m; + +public: + + /** + * @brief Init cache + * @param N Mellin base value + */ + Cache(const cmplx N) : n(N), m() {} + + /** + * Get Mellin base value + * @return Mellin base value + */ + cmplx N() const { return this->n; } + + /** + * @brief Obtain an element + * @param key key + * @return associated element + */ + cmplx get(const K key) { + mp::const_iterator it = this->m.find(key); + if (it != this->m.cend()) + return it->second; + cmplx res; + switch (key) + { + case S1: res = w1::S1(this->n); break; + } + this->m[key] = res; + return res; + } +}; + +} // namespace harmonics +} // namespace ekorepp + +#endif // EKOREPP_HARMONICS_CACHE_HPP_ diff --git a/src/ekorepp/harmonics/polygamma.hpp b/src/ekorepp/harmonics/polygamma.hpp index 7deca2af7..200154531 100644 --- a/src/ekorepp/harmonics/polygamma.hpp +++ b/src/ekorepp/harmonics/polygamma.hpp @@ -3,6 +3,7 @@ #include #include +#include "../types.hpp" namespace ekorepp { namespace harmonics { @@ -14,7 +15,7 @@ namespace harmonics { * @param K order of polygamma function * @return k-th polygamma function :math:`\psi_k(z)` */ -std::complex cern_polygamma(const std::complex Z, const unsigned int K) { +cmplx cern_polygamma(const cmplx Z, const unsigned int K) { const double DELTA = 5e-13; const double R1 = 1; const double HF = R1/2; @@ -64,7 +65,7 @@ std::complex cern_polygamma(const std::complex Z, const unsigned }, {10., -7., 12., -33., 130., -691.} }; - std::complex U=Z; + cmplx U=Z; double X=U.real(); double A=abs(X); if (K < 0 || K > 4) @@ -75,8 +76,8 @@ std::complex cern_polygamma(const std::complex Z, const unsigned const unsigned int K1=K+1; if (X < 0) U=-U; - std::complex V=U; - std::complex H=0; + cmplx V=U; + cmplx H=0; if (A < 15) { H=1./pow(V,K1); for (int i = 1; i < 14 - A_as_int + 1 ; ++i) { @@ -85,8 +86,8 @@ std::complex cern_polygamma(const std::complex Z, const unsigned } V=V+1.; } - std::complex R=1./pow(V,2); - std::complex P=R*C[K][6-1]; + cmplx R=1./pow(V,2); + cmplx P=R*C[K][6-1]; for (int i = 5; i>1-1; i--) P=R*(C[K][i-1]+P); H=SGN[K]*(FCT[K+1]*H+(V*(FCT[K-1+1]+P)+HF*FCT[K+1])/pow(V,K1)); @@ -99,7 +100,7 @@ std::complex cern_polygamma(const std::complex Z, const unsigned const double A=sin(X); const double B=cos(X); const double T=tanh(Y); - P=std::complex(B,-A*T)/std::complex(A,+B*T); + P=cmplx(B,-A*T)/cmplx(A,+B*T); if (0 == K) H=H+1./U+M_PI*P; else if (1 == K) @@ -117,8 +118,8 @@ std::complex cern_polygamma(const std::complex Z, const unsigned return H; } -std::complex recursive_harmonic_sum(const std::complex base_value, const std::complex n, const unsigned int iterations, const unsigned int weight) { - std::complex fact = 0.0; +cmplx recursive_harmonic_sum(const cmplx base_value, const cmplx n, const unsigned int iterations, const unsigned int weight) { + cmplx fact = 0.0; for (unsigned int i = 1; i <= iterations + 1; ++i) fact += pow(1.0 / (n + double(i)), weight); return base_value + fact; diff --git a/src/ekorepp/harmonics/w1.hpp b/src/ekorepp/harmonics/w1.hpp index 59c6f00f7..fce057d4d 100644 --- a/src/ekorepp/harmonics/w1.hpp +++ b/src/ekorepp/harmonics/w1.hpp @@ -1,8 +1,12 @@ -#include +#ifndef EKOREPP_HARMONICS_W1_HPP_ +#define EKOREPP_HARMONICS_W1_HPP_ + +#include "../types.hpp" #include "./polygamma.hpp" namespace ekorepp { namespace harmonics { +namespace w1 { /** * @brief Compute the harmonic sum :math:`S_1(N)`. @@ -14,9 +18,12 @@ namespace harmonics { * @param N Mellin moment * @return (simple) Harmonic sum :math:`S_1(N)` */ -std::complex S1(const std::complex N) { +cmplx S1(const cmplx N) { return cern_polygamma(N + 1.0, 0) + 0.5772156649015328606065120; } +} // namespace w1 } // namespace harmonics } // namespace ekorepp + +#endif // EKOREPP_HARMONICS_W1_HPP_ diff --git a/src/ekorepp/types.hpp b/src/ekorepp/types.hpp new file mode 100644 index 000000000..052d2adff --- /dev/null +++ b/src/ekorepp/types.hpp @@ -0,0 +1,16 @@ +#ifndef EKOREPP_TYPES_HPP_ +#define EKOREPP_TYPES_HPP_ + +#include +#include + +namespace ekorepp { + +typedef std::complex cmplx; +typedef std::vector vctr; +typedef std::vector mtrx; +typedef std::vector tnsr; + +} // namespace ekorepp + +#endif // EKOREPP_TYPES_HPP_ From 442edfb80e7885cc4668b425cc6aff3d60aeba68 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 11 Jul 2023 17:49:44 +0200 Subject: [PATCH 05/12] Start debugging ekorepp --- src/eko/evolution_operator/__init__.py | 9 +++-- .../unpolarized/space_like/__init__.py | 33 +++++------------ src/ekore/cffilib.py | 37 +++++++++++++------ 3 files changed, 40 insertions(+), 39 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index eab3c64de..c50ff126a 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -185,7 +185,7 @@ def integrand( return self.path.prefactor * pj * self.path.jac -@nb.njit(cache=True) +@nb.njit(cache=False) def quad_ker( u, order, @@ -313,7 +313,7 @@ def quad_ker( return np.real(ker * integrand) -@nb.njit(cache=True) +@nb.njit(cache=False) def quad_ker_qcd( ker_base, order, @@ -415,7 +415,7 @@ def quad_ker_qcd( if is_time_like: gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) else: - gamma_ns = ad_us.gamma_ns(order, mode0, ker_base.n, nf) + gamma_ns = ad_us.c_gamma_ns(order, mode0, ker_base.n, nf) if sv_mode == sv.Modes.exponentiated: gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) ker = ns.dispatcher( @@ -863,6 +863,7 @@ def run_op_integration( """ column = [] k, logx = log_grid + labels = self.labels start_time = time.perf_counter() # iterate basis functions for l, bf in enumerate(self.int_disp): @@ -870,7 +871,7 @@ def run_op_integration( continue temp_dict = {} # iterate sectors - for label in self.labels: + for label in labels: res = integrate.quad( self.quad_ker( label=label, logx=logx, areas=bf.areas_representation diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 57545495c..36eb9f017 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -23,11 +23,9 @@ from ....cffilib import ( c_ad_us_gamma_ns, c_ad_us_gamma_singlet, - c_getdouble, - im_double_4_address, - im_double_16_address, - re_double_4_address, - re_double_16_address, + im_double_address, + re_double_address, + read_complex, ) from ....harmonics import cache as c from . import aem1, aem2, as1, as1aem1, as2, as3, as4 @@ -43,20 +41,15 @@ def c_gamma_ns(order, mode, n, nf): np.real(n), np.imag(n), nf, - re_double_4_address, - im_double_4_address, + re_double_address, + im_double_address, ) - gamma_ns = np.zeros(order[0], np.complex_) - for j in range(order[0]): - gamma_ns[j] = ( - c_getdouble(re_double_4_address + j * 8) - + c_getdouble(im_double_4_address + j * 8) * 1j - ) + gamma_ns = read_complex(order[0]) return gamma_ns @nb.njit() -def c_gamma_singlet(order, n, nf): +def c_gamma_singlet(order, n, nf, _n3lo_ad_variation): """Compute gamma_singlet from C.""" # bypass complex and pointers res = c_ad_us_gamma_singlet( @@ -64,17 +57,11 @@ def c_gamma_singlet(order, n, nf): np.real(n), np.imag(n), nf, - re_double_16_address, - im_double_16_address, + re_double_address, + im_double_address, ) size = order[0] * 4 - gamma_s = np.zeros(size, np.complex_) - for j in range(size): - gamma_s[j] = ( - c_getdouble(re_double_16_address + j * 8) - + c_getdouble(im_double_16_address + j * 8) * 1j - ) - + gamma_s = read_complex(size) return gamma_s.reshape(order[0], 2, 2) diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index c6a031e35..e3a5f081a 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -1,8 +1,11 @@ """CFFI interface.""" +import numba as nb +import numpy as np +import numpy.typing as npt from cffi import FFI ffi = FFI() -cernlibc = ffi.dlopen("./ekorepp.so") +ekorepplibc = ffi.dlopen("./ekorepp.so") ffi.cdef( """ @@ -13,17 +16,27 @@ ) # we need to "activate" the actual function first -c_getdouble = cernlibc.c_getdouble -c_ad_us_gamma_ns = cernlibc.c_ad_us_gamma_ns -c_ad_us_gamma_singlet = cernlibc.c_ad_us_gamma_singlet +c_getdouble = ekorepplibc.c_getdouble +c_ad_us_gamma_ns = ekorepplibc.c_ad_us_gamma_ns +c_ad_us_gamma_singlet = ekorepplibc.c_ad_us_gamma_singlet # allocate the pointers and get their addresses -re_double_4 = ffi.new("double[4]") -im_double_4 = ffi.new("double[4]") -re_double_4_address = int(ffi.cast("uintptr_t", re_double_4)) -im_double_4_address = int(ffi.cast("uintptr_t", im_double_4)) +MAX_DOUBLES = 16 +_re_double = ffi.new(f"double[{MAX_DOUBLES}]") +_im_double = ffi.new(f"double[{MAX_DOUBLES}]") +re_double_address = int(ffi.cast("uintptr_t", _re_double)) +im_double_address = int(ffi.cast("uintptr_t", _im_double)) -re_double_16 = ffi.new("double[16]") -im_double_16 = ffi.new("double[16]") -re_double_16_address = int(ffi.cast("uintptr_t", re_double_16)) -im_double_16_address = int(ffi.cast("uintptr_t", im_double_16)) + +@nb.njit() +def read_complex(size: int) -> npt.ArrayLike: + """Read `size` complex numbers from C.""" + if size > MAX_DOUBLES: + raise MemoryError("Not enough memory allocated") + res = np.zeros(size, np.complex_) + for j in range(size): + res[j] = ( + c_getdouble(re_double_address + j * 8) + + c_getdouble(im_double_address + j * 8) * 1j + ) + return res From 6b5f31e53604a2a61dc3edb1fe47bcd7a804b9c7 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 12 Jul 2023 16:25:21 +0200 Subject: [PATCH 06/12] Fix polygamma --- src/ekorepp/cinterface.cpp | 12 ++++++++---- src/ekorepp/harmonics/polygamma.hpp | 17 ++++++++++++----- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/src/ekorepp/cinterface.cpp b/src/ekorepp/cinterface.cpp index c328430d6..fdf2785d1 100644 --- a/src/ekorepp/cinterface.cpp +++ b/src/ekorepp/cinterface.cpp @@ -19,8 +19,10 @@ extern "C" int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int re_out[j] = res[j].real(); im_out[j] = res[j].imag(); } - } catch (std::invalid_argument &e) { - return -1; + } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { + return -2; + } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { + return -3; } return 0; } @@ -41,8 +43,10 @@ extern "C" int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double } } } - } catch (std::invalid_argument &e) { - return -1; + } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { + return -2; + } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { + return -3; } return 0; } diff --git a/src/ekorepp/harmonics/polygamma.hpp b/src/ekorepp/harmonics/polygamma.hpp index 200154531..6d388d404 100644 --- a/src/ekorepp/harmonics/polygamma.hpp +++ b/src/ekorepp/harmonics/polygamma.hpp @@ -8,6 +8,12 @@ namespace ekorepp { namespace harmonics { +/** @brief Wrong order k of :math:`\psi_k(z)` */ +struct InvalidPolygammaOrder : public std::invalid_argument { using std::invalid_argument::invalid_argument; }; + +/** @brief Wrong argument z of :math:`\psi_k(z)` */ +struct InvalidPolygammaArgument : public std::invalid_argument { using std::invalid_argument::invalid_argument; }; + /** * @brief Compute the polygamma functions :math:`\psi_k(z)`. * Reimplementation of ``WPSIPG`` (C317) in `CERNlib `_ :cite:`KOLBIG1972221`. @@ -67,12 +73,13 @@ cmplx cern_polygamma(const cmplx Z, const unsigned int K) { }; cmplx U=Z; double X=U.real(); - double A=abs(X); + double A=fabs(X); if (K < 0 || K > 4) - throw std::invalid_argument("Order K has to be in [0:4]"); - const int A_as_int = A; - if (abs(U.imag()) < DELTA && abs(X+A_as_int) < DELTA) - throw std::invalid_argument("Argument Z equals non-positive integer"); + throw InvalidPolygammaOrder("Order K has to be in [0:4]"); + const int A_as_int = int(A); + if (fabs(U.imag()) < DELTA && fabs(X+A_as_int) < DELTA){ + throw InvalidPolygammaArgument("Argument Z equals non-positive integer"); + } const unsigned int K1=K+1; if (X < 0) U=-U; From 0616a81a8d352a97e8a178e48174b218a36afc79 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 12 Jul 2023 16:47:55 +0200 Subject: [PATCH 07/12] Reactivate c_gamma_singlet --- src/eko/evolution_operator/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index c50ff126a..2f9093b37 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -381,7 +381,7 @@ def quad_ker_qcd( if is_time_like: gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf) else: - gamma_singlet = ad_us.gamma_singlet( + gamma_singlet = ad_us.c_gamma_singlet( order, ker_base.n, nf, n3lo_ad_variation ) # scale var exponentiated is directly applied on gamma From 091072e65a171cf025143e50d1fdb921532edde9 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 12 Jul 2023 18:51:17 +0200 Subject: [PATCH 08/12] Call numba inside C --- src/eko/evolution_operator/__init__.py | 140 +++++++++++++++++++++---- src/ekore/cffilib.py | 5 + src/ekorepp/cinterface.cpp | 30 ++++++ 3 files changed, 157 insertions(+), 18 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 2f9093b37..160969879 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -16,6 +16,12 @@ import ekore.anomalous_dimensions.polarized.space_like as ad_ps import ekore.anomalous_dimensions.unpolarized.space_like as ad_us import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut +from ekore.cffilib import ( + c_quad_ker_qcd_ns, + im_double_address, + re_double_address, + read_complex, +) from .. import basis_rotation as br from .. import interpolation, mellin @@ -209,6 +215,7 @@ def quad_ker( n3lo_ad_variation, is_polarized, is_time_like, + quad_ker_qcd_ns_address, ): """Raw evolution kernel inside quad. @@ -288,6 +295,7 @@ def quad_ker( is_polarized, is_time_like, n3lo_ad_variation, + quad_ker_qcd_ns_address, ) else: ker = quad_ker_qed( @@ -331,6 +339,7 @@ def quad_ker_qcd( is_polarized, is_time_like, n3lo_ad_variation, + quad_ker_qcd_ns_address, ): """Raw evolution kernel inside quad. @@ -406,32 +415,126 @@ def quad_ker_qcd( ) @ np.ascontiguousarray(ker) ker = select_singlet_element(ker, mode0, mode1) else: - if is_polarized: - if is_time_like: - raise NotImplementedError("Polarized, time-like is not implemented") - else: - gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf) - else: - if is_time_like: - gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) - else: - gamma_ns = ad_us.c_gamma_ns(order, mode0, ker_base.n, nf) - if sv_mode == sv.Modes.exponentiated: - gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) - ker = ns.dispatcher( + # if is_polarized: + # if is_time_like: + # raise NotImplementedError("Polarized, time-like is not implemented") + # else: + # gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf) + # else: + # if is_time_like: + # gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) + # else: + # gamma_ns = ad_us.c_gamma_ns(order, mode0, ker_base.n, nf) + ker = quad_ker_qcd_ns( order, - method, - gamma_ns, + mode0, + ker_base.n, + nf, + quad_ker_qcd_ns_address, + L, as1, as0, - nf, ev_op_iterations, + is_threshold, ) - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker return ker +@nb.cfunc( + nb.types.void( + nb.types.CPointer(nb.types.double), + nb.types.CPointer(nb.types.double), + nb.types.uint, + nb.types.uint, + nb.types.double, + nb.types.uint, + nb.types.double, + nb.types.double, + nb.types.uint, + nb.types.uint, + nb.types.bool_, + nb.types.CPointer(nb.types.double), + nb.types.CPointer(nb.types.double), + ), + cache=True, +) +def cb_quad_ker_qcd_ns( + re_gamma_ns_raw, + im_gamma_ns_raw, + order_qcd, + nf, + L, + _method_num, + as1, + as0, + ev_op_iterations, + _sv_mode_num, + is_threshold, + re_ker, + im_ker, +): + """C Callback inside integration kernel.""" + re_gamma_ns = nb.carray(re_gamma_ns_raw, order_qcd) + im_gamma_ns = nb.carray(im_gamma_ns_raw, order_qcd) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + method = "iterate-expanded" + sv_mode = sv.Modes.exponentiated + order = (order_qcd, 0) + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) + ker = ns.dispatcher( + order, + method, + gamma_ns, + as1, + as0, + nf, + ev_op_iterations, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker + re_ker[0] = np.real(ker) + im_ker[0] = np.imag(ker) + + +@nb.njit() +def quad_ker_qcd_ns( + order, + mode, + n, + nf, + quad_ker_qcd_ns_address, + L, + as1, + as0, + ev_op_iterations, + is_threshold, +): + """Compute ker in C.""" + # bypass complex and pointers + method_num = 1 + sv_mode_num = 1 + res = c_quad_ker_qcd_ns( + order[0], + mode, + np.real(n), + np.imag(n), + nf, + quad_ker_qcd_ns_address, + L, + method_num, + as1, + as0, + ev_op_iterations, + sv_mode_num, + is_threshold, + re_double_address, + im_double_address, + ) + ker = read_complex(1) + return ker[0] + + @nb.njit(cache=True) def quad_ker_qed( ker_base, @@ -821,6 +924,7 @@ def quad_ker(self, label, logx, areas): n3lo_ad_variation=self.config["n3lo_ad_variation"], is_polarized=self.config["polarized"], is_time_like=self.config["time_like"], + quad_ker_qcd_ns_address=cb_quad_ker_qcd_ns.address, ) def initialize_op_members(self): diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index e3a5f081a..69b58584b 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -12,6 +12,10 @@ double c_getdouble(void* ptr); int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); +int c_quad_ker_qcd_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, + const unsigned int nf, void* py, const double L, + const unsigned int method_num, const double as1, const double as0, const unsigned int ev_op_iterations, + const unsigned int sv_mode_num, const bool is_threshold, void* re_out, void* im_out); """ ) @@ -19,6 +23,7 @@ c_getdouble = ekorepplibc.c_getdouble c_ad_us_gamma_ns = ekorepplibc.c_ad_us_gamma_ns c_ad_us_gamma_singlet = ekorepplibc.c_ad_us_gamma_singlet +c_quad_ker_qcd_ns = ekorepplibc.c_quad_ker_qcd_ns # allocate the pointers and get their addresses MAX_DOUBLES = 16 diff --git a/src/ekorepp/cinterface.cpp b/src/ekorepp/cinterface.cpp index fdf2785d1..203164c4f 100644 --- a/src/ekorepp/cinterface.cpp +++ b/src/ekorepp/cinterface.cpp @@ -1,6 +1,7 @@ // g++ -shared -o ekorepp.so -fPIC src/ekorepp/cinterface.cpp #include +#include #include "./types.hpp" #include "./harmonics/cache.hpp" #include "./anomalous_dimensions/unpolarized/space_like/init.hpp" @@ -50,3 +51,32 @@ extern "C" int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double } return 0; } + +typedef void (*py_quad_ker_qcd_ns)(double*, double*, + unsigned int, unsigned int, double, unsigned int, + double, double, + unsigned int, unsigned int, bool, double*, double*); + +extern "C" int c_quad_ker_qcd_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, + const unsigned int nf, py_quad_ker_qcd_ns py, const double L, + const unsigned int method_num, const double as1, const double as0, const unsigned int ev_op_iterations, + const unsigned int sv_mode_num, const bool is_threshold, double* re_out, double* im_out) { + const ekorepp::cmplx n (re_in, im_in); + try { + ekorepp::harmonics::Cache c(n); + const ekorepp::vctr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, c, nf); + std::vector re(order_qcd); + std::vector im(order_qcd); + for (unsigned int j = 0; j < res.size(); ++j) { + re[j] = res[j].real(); + im[j] = res[j].imag(); + } + py(re.data(), im.data(), order_qcd, nf, L, method_num, as1, as0, + ev_op_iterations, sv_mode_num, is_threshold, re_out, im_out); + } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { + return -2; + } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { + return -3; + } + return 0; +} From 3e5f8c81be2c9d594c33eeee7397b998e6665001 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 13 Jul 2023 19:05:36 +0200 Subject: [PATCH 09/12] Lift C dirctly under quad --- src/eko/evolution_operator/__init__.py | 231 ++++++++++++++---- src/ekore/cffilib.py | 17 +- .../{space_like/init.hpp => space_like.hpp} | 6 +- src/ekorepp/cinterface.cpp | 130 ++++++++-- 4 files changed, 310 insertions(+), 74 deletions(-) rename src/ekorepp/anomalous_dimensions/unpolarized/{space_like/init.hpp => space_like.hpp} (92%) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 160969879..05608dd90 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -3,6 +3,7 @@ See :doc:`Operator overview `. """ +import ctypes import functools import logging import os @@ -11,13 +12,13 @@ import numba as nb import numpy as np -from scipy import integrate +from scipy import LowLevelCallable, integrate import ekore.anomalous_dimensions.polarized.space_like as ad_ps import ekore.anomalous_dimensions.unpolarized.space_like as ad_us import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut from ekore.cffilib import ( - c_quad_ker_qcd_ns, + c_quad_ker_qcd, im_double_address, re_double_address, read_complex, @@ -440,61 +441,176 @@ def quad_ker_qcd( return ker +cb_quad_ker_qcd_T = ctypes.CFUNCTYPE( + ctypes.c_double, + ctypes.POINTER(ctypes.c_double), # re_gamma_ns_raw + ctypes.POINTER(ctypes.c_double), # im_gamma_ns_raw + ctypes.c_double, # re_n + ctypes.c_double, # im_n + ctypes.c_double, # re_jac + ctypes.c_double, # im_jac + ctypes.c_uint, # order_qcd + ctypes.c_bool, # is_singlet + ctypes.c_uint, # mode0 + ctypes.c_uint, # mode1 + ctypes.c_uint, # nf + ctypes.c_bool, # is_log + ctypes.c_double, # logx + ctypes.c_double * 12, # areas_raw + ctypes.c_uint, # polynomial_degreee + ctypes.c_double, # L + ctypes.c_uint, # method_num + ctypes.c_double, # as1 + ctypes.c_double, # as0 + ctypes.c_uint, # ev_op_iterations + ctypes.c_uint, # ev_op_max_order_qcd + ctypes.c_uint, # sv_mode_num + ctypes.c_bool, # is_threshold +) + +libc = ctypes.CDLL("./ekorepp.so") +c_quad_ker_qcd2 = libc.c_quad_ker_qcd + + +class QuadCargs(ctypes.Structure): + """Arguments to C call.""" + + _fields_ = [ + ("order_qcd", ctypes.c_uint), + ("mode0", ctypes.c_uint), + ("mode1", ctypes.c_uint), + ("is_polarized", ctypes.c_bool), + ("is_time_like", ctypes.c_bool), + ("nf", ctypes.c_uint), + ("py", cb_quad_ker_qcd_T), + ("is_log", ctypes.c_bool), + ("logx", ctypes.c_double), + ("areas", ctypes.POINTER(ctypes.c_double)), + ("areas_x", ctypes.c_uint), + ("areas_y", ctypes.c_uint), + ("L", ctypes.c_double), + ("method_num", ctypes.c_uint), + ("as1", ctypes.c_double), + ("as0", ctypes.c_double), + ("ev_op_iterations", ctypes.c_uint), + ("ev_op_max_order_qcd", ctypes.c_uint), + ("sv_mode_num", ctypes.c_uint), + ("is_threshold", ctypes.c_bool), + ] + + +c_quad_ker_qcd2.argtypes = [ctypes.c_double, ctypes.c_void_p] +c_quad_ker_qcd2.restype = ctypes.c_double + + @nb.cfunc( - nb.types.void( - nb.types.CPointer(nb.types.double), - nb.types.CPointer(nb.types.double), - nb.types.uint, - nb.types.uint, - nb.types.double, - nb.types.uint, - nb.types.double, - nb.types.double, - nb.types.uint, - nb.types.uint, - nb.types.bool_, - nb.types.CPointer(nb.types.double), - nb.types.CPointer(nb.types.double), + nb.types.double( + nb.types.CPointer(nb.types.double), # re_gamma_ns_raw + nb.types.CPointer(nb.types.double), # im_gamma_ns_raw + nb.types.double, # re_n + nb.types.double, # im_n + nb.types.double, # re_jac + nb.types.double, # im_jac + nb.types.uint, # order_qcd + nb.types.bool_, # is_singlet + nb.types.uint, # mode0 + nb.types.uint, # mode1 + nb.types.uint, # nf + nb.types.bool_, # is_log + nb.types.double, # logx + nb.types.CPointer(nb.types.double), # areas_raw + nb.types.uint, # areas_x + nb.types.uint, # areas_y + nb.types.double, # L + nb.types.uint, # method_num + nb.types.double, # as1 + nb.types.double, # as0 + nb.types.uint, # ev_op_iterations + nb.types.uint, # ev_op_max_order_qcd + nb.types.uint, # sv_mode_num + nb.types.bool_, # is_threshold ), - cache=True, + cache=False, ) -def cb_quad_ker_qcd_ns( - re_gamma_ns_raw, - im_gamma_ns_raw, +def cb_quad_ker_qcd( + re_gamma_raw, + im_gamma_raw, + re_n, + im_n, + re_jac, + im_jac, order_qcd, + is_singlet, + mode0, + mode1, nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, L, _method_num, as1, as0, ev_op_iterations, + ev_op_max_order_qcd, _sv_mode_num, is_threshold, - re_ker, - im_ker, ): """C Callback inside integration kernel.""" - re_gamma_ns = nb.carray(re_gamma_ns_raw, order_qcd) - im_gamma_ns = nb.carray(im_gamma_ns_raw, order_qcd) - gamma_ns = re_gamma_ns + im_gamma_ns * 1j + n = re_n + im_n * 1j + jac = re_jac + im_jac * 1j + areas = nb.carray(areas_raw, (areas_x, areas_y)) + pj = interpolation.evaluate_grid(n, is_log, logx, areas) method = "iterate-expanded" sv_mode = sv.Modes.exponentiated order = (order_qcd, 0) - if sv_mode == sv.Modes.exponentiated: - gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) - ker = ns.dispatcher( - order, - method, - gamma_ns, - as1, - as0, - nf, - ev_op_iterations, - ) - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker - re_ker[0] = np.real(ker) - im_ker[0] = np.imag(ker) + ev_op_max_order = (ev_op_max_order_qcd, 0) + if is_singlet: + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2)) + gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation( + gamma_singlet, order, nf, L + ) + ker = s.dispatcher( + order, + method, + gamma_singlet, + as1, + as0, + nf, + ev_op_iterations, + ev_op_max_order, + ) + # scale var expanded is applied on the kernel + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) + ) @ np.ascontiguousarray(ker) + ker = select_singlet_element(ker, mode0, mode1) + else: + re_gamma_ns = nb.carray(re_gamma_raw, order_qcd) + im_gamma_ns = nb.carray(im_gamma_raw, order_qcd) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) + ker = ns.dispatcher( + order, + method, + gamma_ns, + as1, + as0, + nf, + ev_op_iterations, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker + # recombine everything + res = ker * pj * jac + return np.real(res) @nb.njit() @@ -924,7 +1040,7 @@ def quad_ker(self, label, logx, areas): n3lo_ad_variation=self.config["n3lo_ad_variation"], is_polarized=self.config["polarized"], is_time_like=self.config["time_like"], - quad_ker_qcd_ns_address=cb_quad_ker_qcd_ns.address, + quad_ker_qcd_ns_address=cb_quad_ker_qcd.address, ) def initialize_op_members(self): @@ -970,16 +1086,43 @@ def run_op_integration( labels = self.labels start_time = time.perf_counter() # iterate basis functions + cfg = QuadCargs( + order_qcd=self.order[0], + is_polarized=self.config["polarized"], + is_time_like=self.config["time_like"], + nf=self.nf, + py=ctypes.cast(cb_quad_ker_qcd.address, cb_quad_ker_qcd_T), + is_log=self.int_disp.log, + logx=logx, + L=np.log(self.xif2), + method_num=1, + as1=self.as_list[1], + as0=self.as_list[0], + ev_op_iterations=self.config["ev_op_iterations"], + ev_op_max_order_qcd=self.config["ev_op_max_order"][0], + sv_mode_num=1, + is_threshold=self.is_threshold, + ) for l, bf in enumerate(self.int_disp): if k == l and l == self.grid_size - 1: continue + if bf.is_below_x(np.exp(logx)): + continue temp_dict = {} + curareas = bf.areas_representation + cfg.areas = (ctypes.c_double * (curareas.shape[0] * curareas.shape[1]))( + *curareas.flatten().tolist() + ) + cfg.areas_x = curareas.shape[0] + cfg.areas_y = curareas.shape[1] # iterate sectors for label in labels: + cfg.mode0 = label[0] + cfg.mode1 = label[1] + user_data = ctypes.cast(ctypes.pointer(cfg), ctypes.c_void_p) + func = LowLevelCallable(c_quad_ker_qcd2, user_data) res = integrate.quad( - self.quad_ker( - label=label, logx=logx, areas=bf.areas_representation - ), + func, 0.5, 1.0 - self._mellin_cut, epsabs=1e-12, diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py index 69b58584b..6f9d4c2cb 100644 --- a/src/ekore/cffilib.py +++ b/src/ekore/cffilib.py @@ -12,10 +12,17 @@ double c_getdouble(void* ptr); int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); -int c_quad_ker_qcd_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, - const unsigned int nf, void* py, const double L, - const unsigned int method_num, const double as1, const double as0, const unsigned int ev_op_iterations, - const unsigned int sv_mode_num, const bool is_threshold, void* re_out, void* im_out); +double c_quad_ker_qcd(const double u, + const unsigned int order_qcd, + const unsigned int mode0, const unsigned int mode1, + const bool is_polarized, const bool is_time_like, + const unsigned int nf, + void* py, + const bool is_log, const double logx, void* areas, const unsigned int polynomial_degree, + const double L, + const unsigned int method_num, const double as1, const double as0, + const unsigned int ev_op_iterations, const unsigned int ev_op_max_order, + const unsigned int sv_mode_num, const bool is_threshold); """ ) @@ -23,7 +30,7 @@ c_getdouble = ekorepplibc.c_getdouble c_ad_us_gamma_ns = ekorepplibc.c_ad_us_gamma_ns c_ad_us_gamma_singlet = ekorepplibc.c_ad_us_gamma_singlet -c_quad_ker_qcd_ns = ekorepplibc.c_quad_ker_qcd_ns +c_quad_ker_qcd = ekorepplibc.c_quad_ker_qcd # allocate the pointers and get their addresses MAX_DOUBLES = 16 diff --git a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp similarity index 92% rename from src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp rename to src/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp index 77134747b..31f5ff572 100644 --- a/src/ekorepp/anomalous_dimensions/unpolarized/space_like/init.hpp +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp @@ -1,6 +1,6 @@ -#include "../../../types.hpp" -#include "../../../harmonics/cache.hpp" -#include "./as1.hpp" +#include "../../types.hpp" +#include "../../harmonics/cache.hpp" +#include "./space_like/as1.hpp" namespace ekorepp { diff --git a/src/ekorepp/cinterface.cpp b/src/ekorepp/cinterface.cpp index 203164c4f..a65bee980 100644 --- a/src/ekorepp/cinterface.cpp +++ b/src/ekorepp/cinterface.cpp @@ -4,7 +4,7 @@ #include #include "./types.hpp" #include "./harmonics/cache.hpp" -#include "./anomalous_dimensions/unpolarized/space_like/init.hpp" +#include "./anomalous_dimensions/unpolarized/space_like.hpp" extern "C" double c_getdouble(double* ptr) { return *ptr; @@ -52,31 +52,117 @@ extern "C" int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double return 0; } -typedef void (*py_quad_ker_qcd_ns)(double*, double*, - unsigned int, unsigned int, double, unsigned int, - double, double, - unsigned int, unsigned int, bool, double*, double*); +class TalbotPath { + double t; + double r; + double o; +public: + TalbotPath(const double t, const double logx, const bool is_singlet) : t(t), r(0.4 * 16.0 / (1.0 - logx)), o(is_singlet ? 1. : 0.) {} + double theta() const { return M_PI * (2.0 * t - 1.0); } + ekorepp::cmplx n() const { + const double theta = this->theta(); + // treat singular point separately + const double re = 0.5 == t ? 1.0 : theta / tan(theta); + const double im = theta; + return o + r * ekorepp::cmplx(re, im); + } + ekorepp::cmplx jac() const { + const double theta = this->theta(); + // treat singular point separately + const double re = 0.5 == t ? 0.0 : 1.0 / tan(theta) - theta / pow(sin(theta), 2); + const double im = 1.0; + return r * M_PI * 2.0 * ekorepp::cmplx(re, im); + } + ekorepp::cmplx prefactor() const { + return ekorepp::cmplx(0,-1./M_PI); + } +}; -extern "C" int c_quad_ker_qcd_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, - const unsigned int nf, py_quad_ker_qcd_ns py, const double L, - const unsigned int method_num, const double as1, const double as0, const unsigned int ev_op_iterations, - const unsigned int sv_mode_num, const bool is_threshold, double* re_out, double* im_out) { - const ekorepp::cmplx n (re_in, im_in); - try { - ekorepp::harmonics::Cache c(n); - const ekorepp::vctr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, c, nf); - std::vector re(order_qcd); - std::vector im(order_qcd); +#include + +typedef double (*py_quad_ker_qcd)( + double* re_gamma, double* im_gamma, + double re_n, double im_n, + double re_jac, double im_jac, + unsigned int order_qcd, + bool is_singlet, + unsigned int mode0, + unsigned int mode1, + unsigned int nf, + bool is_log, + double logx, + double* areas, + unsigned int areas_x, + unsigned int areas_y, + double L, + unsigned int method_num, + double as1, double as0, + unsigned int ev_op_iterations, unsigned int ev_op_max_order_qcd, + unsigned int sv_mode_num, bool is_threshold); + +struct QuadCargs { + unsigned int order_qcd; + unsigned int mode0; + unsigned int mode1; + bool is_polarized; + bool is_time_like; + unsigned int nf; + py_quad_ker_qcd py; + bool is_log; + double logx; + double* areas; + unsigned int areas_x; + unsigned int areas_y; + double L; + unsigned int method_num; + double as1; + double as0; + unsigned int ev_op_iterations; + unsigned int ev_op_max_order_qcd; + unsigned int sv_mode_num; + bool is_threshold; +}; + +extern "C" double c_quad_ker_qcd(const double u, void* rargs) { + QuadCargs args = *(QuadCargs* )rargs; + const bool is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0); + // prepare gamma + const TalbotPath path(u, args.logx, is_singlet); + const ekorepp::cmplx jac = path.jac() * path.prefactor(); + ekorepp::harmonics::Cache c(path.n()); + const size_t size = args.order_qcd * (is_singlet ? 4 : 1); + std::vector re(size); + std::vector im(size); + if (is_singlet) { + const ekorepp::tnsr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_singlet(args.order_qcd, c, args.nf); + unsigned int j = 0; + for (unsigned int k = 0; k < res.size(); ++k) { + for (unsigned int l = 0; l < res[k].size(); ++l) { + for (unsigned int m = 0; m < res[k][l].size(); ++m) { + re[j] = res[k][l][m].real(); + im[j] = res[k][l][m].imag(); + ++j; + } + } + } + } else { + const ekorepp::vctr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(args.order_qcd, args.mode0, c, args.nf); for (unsigned int j = 0; j < res.size(); ++j) { re[j] = res[j].real(); im[j] = res[j].imag(); } - py(re.data(), im.data(), order_qcd, nf, L, method_num, as1, as0, - ev_op_iterations, sv_mode_num, is_threshold, re_out, im_out); - } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { - return -2; - } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { - return -3; } - return 0; + // pass on + return args.py(re.data(), im.data(), + c.N().real(), c.N().imag(), + jac.real(),jac.real(), + args.order_qcd, + is_singlet, + args.mode0, args.mode1, + args.nf, + args.is_log, args.logx, + args.areas, args.areas_x, args.areas_y, + args.L, args.method_num, args.as1, args.as0, + args.ev_op_iterations, args.ev_op_max_order_qcd, + args.sv_mode_num, args.is_threshold); } From 954bdd635d544333b9724ca5788e16b0d7fe69e5 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 14 Jul 2023 10:57:11 +0200 Subject: [PATCH 10/12] Reshuffle files --- src/eko/evolution_operator/__init__.py | 697 +----------------- .../evolution_operator/c_quad_ker.cpp} | 56 +- src/eko/evolution_operator/quad_ker.py | 426 +++++++++++ .../unpolarized/space_like/__init__.py | 41 -- src/ekore/cffilib.py | 54 -- 5 files changed, 432 insertions(+), 842 deletions(-) rename src/{ekorepp/cinterface.cpp => eko/evolution_operator/c_quad_ker.cpp} (68%) create mode 100644 src/eko/evolution_operator/quad_ker.py delete mode 100644 src/ekore/cffilib.py diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 05608dd90..9e1f6efbd 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -17,12 +17,6 @@ import ekore.anomalous_dimensions.polarized.space_like as ad_ps import ekore.anomalous_dimensions.unpolarized.space_like as ad_us import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut -from ekore.cffilib import ( - c_quad_ker_qcd, - im_double_address, - re_double_address, - read_complex, -) from .. import basis_rotation as br from .. import interpolation, mellin @@ -35,92 +29,10 @@ from ..kernels import valence_qed as qed_v from ..matchings import Segment from ..member import OpMember +from .quad_ker import QuadCargs, c_quad_ker_qcd, cb_quad_ker_qcd, cb_quad_ker_qcd_T logger = logging.getLogger(__name__) - -@nb.njit(cache=True) -def select_singlet_element(ker, mode0, mode1): - """Select element of the singlet matrix. - - Parameters - ---------- - ker : numpy.ndarray - singlet integration kernel - mode0 : int - id for first sector element - mode1 : int - id for second sector element - - Returns - ------- - complex - singlet integration kernel element - """ - k = 0 if mode0 == 100 else 1 - l = 0 if mode1 == 100 else 1 - return ker[k, l] - - -@nb.njit(cache=True) -def select_QEDsinglet_element(ker, mode0, mode1): - """Select element of the QEDsinglet matrix. - - Parameters - ---------- - ker : numpy.ndarray - QEDsinglet integration kernel - mode0 : int - id for first sector element - mode1 : int - id for second sector element - Returns - ------- - ker : complex - QEDsinglet integration kernel element - """ - if mode0 == 21: - index1 = 0 - elif mode0 == 22: - index1 = 1 - elif mode0 == 100: - index1 = 2 - else: - index1 = 3 - if mode1 == 21: - index2 = 0 - elif mode1 == 22: - index2 = 1 - elif mode1 == 100: - index2 = 2 - else: - index2 = 3 - return ker[index1, index2] - - -@nb.njit(cache=True) -def select_QEDvalence_element(ker, mode0, mode1): - """ - Select element of the QEDvalence matrix. - - Parameters - ---------- - ker : numpy.ndarray - QEDvalence integration kernel - mode0 : int - id for first sector element - mode1 : int - id for second sector element - Returns - ------- - ker : complex - QEDvalence integration kernel element - """ - index1 = 0 if mode0 == 10200 else 1 - index2 = 0 if mode1 == 10200 else 1 - return ker[index1, index2] - - spec = [ ("is_singlet", nb.boolean), ("is_QEDsinglet", nb.boolean), @@ -192,611 +104,6 @@ def integrand( return self.path.prefactor * pj * self.path.jac -@nb.njit(cache=False) -def quad_ker( - u, - order, - mode0, - mode1, - method, - is_log, - logx, - areas, - as_list, - mu2_from, - mu2_to, - a_half, - alphaem_running, - nf, - L, - ev_op_iterations, - ev_op_max_order, - sv_mode, - is_threshold, - n3lo_ad_variation, - is_polarized, - is_time_like, - quad_ker_qcd_ns_address, -): - """Raw evolution kernel inside quad. - - Parameters - ---------- - u : float - quad argument - order : int - perturbation order - mode0: int - pid for first sector element - mode1 : int - pid for second sector element - method : str - method - is_log : boolean - is a logarithmic interpolation - logx : float - Mellin inversion point - areas : tuple - basis function configuration - as1 : float - target coupling value - as0 : float - initial coupling value - mu2_from : float - initial value of mu2 - mu2_from : float - final value of mu2 - aem_list : list - list of electromagnetic coupling values - alphaem_running : bool - whether alphaem is running or not - nf : int - number of active flavors - L : float - logarithm of the squared ratio of factorization and renormalization scale - ev_op_iterations : int - number of evolution steps - ev_op_max_order : int - perturbative expansion order of U - sv_mode: int, `enum.IntEnum` - scale variation mode, see `eko.scale_variations.Modes` - is_threshold : boolean - is this an intermediate threshold operator? - n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` - is_polarized : boolean - is polarized evolution ? - is_time_like : boolean - is time-like evolution ? - - Returns - ------- - float - evaluated integration kernel - """ - ker_base = QuadKerBase(u, is_log, logx, mode0) - integrand = ker_base.integrand(areas) - if integrand == 0.0: - return 0.0 - if order[1] == 0: - ker = quad_ker_qcd( - ker_base, - order, - mode0, - mode1, - method, - as_list[-1], - as_list[0], - nf, - L, - ev_op_iterations, - ev_op_max_order, - sv_mode, - is_threshold, - is_polarized, - is_time_like, - n3lo_ad_variation, - quad_ker_qcd_ns_address, - ) - else: - ker = quad_ker_qed( - ker_base, - order, - mode0, - mode1, - method, - as_list, - mu2_from, - mu2_to, - a_half, - alphaem_running, - nf, - L, - ev_op_iterations, - ev_op_max_order, - sv_mode, - is_threshold, - ) - - # recombine everything - return np.real(ker * integrand) - - -@nb.njit(cache=False) -def quad_ker_qcd( - ker_base, - order, - mode0, - mode1, - method, - as1, - as0, - nf, - L, - ev_op_iterations, - ev_op_max_order, - sv_mode, - is_threshold, - is_polarized, - is_time_like, - n3lo_ad_variation, - quad_ker_qcd_ns_address, -): - """Raw evolution kernel inside quad. - - Parameters - ---------- - quad_ker : float - quad argument - order : int - perturbation order - mode0: int - pid for first sector element - mode1 : int - pid for second sector element - method : str - method - as1 : float - target coupling value - as0 : float - initial coupling value - nf : int - number of active flavors - L : float - logarithm of the squared ratio of factorization and renormalization scale - ev_op_iterations : int - number of evolution steps - ev_op_max_order : int - perturbative expansion order of U - sv_mode: int, `enum.IntEnum` - scale variation mode, see `eko.scale_variations.Modes` - is_threshold : boolean - is this an itermediate threshold operator? - n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` - - Returns - ------- - float - evaluated integration kernel - """ - # compute the actual evolution kernel for pure QCD - if ker_base.is_singlet: - if is_polarized: - if is_time_like: - raise NotImplementedError("Polarized, time-like is not implemented") - else: - gamma_singlet = ad_ps.gamma_singlet(order, ker_base.n, nf) - else: - if is_time_like: - gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf) - else: - gamma_singlet = ad_us.c_gamma_singlet( - order, ker_base.n, nf, n3lo_ad_variation - ) - # scale var exponentiated is directly applied on gamma - if sv_mode == sv.Modes.exponentiated: - gamma_singlet = sv.exponentiated.gamma_variation( - gamma_singlet, order, nf, L - ) - ker = s.dispatcher( - order, - method, - gamma_singlet, - as1, - as0, - nf, - ev_op_iterations, - ev_op_max_order, - ) - # scale var expanded is applied on the kernel - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = np.ascontiguousarray( - sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) - ) @ np.ascontiguousarray(ker) - ker = select_singlet_element(ker, mode0, mode1) - else: - # if is_polarized: - # if is_time_like: - # raise NotImplementedError("Polarized, time-like is not implemented") - # else: - # gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf) - # else: - # if is_time_like: - # gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) - # else: - # gamma_ns = ad_us.c_gamma_ns(order, mode0, ker_base.n, nf) - ker = quad_ker_qcd_ns( - order, - mode0, - ker_base.n, - nf, - quad_ker_qcd_ns_address, - L, - as1, - as0, - ev_op_iterations, - is_threshold, - ) - return ker - - -cb_quad_ker_qcd_T = ctypes.CFUNCTYPE( - ctypes.c_double, - ctypes.POINTER(ctypes.c_double), # re_gamma_ns_raw - ctypes.POINTER(ctypes.c_double), # im_gamma_ns_raw - ctypes.c_double, # re_n - ctypes.c_double, # im_n - ctypes.c_double, # re_jac - ctypes.c_double, # im_jac - ctypes.c_uint, # order_qcd - ctypes.c_bool, # is_singlet - ctypes.c_uint, # mode0 - ctypes.c_uint, # mode1 - ctypes.c_uint, # nf - ctypes.c_bool, # is_log - ctypes.c_double, # logx - ctypes.c_double * 12, # areas_raw - ctypes.c_uint, # polynomial_degreee - ctypes.c_double, # L - ctypes.c_uint, # method_num - ctypes.c_double, # as1 - ctypes.c_double, # as0 - ctypes.c_uint, # ev_op_iterations - ctypes.c_uint, # ev_op_max_order_qcd - ctypes.c_uint, # sv_mode_num - ctypes.c_bool, # is_threshold -) - -libc = ctypes.CDLL("./ekorepp.so") -c_quad_ker_qcd2 = libc.c_quad_ker_qcd - - -class QuadCargs(ctypes.Structure): - """Arguments to C call.""" - - _fields_ = [ - ("order_qcd", ctypes.c_uint), - ("mode0", ctypes.c_uint), - ("mode1", ctypes.c_uint), - ("is_polarized", ctypes.c_bool), - ("is_time_like", ctypes.c_bool), - ("nf", ctypes.c_uint), - ("py", cb_quad_ker_qcd_T), - ("is_log", ctypes.c_bool), - ("logx", ctypes.c_double), - ("areas", ctypes.POINTER(ctypes.c_double)), - ("areas_x", ctypes.c_uint), - ("areas_y", ctypes.c_uint), - ("L", ctypes.c_double), - ("method_num", ctypes.c_uint), - ("as1", ctypes.c_double), - ("as0", ctypes.c_double), - ("ev_op_iterations", ctypes.c_uint), - ("ev_op_max_order_qcd", ctypes.c_uint), - ("sv_mode_num", ctypes.c_uint), - ("is_threshold", ctypes.c_bool), - ] - - -c_quad_ker_qcd2.argtypes = [ctypes.c_double, ctypes.c_void_p] -c_quad_ker_qcd2.restype = ctypes.c_double - - -@nb.cfunc( - nb.types.double( - nb.types.CPointer(nb.types.double), # re_gamma_ns_raw - nb.types.CPointer(nb.types.double), # im_gamma_ns_raw - nb.types.double, # re_n - nb.types.double, # im_n - nb.types.double, # re_jac - nb.types.double, # im_jac - nb.types.uint, # order_qcd - nb.types.bool_, # is_singlet - nb.types.uint, # mode0 - nb.types.uint, # mode1 - nb.types.uint, # nf - nb.types.bool_, # is_log - nb.types.double, # logx - nb.types.CPointer(nb.types.double), # areas_raw - nb.types.uint, # areas_x - nb.types.uint, # areas_y - nb.types.double, # L - nb.types.uint, # method_num - nb.types.double, # as1 - nb.types.double, # as0 - nb.types.uint, # ev_op_iterations - nb.types.uint, # ev_op_max_order_qcd - nb.types.uint, # sv_mode_num - nb.types.bool_, # is_threshold - ), - cache=False, -) -def cb_quad_ker_qcd( - re_gamma_raw, - im_gamma_raw, - re_n, - im_n, - re_jac, - im_jac, - order_qcd, - is_singlet, - mode0, - mode1, - nf, - is_log, - logx, - areas_raw, - areas_x, - areas_y, - L, - _method_num, - as1, - as0, - ev_op_iterations, - ev_op_max_order_qcd, - _sv_mode_num, - is_threshold, -): - """C Callback inside integration kernel.""" - n = re_n + im_n * 1j - jac = re_jac + im_jac * 1j - areas = nb.carray(areas_raw, (areas_x, areas_y)) - pj = interpolation.evaluate_grid(n, is_log, logx, areas) - method = "iterate-expanded" - sv_mode = sv.Modes.exponentiated - order = (order_qcd, 0) - ev_op_max_order = (ev_op_max_order_qcd, 0) - if is_singlet: - re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2)) - im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2)) - gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j - if sv_mode == sv.Modes.exponentiated: - gamma_singlet = sv.exponentiated.gamma_variation( - gamma_singlet, order, nf, L - ) - ker = s.dispatcher( - order, - method, - gamma_singlet, - as1, - as0, - nf, - ev_op_iterations, - ev_op_max_order, - ) - # scale var expanded is applied on the kernel - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = np.ascontiguousarray( - sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) - ) @ np.ascontiguousarray(ker) - ker = select_singlet_element(ker, mode0, mode1) - else: - re_gamma_ns = nb.carray(re_gamma_raw, order_qcd) - im_gamma_ns = nb.carray(im_gamma_raw, order_qcd) - gamma_ns = re_gamma_ns + im_gamma_ns * 1j - if sv_mode == sv.Modes.exponentiated: - gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) - ker = ns.dispatcher( - order, - method, - gamma_ns, - as1, - as0, - nf, - ev_op_iterations, - ) - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker - # recombine everything - res = ker * pj * jac - return np.real(res) - - -@nb.njit() -def quad_ker_qcd_ns( - order, - mode, - n, - nf, - quad_ker_qcd_ns_address, - L, - as1, - as0, - ev_op_iterations, - is_threshold, -): - """Compute ker in C.""" - # bypass complex and pointers - method_num = 1 - sv_mode_num = 1 - res = c_quad_ker_qcd_ns( - order[0], - mode, - np.real(n), - np.imag(n), - nf, - quad_ker_qcd_ns_address, - L, - method_num, - as1, - as0, - ev_op_iterations, - sv_mode_num, - is_threshold, - re_double_address, - im_double_address, - ) - ker = read_complex(1) - return ker[0] - - -@nb.njit(cache=True) -def quad_ker_qed( - ker_base, - order, - mode0, - mode1, - method, - as_list, - mu2_from, - mu2_to, - a_half, - alphaem_running, - nf, - L, - ev_op_iterations, - ev_op_max_order, - sv_mode, - is_threshold, -): - """Raw evolution kernel inside quad. - - Parameters - ---------- - ker_base : QuadKerBase - quad argument - order : int - perturbation order - mode0: int - pid for first sector element - mode1 : int - pid for second sector element - method : str - method - as1 : float - target coupling value - as0 : float - initial coupling value - mu2_from : float - initial value of mu2 - mu2_from : float - final value of mu2 - aem_list : list - list of electromagnetic coupling values - alphaem_running : bool - whether alphaem is running or not - nf : int - number of active flavors - L : float - logarithm of the squared ratio of factorization and renormalization scale - ev_op_iterations : int - number of evolution steps - ev_op_max_order : int - perturbative expansion order of U - sv_mode: int, `enum.IntEnum` - scale variation mode, see `eko.scale_variations.Modes` - is_threshold : boolean - is this an itermediate threshold operator? - - Returns - ------- - float - evaluated integration kernel - """ - # compute the actual evolution kernel for QEDxQCD - if ker_base.is_QEDsinglet: - gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf) - # scale var exponentiated is directly applied on gamma - if sv_mode == sv.Modes.exponentiated: - gamma_s = sv.exponentiated.gamma_variation_qed( - gamma_s, order, nf, L, alphaem_running - ) - ker = qed_s.dispatcher( - order, - method, - gamma_s, - as_list, - a_half, - nf, - ev_op_iterations, - ev_op_max_order, - ) - # scale var expanded is applied on the kernel - # TODO : in this way a_half[-1][1] is the aem value computed in - # the middle point of the last step. Instead we want aem computed in mu2_final. - # However the distance between the two is very small and affects only the running aem - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = np.ascontiguousarray( - sv.expanded.singlet_variation_qed( - gamma_s, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L - ) - ) @ np.ascontiguousarray(ker) - ker = select_QEDsinglet_element(ker, mode0, mode1) - elif ker_base.is_QEDvalence: - gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) - # scale var exponentiated is directly applied on gamma - if sv_mode == sv.Modes.exponentiated: - gamma_v = sv.exponentiated.gamma_variation_qed( - gamma_v, order, nf, L, alphaem_running - ) - ker = qed_v.dispatcher( - order, - method, - gamma_v, - as_list, - a_half, - nf, - ev_op_iterations, - ev_op_max_order, - ) - # scale var expanded is applied on the kernel - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = np.ascontiguousarray( - sv.expanded.valence_variation_qed( - gamma_v, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L - ) - ) @ np.ascontiguousarray(ker) - ker = select_QEDvalence_element(ker, mode0, mode1) - else: - gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) - # scale var exponentiated is directly applied on gamma - if sv_mode == sv.Modes.exponentiated: - gamma_ns = sv.exponentiated.gamma_variation_qed( - gamma_ns, order, nf, L, alphaem_running - ) - ker = qed_ns.dispatcher( - order, - method, - gamma_ns, - as_list, - a_half[:, 1], - alphaem_running, - nf, - ev_op_iterations, - mu2_from, - mu2_to, - ) - if sv_mode == sv.Modes.expanded and not is_threshold: - ker = ( - sv.expanded.non_singlet_variation_qed( - gamma_ns, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L - ) - * ker - ) - return ker - - class Operator(sv.ModeMixin): """Internal representation of a single EKO. @@ -1120,7 +427,7 @@ def run_op_integration( cfg.mode0 = label[0] cfg.mode1 = label[1] user_data = ctypes.cast(ctypes.pointer(cfg), ctypes.c_void_p) - func = LowLevelCallable(c_quad_ker_qcd2, user_data) + func = LowLevelCallable(c_quad_ker_qcd, user_data) res = integrate.quad( func, 0.5, diff --git a/src/ekorepp/cinterface.cpp b/src/eko/evolution_operator/c_quad_ker.cpp similarity index 68% rename from src/ekorepp/cinterface.cpp rename to src/eko/evolution_operator/c_quad_ker.cpp index a65bee980..4519df156 100644 --- a/src/ekorepp/cinterface.cpp +++ b/src/eko/evolution_operator/c_quad_ker.cpp @@ -1,56 +1,10 @@ -// g++ -shared -o ekorepp.so -fPIC src/ekorepp/cinterface.cpp +// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp #include #include -#include "./types.hpp" -#include "./harmonics/cache.hpp" -#include "./anomalous_dimensions/unpolarized/space_like.hpp" - -extern "C" double c_getdouble(double* ptr) { - return *ptr; -} - -extern "C" int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, - const double re_in, const double im_in, const unsigned int nf, double* re_out, double* im_out) { - const ekorepp::cmplx n (re_in, im_in); - try { - ekorepp::harmonics::Cache c(n); - const ekorepp::vctr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_ns(order_qcd, mode, c, nf); - for (unsigned int j = 0; j < res.size(); ++j) { - re_out[j] = res[j].real(); - im_out[j] = res[j].imag(); - } - } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { - return -2; - } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { - return -3; - } - return 0; -} - -extern "C" int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, - double* re_out, double* im_out) { - const ekorepp::cmplx n (re_in, im_in); - try { - ekorepp::harmonics::Cache c(n); - const ekorepp::tnsr res = ekorepp::anomalous_dimensions::unpolarized::spacelike::gamma_singlet(order_qcd, c, nf); - unsigned int j = 0; - for (unsigned int k = 0; k < res.size(); ++k) { - for (unsigned int l = 0; l < res[k].size(); ++l) { - for (unsigned int m = 0; m < res[k][l].size(); ++m) { - re_out[j] = res[k][l][m].real(); - im_out[j] = res[k][l][m].imag(); - ++j; - } - } - } - } catch (ekorepp::harmonics::InvalidPolygammaOrder &e) { - return -2; - } catch (ekorepp::harmonics::InvalidPolygammaArgument &e) { - return -3; - } - return 0; -} +#include "./../../ekorepp/types.hpp" +#include "./../../ekorepp/harmonics/cache.hpp" +#include "./../../ekorepp/anomalous_dimensions/unpolarized/space_like.hpp" class TalbotPath { double t; @@ -78,8 +32,6 @@ class TalbotPath { } }; -#include - typedef double (*py_quad_ker_qcd)( double* re_gamma, double* im_gamma, double re_n, double im_n, diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py new file mode 100644 index 000000000..599934faf --- /dev/null +++ b/src/eko/evolution_operator/quad_ker.py @@ -0,0 +1,426 @@ +r"""Integration kernels.""" + +import ctypes +import logging + +import numba as nb +import numpy as np +from scipy import LowLevelCallable, integrate + +import ekore.anomalous_dimensions.polarized.space_like as ad_ps +import ekore.anomalous_dimensions.unpolarized.space_like as ad_us +import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut + +from .. import basis_rotation as br +from .. import interpolation, mellin +from .. import scale_variations as sv +from ..kernels import non_singlet as ns +from ..kernels import non_singlet_qed as qed_ns +from ..kernels import singlet as s +from ..kernels import singlet_qed as qed_s +from ..kernels import utils +from ..kernels import valence_qed as qed_v +from ..matchings import Segment +from ..member import OpMember + +logger = logging.getLogger(__name__) + + +@nb.njit(cache=True) +def select_singlet_element(ker, mode0, mode1): + """Select element of the singlet matrix. + + Parameters + ---------- + ker : numpy.ndarray + singlet integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + + Returns + ------- + complex + singlet integration kernel element + """ + k = 0 if mode0 == 100 else 1 + l = 0 if mode1 == 100 else 1 + return ker[k, l] + + +@nb.njit(cache=True) +def select_QEDsinglet_element(ker, mode0, mode1): + """Select element of the QEDsinglet matrix. + + Parameters + ---------- + ker : numpy.ndarray + QEDsinglet integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + Returns + ------- + ker : complex + QEDsinglet integration kernel element + """ + if mode0 == 21: + index1 = 0 + elif mode0 == 22: + index1 = 1 + elif mode0 == 100: + index1 = 2 + else: + index1 = 3 + if mode1 == 21: + index2 = 0 + elif mode1 == 22: + index2 = 1 + elif mode1 == 100: + index2 = 2 + else: + index2 = 3 + return ker[index1, index2] + + +@nb.njit(cache=True) +def select_QEDvalence_element(ker, mode0, mode1): + """ + Select element of the QEDvalence matrix. + + Parameters + ---------- + ker : numpy.ndarray + QEDvalence integration kernel + mode0 : int + id for first sector element + mode1 : int + id for second sector element + Returns + ------- + ker : complex + QEDvalence integration kernel element + """ + index1 = 0 if mode0 == 10200 else 1 + index2 = 0 if mode1 == 10200 else 1 + return ker[index1, index2] + + +cb_quad_ker_qcd_T = ctypes.CFUNCTYPE( + ctypes.c_double, + ctypes.POINTER(ctypes.c_double), # re_gamma_ns_raw + ctypes.POINTER(ctypes.c_double), # im_gamma_ns_raw + ctypes.c_double, # re_n + ctypes.c_double, # im_n + ctypes.c_double, # re_jac + ctypes.c_double, # im_jac + ctypes.c_uint, # order_qcd + ctypes.c_bool, # is_singlet + ctypes.c_uint, # mode0 + ctypes.c_uint, # mode1 + ctypes.c_uint, # nf + ctypes.c_bool, # is_log + ctypes.c_double, # logx + ctypes.c_double * 12, # areas_raw + ctypes.c_uint, # polynomial_degreee + ctypes.c_double, # L + ctypes.c_uint, # method_num + ctypes.c_double, # as1 + ctypes.c_double, # as0 + ctypes.c_uint, # ev_op_iterations + ctypes.c_uint, # ev_op_max_order_qcd + ctypes.c_uint, # sv_mode_num + ctypes.c_bool, # is_threshold +) + +libc = ctypes.CDLL("./c_quad_ker.so") +c_quad_ker_qcd = libc.c_quad_ker_qcd + + +class QuadCargs(ctypes.Structure): + """Arguments to C call.""" + + _fields_ = [ + ("order_qcd", ctypes.c_uint), + ("mode0", ctypes.c_uint), + ("mode1", ctypes.c_uint), + ("is_polarized", ctypes.c_bool), + ("is_time_like", ctypes.c_bool), + ("nf", ctypes.c_uint), + ("py", cb_quad_ker_qcd_T), + ("is_log", ctypes.c_bool), + ("logx", ctypes.c_double), + ("areas", ctypes.POINTER(ctypes.c_double)), + ("areas_x", ctypes.c_uint), + ("areas_y", ctypes.c_uint), + ("L", ctypes.c_double), + ("method_num", ctypes.c_uint), + ("as1", ctypes.c_double), + ("as0", ctypes.c_double), + ("ev_op_iterations", ctypes.c_uint), + ("ev_op_max_order_qcd", ctypes.c_uint), + ("sv_mode_num", ctypes.c_uint), + ("is_threshold", ctypes.c_bool), + ] + + +c_quad_ker_qcd.argtypes = [ctypes.c_double, ctypes.c_void_p] +c_quad_ker_qcd.restype = ctypes.c_double + + +@nb.cfunc( + nb.types.double( + nb.types.CPointer(nb.types.double), # re_gamma_ns_raw + nb.types.CPointer(nb.types.double), # im_gamma_ns_raw + nb.types.double, # re_n + nb.types.double, # im_n + nb.types.double, # re_jac + nb.types.double, # im_jac + nb.types.uint, # order_qcd + nb.types.bool_, # is_singlet + nb.types.uint, # mode0 + nb.types.uint, # mode1 + nb.types.uint, # nf + nb.types.bool_, # is_log + nb.types.double, # logx + nb.types.CPointer(nb.types.double), # areas_raw + nb.types.uint, # areas_x + nb.types.uint, # areas_y + nb.types.double, # L + nb.types.uint, # method_num + nb.types.double, # as1 + nb.types.double, # as0 + nb.types.uint, # ev_op_iterations + nb.types.uint, # ev_op_max_order_qcd + nb.types.uint, # sv_mode_num + nb.types.bool_, # is_threshold + ), + cache=False, +) +def cb_quad_ker_qcd( + re_gamma_raw, + im_gamma_raw, + re_n, + im_n, + re_jac, + im_jac, + order_qcd, + is_singlet, + mode0, + mode1, + nf, + is_log, + logx, + areas_raw, + areas_x, + areas_y, + L, + _method_num, + as1, + as0, + ev_op_iterations, + ev_op_max_order_qcd, + _sv_mode_num, + is_threshold, +): + """C Callback inside integration kernel.""" + n = re_n + im_n * 1j + jac = re_jac + im_jac * 1j + areas = nb.carray(areas_raw, (areas_x, areas_y)) + pj = interpolation.evaluate_grid(n, is_log, logx, areas) + method = "iterate-expanded" + sv_mode = sv.Modes.exponentiated + order = (order_qcd, 0) + ev_op_max_order = (ev_op_max_order_qcd, 0) + if is_singlet: + re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2)) + im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2)) + gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation( + gamma_singlet, order, nf, L + ) + ker = s.dispatcher( + order, + method, + gamma_singlet, + as1, + as0, + nf, + ev_op_iterations, + ev_op_max_order, + ) + # scale var expanded is applied on the kernel + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = np.ascontiguousarray( + sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L, dim=2) + ) @ np.ascontiguousarray(ker) + ker = select_singlet_element(ker, mode0, mode1) + else: + re_gamma_ns = nb.carray(re_gamma_raw, order_qcd) + im_gamma_ns = nb.carray(im_gamma_raw, order_qcd) + gamma_ns = re_gamma_ns + im_gamma_ns * 1j + if sv_mode == sv.Modes.exponentiated: + gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) + ker = ns.dispatcher( + order, + method, + gamma_ns, + as1, + as0, + nf, + ev_op_iterations, + ) + if sv_mode == sv.Modes.expanded and not is_threshold: + ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker + # recombine everything + res = ker * pj * jac + return np.real(res) + + +# @nb.njit(cache=True) +# def quad_ker_qed( +# ker_base, +# order, +# mode0, +# mode1, +# method, +# as_list, +# mu2_from, +# mu2_to, +# a_half, +# alphaem_running, +# nf, +# L, +# ev_op_iterations, +# ev_op_max_order, +# sv_mode, +# is_threshold, +# ): +# """Raw evolution kernel inside quad. + +# Parameters +# ---------- +# ker_base : QuadKerBase +# quad argument +# order : int +# perturbation order +# mode0: int +# pid for first sector element +# mode1 : int +# pid for second sector element +# method : str +# method +# as1 : float +# target coupling value +# as0 : float +# initial coupling value +# mu2_from : float +# initial value of mu2 +# mu2_from : float +# final value of mu2 +# aem_list : list +# list of electromagnetic coupling values +# alphaem_running : bool +# whether alphaem is running or not +# nf : int +# number of active flavors +# L : float +# logarithm of the squared ratio of factorization and renormalization scale +# ev_op_iterations : int +# number of evolution steps +# ev_op_max_order : int +# perturbative expansion order of U +# sv_mode: int, `enum.IntEnum` +# scale variation mode, see `eko.scale_variations.Modes` +# is_threshold : boolean +# is this an itermediate threshold operator? + +# Returns +# ------- +# float +# evaluated integration kernel +# """ +# # compute the actual evolution kernel for QEDxQCD +# if ker_base.is_QEDsinglet: +# gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_s = sv.exponentiated.gamma_variation_qed( +# gamma_s, order, nf, L, alphaem_running +# ) +# ker = qed_s.dispatcher( +# order, +# method, +# gamma_s, +# as_list, +# a_half, +# nf, +# ev_op_iterations, +# ev_op_max_order, +# ) +# # scale var expanded is applied on the kernel +# # TODO : in this way a_half[-1][1] is the aem value computed in +# # the middle point of the last step. Instead we want aem computed in mu2_final. +# # However the distance between the two is very small and affects only the running aem +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = np.ascontiguousarray( +# sv.expanded.singlet_variation_qed( +# gamma_s, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# ) @ np.ascontiguousarray(ker) +# ker = select_QEDsinglet_element(ker, mode0, mode1) +# elif ker_base.is_QEDvalence: +# gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_v = sv.exponentiated.gamma_variation_qed( +# gamma_v, order, nf, L, alphaem_running +# ) +# ker = qed_v.dispatcher( +# order, +# method, +# gamma_v, +# as_list, +# a_half, +# nf, +# ev_op_iterations, +# ev_op_max_order, +# ) +# # scale var expanded is applied on the kernel +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = np.ascontiguousarray( +# sv.expanded.valence_variation_qed( +# gamma_v, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# ) @ np.ascontiguousarray(ker) +# ker = select_QEDvalence_element(ker, mode0, mode1) +# else: +# gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) +# # scale var exponentiated is directly applied on gamma +# if sv_mode == sv.Modes.exponentiated: +# gamma_ns = sv.exponentiated.gamma_variation_qed( +# gamma_ns, order, nf, L, alphaem_running +# ) +# ker = qed_ns.dispatcher( +# order, +# method, +# gamma_ns, +# as_list, +# a_half[:, 1], +# alphaem_running, +# nf, +# ev_op_iterations, +# mu2_from, +# mu2_to, +# ) +# if sv_mode == sv.Modes.expanded and not is_threshold: +# ker = ( +# sv.expanded.non_singlet_variation_qed( +# gamma_ns, as_list[-1], a_half[-1][1], alphaem_running, order, nf, L +# ) +# * ker +# ) +# return ker diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 36eb9f017..6a3c6cf06 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -20,51 +20,10 @@ from eko import constants -from ....cffilib import ( - c_ad_us_gamma_ns, - c_ad_us_gamma_singlet, - im_double_address, - re_double_address, - read_complex, -) from ....harmonics import cache as c from . import aem1, aem2, as1, as1aem1, as2, as3, as4 -@nb.njit() -def c_gamma_ns(order, mode, n, nf): - """Compute gamma_ns from C.""" - # bypass complex and pointers - res = c_ad_us_gamma_ns( - order[0], - mode, - np.real(n), - np.imag(n), - nf, - re_double_address, - im_double_address, - ) - gamma_ns = read_complex(order[0]) - return gamma_ns - - -@nb.njit() -def c_gamma_singlet(order, n, nf, _n3lo_ad_variation): - """Compute gamma_singlet from C.""" - # bypass complex and pointers - res = c_ad_us_gamma_singlet( - order[0], - np.real(n), - np.imag(n), - nf, - re_double_address, - im_double_address, - ) - size = order[0] * 4 - gamma_s = read_complex(size) - return gamma_s.reshape(order[0], 2, 2) - - @nb.njit(cache=True) def gamma_ns(order, mode, n, nf): r"""Compute the tower of the non-singlet anomalous dimensions. diff --git a/src/ekore/cffilib.py b/src/ekore/cffilib.py deleted file mode 100644 index 6f9d4c2cb..000000000 --- a/src/ekore/cffilib.py +++ /dev/null @@ -1,54 +0,0 @@ -"""CFFI interface.""" -import numba as nb -import numpy as np -import numpy.typing as npt -from cffi import FFI - -ffi = FFI() -ekorepplibc = ffi.dlopen("./ekorepp.so") - -ffi.cdef( - """ -double c_getdouble(void* ptr); -int c_ad_us_gamma_ns(const unsigned int order_qcd, const unsigned int mode, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); -int c_ad_us_gamma_singlet(const unsigned int order_qcd, const double re_in, const double im_in, const unsigned int nf, void* re_out, void* im_out); -double c_quad_ker_qcd(const double u, - const unsigned int order_qcd, - const unsigned int mode0, const unsigned int mode1, - const bool is_polarized, const bool is_time_like, - const unsigned int nf, - void* py, - const bool is_log, const double logx, void* areas, const unsigned int polynomial_degree, - const double L, - const unsigned int method_num, const double as1, const double as0, - const unsigned int ev_op_iterations, const unsigned int ev_op_max_order, - const unsigned int sv_mode_num, const bool is_threshold); -""" -) - -# we need to "activate" the actual function first -c_getdouble = ekorepplibc.c_getdouble -c_ad_us_gamma_ns = ekorepplibc.c_ad_us_gamma_ns -c_ad_us_gamma_singlet = ekorepplibc.c_ad_us_gamma_singlet -c_quad_ker_qcd = ekorepplibc.c_quad_ker_qcd - -# allocate the pointers and get their addresses -MAX_DOUBLES = 16 -_re_double = ffi.new(f"double[{MAX_DOUBLES}]") -_im_double = ffi.new(f"double[{MAX_DOUBLES}]") -re_double_address = int(ffi.cast("uintptr_t", _re_double)) -im_double_address = int(ffi.cast("uintptr_t", _im_double)) - - -@nb.njit() -def read_complex(size: int) -> npt.ArrayLike: - """Read `size` complex numbers from C.""" - if size > MAX_DOUBLES: - raise MemoryError("Not enough memory allocated") - res = np.zeros(size, np.complex_) - for j in range(size): - res[j] = ( - c_getdouble(re_double_address + j * 8) - + c_getdouble(im_double_address + j * 8) * 1j - ) - return res From f7075339548be399d608f9b5e0327df964d071d9 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 14 Jul 2023 15:32:00 +0200 Subject: [PATCH 11/12] Fix 2 bugs --- src/eko/evolution_operator/__init__.py | 1 + src/eko/evolution_operator/c_quad_ker.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 9e1f6efbd..0205e3ff8 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -414,6 +414,7 @@ def run_op_integration( if k == l and l == self.grid_size - 1: continue if bf.is_below_x(np.exp(logx)): + column.append({label: (0.0, 0.0) for label in labels}) continue temp_dict = {} curareas = bf.areas_representation diff --git a/src/eko/evolution_operator/c_quad_ker.cpp b/src/eko/evolution_operator/c_quad_ker.cpp index 4519df156..fc3dc5aa1 100644 --- a/src/eko/evolution_operator/c_quad_ker.cpp +++ b/src/eko/evolution_operator/c_quad_ker.cpp @@ -107,7 +107,7 @@ extern "C" double c_quad_ker_qcd(const double u, void* rargs) { // pass on return args.py(re.data(), im.data(), c.N().real(), c.N().imag(), - jac.real(),jac.real(), + jac.real(),jac.imag(), args.order_qcd, is_singlet, args.mode0, args.mode1, From 765582b7858d11714db7b9250e69ec0e8ac5fb1c Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 14 Jul 2023 15:56:13 +0200 Subject: [PATCH 12/12] Update docs, remove temp hacks --- src/eko/evolution_operator/__init__.py | 60 ++++------------------- src/eko/evolution_operator/c_quad_ker.cpp | 37 ++++++++++++-- src/eko/evolution_operator/quad_ker.py | 13 +++-- 3 files changed, 52 insertions(+), 58 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 0205e3ff8..361cef5eb 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -29,7 +29,7 @@ from ..kernels import valence_qed as qed_v from ..matchings import Segment from ..member import OpMember -from .quad_ker import QuadCargs, c_quad_ker_qcd, cb_quad_ker_qcd, cb_quad_ker_qcd_T +from .quad_ker import QuadQCDargs, c_quad_ker_qcd, cb_quad_ker_qcd, cb_quad_ker_qcd_T logger = logging.getLogger(__name__) @@ -306,50 +306,6 @@ def labels(self): labels.extend(br.singlet_unified_labels) return labels - def quad_ker(self, label, logx, areas): - """Return partially initialized integrand function. - - Parameters - ---------- - label: tuple - operator element pids - logx: float - Mellin inversion point - areas : tuple - basis function configuration - - Returns - ------- - functools.partial - partially initialized integration kernel - - """ - return functools.partial( - quad_ker, - order=self.order, - mode0=label[0], - mode1=label[1], - method=self.config["method"], - is_log=self.int_disp.log, - logx=logx, - areas=areas, - as_list=self.as_list, - mu2_from=self.q2_from, - mu2_to=self.q2_to, - a_half=self.a_half_list, - alphaem_running=self.alphaem_running, - nf=self.nf, - L=np.log(self.xif2), - ev_op_iterations=self.config["ev_op_iterations"], - ev_op_max_order=tuple(self.config["ev_op_max_order"]), - sv_mode=self.sv_mode, - is_threshold=self.is_threshold, - n3lo_ad_variation=self.config["n3lo_ad_variation"], - is_polarized=self.config["polarized"], - is_time_like=self.config["time_like"], - quad_ker_qcd_ns_address=cb_quad_ker_qcd.address, - ) - def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( @@ -372,10 +328,7 @@ def initialize_op_members(self): else: self.op_members[n] = zero.copy() - def run_op_integration( - self, - log_grid, - ): + def run_op_integration(self, log_grid): """Run the integration for each grid point. Parameters @@ -390,10 +343,11 @@ def run_op_integration( """ column = [] k, logx = log_grid + # call(!) self.labels only once labels = self.labels start_time = time.perf_counter() - # iterate basis functions - cfg = QuadCargs( + # start preparing C arguments + cfg = QuadQCDargs( order_qcd=self.order[0], is_polarized=self.config["polarized"], is_time_like=self.config["time_like"], @@ -410,13 +364,16 @@ def run_op_integration( sv_mode_num=1, is_threshold=self.is_threshold, ) + # iterate basis functions for l, bf in enumerate(self.int_disp): if k == l and l == self.grid_size - 1: continue + # add emtpy labels with 0s if bf.is_below_x(np.exp(logx)): column.append({label: (0.0, 0.0) for label in labels}) continue temp_dict = {} + # prepare areas for C curareas = bf.areas_representation cfg.areas = (ctypes.c_double * (curareas.shape[0] * curareas.shape[1]))( *curareas.flatten().tolist() @@ -427,6 +384,7 @@ def run_op_integration( for label in labels: cfg.mode0 = label[0] cfg.mode1 = label[1] + # construct the low level object user_data = ctypes.cast(ctypes.pointer(cfg), ctypes.c_void_p) func = LowLevelCallable(c_quad_ker_qcd, user_data) res = integrate.quad( diff --git a/src/eko/evolution_operator/c_quad_ker.cpp b/src/eko/evolution_operator/c_quad_ker.cpp index fc3dc5aa1..de4cb47e7 100644 --- a/src/eko/evolution_operator/c_quad_ker.cpp +++ b/src/eko/evolution_operator/c_quad_ker.cpp @@ -1,4 +1,4 @@ -// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp +// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp -O2 #include #include @@ -6,13 +6,30 @@ #include "./../../ekorepp/harmonics/cache.hpp" #include "./../../ekorepp/anomalous_dimensions/unpolarized/space_like.hpp" +/** @brief Talbot path in Mellin inversion */ class TalbotPath { + /** @brief integration variable */ double t; + + /** @brief bending variable */ double r; + + /** @brief real offset */ double o; + + /** @brief auxilary angle */ + double theta() const { return M_PI * (2.0 * t - 1.0); } + public: + /** + * @brief Constructor from parameters + * @param t integration variable + * @param logx log of inversion point + * @param is_singlet add real offset? + */ TalbotPath(const double t, const double logx, const bool is_singlet) : t(t), r(0.4 * 16.0 / (1.0 - logx)), o(is_singlet ? 1. : 0.) {} - double theta() const { return M_PI * (2.0 * t - 1.0); } + + /** @brief Mellin-N */ ekorepp::cmplx n() const { const double theta = this->theta(); // treat singular point separately @@ -20,6 +37,8 @@ class TalbotPath { const double im = theta; return o + r * ekorepp::cmplx(re, im); } + + /** @brief transformation jacobian */ ekorepp::cmplx jac() const { const double theta = this->theta(); // treat singular point separately @@ -27,11 +46,14 @@ class TalbotPath { const double im = 1.0; return r * M_PI * 2.0 * ekorepp::cmplx(re, im); } + + /** @brief Mellin inversion prefactor */ ekorepp::cmplx prefactor() const { return ekorepp::cmplx(0,-1./M_PI); } }; +/** @brief Python callback signature */ typedef double (*py_quad_ker_qcd)( double* re_gamma, double* im_gamma, double re_n, double im_n, @@ -52,7 +74,8 @@ typedef double (*py_quad_ker_qcd)( unsigned int ev_op_iterations, unsigned int ev_op_max_order_qcd, unsigned int sv_mode_num, bool is_threshold); -struct QuadCargs { +/** @brief additional integration parameters */ +struct QuadQCDargs { unsigned int order_qcd; unsigned int mode0; unsigned int mode1; @@ -75,8 +98,14 @@ struct QuadCargs { bool is_threshold; }; +/** + * @brief intergration kernel inside quad + * @param u interation variable + * @param rargs raw args, to be casted to QuadQCDargs + * @return intergration value + */ extern "C" double c_quad_ker_qcd(const double u, void* rargs) { - QuadCargs args = *(QuadCargs* )rargs; + QuadQCDargs args = *(QuadQCDargs* )rargs; const bool is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0); // prepare gamma const TalbotPath path(u, args.logx, is_singlet); diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py index 599934faf..93011b131 100644 --- a/src/eko/evolution_operator/quad_ker.py +++ b/src/eko/evolution_operator/quad_ker.py @@ -123,7 +123,7 @@ def select_QEDvalence_element(ker, mode0, mode1): ctypes.c_uint, # nf ctypes.c_bool, # is_log ctypes.c_double, # logx - ctypes.c_double * 12, # areas_raw + ctypes.POINTER(ctypes.c_double), # areas_raw ctypes.c_uint, # polynomial_degreee ctypes.c_double, # L ctypes.c_uint, # method_num @@ -139,7 +139,7 @@ def select_QEDvalence_element(ker, mode0, mode1): c_quad_ker_qcd = libc.c_quad_ker_qcd -class QuadCargs(ctypes.Structure): +class QuadQCDargs(ctypes.Structure): """Arguments to C call.""" _fields_ = [ @@ -197,7 +197,7 @@ class QuadCargs(ctypes.Structure): nb.types.uint, # sv_mode_num nb.types.bool_, # is_threshold ), - cache=False, + cache=True, ) def cb_quad_ker_qcd( re_gamma_raw, @@ -226,15 +226,19 @@ def cb_quad_ker_qcd( is_threshold, ): """C Callback inside integration kernel.""" + # recover complex variables n = re_n + im_n * 1j jac = re_jac + im_jac * 1j + # combute basis functions areas = nb.carray(areas_raw, (areas_x, areas_y)) pj = interpolation.evaluate_grid(n, is_log, logx, areas) + # TODO recover parameters method = "iterate-expanded" sv_mode = sv.Modes.exponentiated order = (order_qcd, 0) ev_op_max_order = (ev_op_max_order_qcd, 0) if is_singlet: + # reconstruct singlet matrices re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2)) im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2)) gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j @@ -242,6 +246,7 @@ def cb_quad_ker_qcd( gamma_singlet = sv.exponentiated.gamma_variation( gamma_singlet, order, nf, L ) + # construct eko ker = s.dispatcher( order, method, @@ -259,11 +264,13 @@ def cb_quad_ker_qcd( ) @ np.ascontiguousarray(ker) ker = select_singlet_element(ker, mode0, mode1) else: + # construct non-singlet matrices re_gamma_ns = nb.carray(re_gamma_raw, order_qcd) im_gamma_ns = nb.carray(im_gamma_raw, order_qcd) gamma_ns = re_gamma_ns + im_gamma_ns * 1j if sv_mode == sv.Modes.exponentiated: gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) + # construct eko ker = ns.dispatcher( order, method,