Caching harmonic sums

benchmarks/eko/benchmark_evol_to_unity.py

@@ -11,7 +11,7 @@
 from eko.io.runcards import OperatorCard, TheoryCard
 from eko.runner.legacy import Runner
 
-# from eko.matching_conditions.operator_matrix_element import OperatorMatrixElement
+# from ekore.matching_conditions.operator_matrix_element import OperatorMatrixElement
 
 
 def update_cards(theory: TheoryCard, operator: OperatorCard):

benchmarks/eko/benchmark_ad.py → benchmarks/ekore/benchmark_ad.py

@@ -2,8 +2,8 @@
 import numpy as np
 import pytest
 
-import eko.anomalous_dimensions.as2 as ad_as2
-import eko.harmonics as h
+import ekore.anomalous_dimensions.unpolarized.space_like.as2 as ad_as2
+import ekore.harmonics as h
 from eko.constants import CA, CF, TR
 
 

benchmarks/performance/harmonics.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from eko import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.unpolarized.space_like as ad
 from eko.mellin import Path
 
 NF = 5

doc/source/overview/features.rst

@@ -1,8 +1,8 @@
 Features
 ========
 
-- perturbation orders: :mod:`LO <eko.anomalous_dimensions.lo>` + :mod:`NLO <eko.anomalous_dimensions.nlo>`
-  + :mod:`NNLO <eko.anomalous_dimensions.nnlo>`
+- perturbation orders: :mod:`LO <ekore.anomalous_dimensions.lo>` + :mod:`NLO <ekore.anomalous_dimensions.nlo>`
+  + :mod:`NNLO <ekore.anomalous_dimensions.nnlo>`
 - evolution in an (almost) arbitrary sequence of |FNS| (:class:`~eko.thresholds.ThresholdsAtlas`)
 - :ref:`theory/pQCD:scale variations`
 - different :doc:`solutions </theory/DGLAP>` of |DGLAP| differential equation

doc/source/theory/Mellin.rst

