diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index eab3c64de..361cef5eb 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,7 +12,7 @@ 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 @@ -28,92 +29,10 @@ from ..kernels import valence_qed as qed_v from ..matchings import Segment from ..member import OpMember +from .quad_ker import QuadQCDargs, 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), @@ -185,399 +104,6 @@ def integrand( return self.path.prefactor * pj * self.path.jac -@nb.njit(cache=True) -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, -): - """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, - ) - 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=True) -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, -): - """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.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.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( - 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 - return ker - - -@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. @@ -780,49 +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"], - ) - def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( @@ -845,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 @@ -863,18 +343,52 @@ def run_op_integration( """ column = [] k, logx = log_grid + # call(!) self.labels only once + labels = self.labels start_time = time.perf_counter() + # start preparing C arguments + cfg = QuadQCDargs( + 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, + ) # 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() + ) + cfg.areas_x = curareas.shape[0] + cfg.areas_y = curareas.shape[1] # iterate sectors - for label in self.labels: + 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( - 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/eko/evolution_operator/c_quad_ker.cpp b/src/eko/evolution_operator/c_quad_ker.cpp new file mode 100644 index 000000000..de4cb47e7 --- /dev/null +++ b/src/eko/evolution_operator/c_quad_ker.cpp @@ -0,0 +1,149 @@ +// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp -O2 + +#include +#include +#include "./../../ekorepp/types.hpp" +#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.) {} + + /** @brief Mellin-N */ + 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); + } + + /** @brief transformation jacobian */ + 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); + } + + /** @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, + 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); + +/** @brief additional integration parameters */ +struct QuadQCDargs { + 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; +}; + +/** + * @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) { + 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); + 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(); + } + } + // pass on + return args.py(re.data(), im.data(), + c.N().real(), c.N().imag(), + jac.real(),jac.imag(), + 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); +} diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py new file mode 100644 index 000000000..93011b131 --- /dev/null +++ b/src/eko/evolution_operator/quad_ker.py @@ -0,0 +1,433 @@ +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.POINTER(ctypes.c_double), # 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 QuadQCDargs(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=True, +) +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.""" + # 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 + if sv_mode == sv.Modes.exponentiated: + gamma_singlet = sv.exponentiated.gamma_variation( + gamma_singlet, order, nf, L + ) + # construct eko + 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: + # 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, + 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/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp b/src/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp new file mode 100644 index 000000000..31f5ff572 --- /dev/null +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like.hpp @@ -0,0 +1,41 @@ +#include "../../types.hpp" +#include "../../harmonics/cache.hpp" +#include "./space_like/as1.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 order + * @param mode sector identifier + * @param c Mellin cache + * @param nf number of active flavors + * @return Unpolarized, space-like, non-singlet anomalous dimensions + */ +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 +} // namespace ekorepp 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..e85eeadf4 --- /dev/null +++ b/src/ekorepp/anomalous_dimensions/unpolarized/space_like/as1.hpp @@ -0,0 +1,81 @@ +#include "../../../constants.hpp" +#include "../../../types.hpp" +#include "../../../harmonics/cache.hpp" + + +namespace ekorepp { +namespace anomalous_dimensions { +namespace unpolarized { +namespace spacelike { +namespace as1 { + +/** + * @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 +} // namespace anomalous_dimensions +} // namespace ekorepp 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/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 new file mode 100644 index 000000000..6d388d404 --- /dev/null +++ b/src/ekorepp/harmonics/polygamma.hpp @@ -0,0 +1,138 @@ +#ifndef EKOREPP_HARMONICS_POLYGAMMA_HPP_ +#define EKOREPP_HARMONICS_POLYGAMMA_HPP_ + +#include +#include +#include "../types.hpp" + +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`. + * @param Z argument of polygamma function + * @param K order of polygamma function + * @return k-th polygamma function :math:`\psi_k(z)` + */ +cmplx cern_polygamma(const cmplx Z, const unsigned int K) { + 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.} + }; + cmplx U=Z; + double X=U.real(); + double A=fabs(X); + if (K < 0 || K > 4) + 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; + 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) { + V=V+1.; + H=H+1./pow(V,K1); + } + V=V+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)); + if (0 == K) + H=H+log(V); + if (X < 0){ + V=M_PI*U; + X=V.real(); + const double Y=V.imag(); + const double A=sin(X); + const double B=cos(X); + const double T=tanh(Y); + P=cmplx(B,-A*T)/cmplx(A,+B*T); + if (0 == K) + H=H+1./U+M_PI*P; + else if (1 == K) + H=-H+1./pow(U,2)+C1*(pow(P,2)+1.); + else if (2 == K) + H=H+2./pow(U,3)+C2*P*(pow(P,2)+1.); + else if (3 == K) { + R=pow(P,2); + H=-H+6./pow(U,4)+C3*((3.*R+4.)*R+1.); + } else if (4 == K) { + R=pow(P,2); + H=H+24./pow(U,5)+C4*P*((3.*R+5.)*R+2.); + } + } + return H; +} + +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; +} + +} // 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..fce057d4d --- /dev/null +++ b/src/ekorepp/harmonics/w1.hpp @@ -0,0 +1,29 @@ +#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)`. + * .. 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)` + */ +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_