@@ -124,7 +124,7 @@ the Mellin inverse.
     S_m(N) = \sum\limits_{j=1}^N \frac 1 {j^m} = \frac{(-1)^{m-1}}{(m-1)!} \psi_{m-1}(N+1)+c_m \quad
     \text{with},\quad c_m = \left\{\begin{array}{ll} \gamma_E, & m=1\\ \zeta(m), & m>1\end{array} \right.
 
-  where :math:`\psi_k(N)` is the :math:`k`-th polygamma function (implemented as :meth:`~eko.harmonics.polygamma.cern_polygamma`)
+  where :math:`\psi_k(N)` is the :math:`k`-th polygamma function (implemented as :meth:`~ekore.harmonics.polygamma.cern_polygamma`)
   and :math:`\zeta` the Riemann zeta function (using :func:`scipy.special.zeta`).
 
 - for the sums :math:`S_{-m}(N)` and m > 0 we use :cite:`Gluck:1989ze`:
@@ -149,7 +149,7 @@ the Mellin inverse.
   which express higher weight sums in terms of simple one :math:`S_{m}, S_{-m}` and some irreducible integrals.
   The above prescription on the analytical continuation of :math:`\eta` is applied.
 
-The complete list of harmonics sums available in :mod:`eko.harmonics` is:
+The complete list of harmonics sums available in :mod:`ekore.harmonics` is:
 
     - weight 1:
 
@@ -166,20 +166,20 @@ The complete list of harmonics sums available in :mod:`eko.harmonics` is:
         .. math::
             S_{3}, S_{2,1}, S_{2,-1}, S_{-2,1}, S_{-2,-1}, S_{-3}
 
-        these sums relies on the integrals :mod:`eko.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
+        these sums relies on the integrals :mod:`ekore.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
 
     - weight 4:
 
         .. math ::
             S_{4}, S_{3,1}, S_{2,1,1}, S_{-2,-2}, S_{-3, 1}, S_{-4}
 
-        these sums relies on the integrals :mod:`eko.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
+        these sums relies on the integrals :mod:`ekore.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
 
     - weight 5:
 
         .. math ::
             S_{5}, S_{4,1}, S_{3,1,1}, S_{2,3}, S_{2,2,1}, S_{2,1,1,1}, S_{2,1,-2}, S_{2,-3}, S_{-2,3}, S_{-2,2,1}, S_{-2,1,1,1}, S_{-5}
 
-        these sums relies on the integrals :mod:`eko.harmonics.f_functions` :cite:`Blumlein:2009ta`
+        these sums relies on the integrals :mod:`ekore.harmonics.f_functions` :cite:`Blumlein:2009ta`
 
-We have also implemented a recursive computation of simple harmonics (single index), see :func:`eko.harmonics.polygamma.recursive_harmonic_sum`
+We have also implemented a recursive computation of simple harmonics (single index), see :func:`ekore.harmonics.polygamma.recursive_harmonic_sum`

pyproject.toml

@@ -133,6 +133,7 @@ python_functions = ['test_*', 'benchmark_*']
 addopts = [
   '--cov=eko',
   '--cov=ekobox',
+  '--cov=ekore',
   '--cov-report=html',
   '--cov-report=xml',
   '--strict-markers',

src/eko/beta.py

@@ -7,7 +7,7 @@
 import numba as nb
 
 from . import constants
-from .harmonics.constants import zeta3
+from ekore.harmonics.constants import zeta3
 
 
 @nb.njit(cache=True)

src/eko/evolution_operator/__init__.py

@@ -13,7 +13,11 @@
 import numpy as np
 from scipy import integrate
 
-from .. import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.polarized.space_like as ad_ps
+import ekore.anomalous_dimensions.polarized.time_like as ad_pt
+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
@@ -130,6 +134,8 @@ def quad_ker(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    is_polarized,
+    is_time_like,
 ):
     """Raw evolution kernel inside quad.
 
@@ -168,7 +174,11 @@ def quad_ker(
     sv_mode: int, `enum.IntEnum`
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
-        is this an itermediate threshold operator?
+        is this an intermediate threshold operator?
+    is_polarized : boolean
+        is polarized evolution ?
+    is_time_like : boolean
+        is time-like evolution ?
 
     Returns
     -------
@@ -182,7 +192,16 @@ def quad_ker(
 
     # compute the actual evolution kernel
     if ker_base.is_singlet:
-        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf)
+        if is_polarized:
+            if is_time_like:
+                gamma_singlet = ad_pt.gamma_singlet(order, ker_base.n, nf)
+            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)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_singlet = sv.exponentiated.gamma_variation(
@@ -205,7 +224,16 @@ def quad_ker(
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
-        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf)
+        if is_polarized:
+            if is_time_like:
+                gamma_ns = ad_pt.gamma_ns(order, mode0, ker_base.n, nf)
+            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(
@@ -382,6 +410,8 @@ def quad_ker(self, label, logx, areas):
             ev_op_max_order=tuple(self.config["ev_op_max_order"]),
             sv_mode=self.sv_mode,
             is_threshold=self.is_threshold,
+            is_polarized=self.config["polarized"],
+            is_time_like=self.config["time_like"],
         )
 
     def initialize_op_members(self):

src/eko/evolution_operator/grid.py

@@ -11,12 +11,12 @@
 import numpy as np
 import numpy.typing as npt
 
-from .. import matching_conditions, member
+from .. import member
 from .. import scale_variations as sv
 from ..io.runcards import Configs, Debug
-from ..matching_conditions.operator_matrix_element import OperatorMatrixElement
 from ..thresholds import flavor_shift, is_downward_path
-from . import Operator, flavors, physical
+from . import Operator, flavors, matching_condition, physical
+from .operator_matrix_element import OperatorMatrixElement
 
 logger = logging.getLogger(__name__)
 
@@ -89,6 +89,8 @@ def __init__(
         config["n_integration_cores"] = configs.n_integration_cores
         config["debug_skip_singlet"] = debug.skip_singlet
         config["debug_skip_non_singlet"] = debug.skip_non_singlet
+        config["polarized"] = configs.polarized
+        config["time_like"] = configs.time_like
 
         if method not in [
             "iterate-exact",
@@ -244,7 +246,7 @@ def generate(self, q2):
 
             # join with the basis rotation, since matching requires c+ (or likewise)
             if is_downward:
-                matching = matching_conditions.MatchingCondition.split_ad_to_evol_map(
+                matching = matching_condition.MatchingCondition.split_ad_to_evol_map(
                     self._matching_operators[op.q2_to],
                     op.nf - 1,
                     op.q2_to,
@@ -255,7 +257,7 @@ def generate(self, q2):
                 )
                 final_op = final_op @ matching @ invrot @ phys_op
             else:
-                matching = matching_conditions.MatchingCondition.split_ad_to_evol_map(
+                matching = matching_condition.MatchingCondition.split_ad_to_evol_map(
                     self._matching_operators[op.q2_to],
                     op.nf,
                     op.q2_to,