Caching harmonic sums

doc/source/overview/features.rst

@@ -1,7 +1,12 @@
 Features
 ========
 
+<<<<<<< HEAD
+- perturbation orders: :mod:`LO <ekore.anomalous_dimensions.lo>` + :mod:`NLO <ekore.anomalous_dimensions.nlo>`
+  + :mod:`NNLO <ekore.anomalous_dimensions.nnlo>`
+=======
 - perturbation orders of anomalous dimensions: :math:`\alpha_s,\alpha_{em},\alpha_s\alpha_{em},\alpha_s^2,\alpha_s^3`
+>>>>>>> master
 - 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

src/eko/beta.py

@@ -7,8 +7,15 @@
 import numba as nb
 
 from . import constants
+
+<<<<<<< HEAD
+from eko.constants import zeta3
+
+=======
 from ekore.harmonics.constants import zeta3
 
+>>>>>>> master
+
 
 @nb.njit(cache=True)
 def beta_qcd_as2(nf):

src/eko/constants.py

@@ -1,6 +1,8 @@
 """This files sets the physical constants."""
 
 import numba as nb
+import numpy as np
+from scipy.special import zeta
 
 NC = 3
 """The number of colors."""
@@ -29,6 +31,23 @@
 ed2 = 1.0 / 9
 """Down quarks charge squared."""
 
+zeta2 = zeta(2)
+r""":math:`\zeta(2)`"""
+
+zeta3 = zeta(3)
+r""":math:`\zeta(3)`"""
+
+zeta4 = zeta(4)
+r""":math:`\zeta(4)`"""
+
+zeta5 = zeta(5)
+r""":math:`\zeta(5)`"""
+
+log2 = np.log(2)
+r""":math:`\ln(2)`"""
+
+li4half = 0.517479
+""":math:`Li_{4}(1/2)`"""
 
 def update_colors(nc):
     """Updates the number of colors to :math:`NC = nc`.

src/eko/evolution_operator/__init__.py

@@ -14,6 +14,12 @@
 from scipy import integrate
 
 import ekore.anomalous_dimensions.polarized.space_like as ad_ps
+
+<<<<<<< HEAD
+import ekore.anomalous_dimensions.polarized.time_like as ad_pt
+
+=======
+>>>>>>> master
 import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
 import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut
 
@@ -222,7 +228,11 @@ def quad_ker(
     else:
         if is_polarized:
             if is_time_like:
+<<<<<<< HEAD
+                gamma_ns = ad_pt.gamma_ns(order, mode0, ker_base.n, nf)
+=======
                 raise NotImplementedError("Polarized, time-like is not implemented")
+>>>>>>> master
             else:
                 gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf)
         else:

src/eko/evolution_operator/operator_matrix_element.py

@@ -36,7 +36,10 @@ def build_ome(A, matching_order, a_s, backward_method):
     -------
     ome : numpy.ndarray
         matching operator matrix
+<<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py
+=======
 
+>>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py
     """
     # to get the inverse one can use this FORM snippet
     # Symbol a;
@@ -123,7 +126,10 @@ def quad_ker(
     -------
     ker : float
         evaluated integration kernel
+<<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py
+=======
 
+>>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py
     """
     ker_base = QuadKerBase(u, is_log, logx, mode0)
     integrand = ker_base.integrand(areas)
@@ -158,7 +164,11 @@ def quad_ker(
         indices = {21: 0, 100: 1, 90: 2}
         if is_polarized:
             if is_time_like:
+<<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py
+                A = ome_pt.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
+=======
                 raise NotImplementedError("Polarized, time-like is not implemented")
+>>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py
             else:
                 A = ome_ps.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
         else:
@@ -170,12 +180,21 @@ def quad_ker(
         indices = {200: 0, 91: 1}
         if is_polarized:
             if is_time_like:
+<<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
+            else:
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
+        else:
+            if is_time_like:
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
+=======
                 raise NotImplementedError("Polarized, time-like is not implemented")
             else:
                 A = ome_ps.A_non_singlet(order, ker_base.n, sx, nf, L)
         else:
             if is_time_like:
                 A = ome_ut.A_non_singlet(order, ker_base.n, sx, nf, L)
+>>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py
             else:
                 A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
 

src/eko/gamma.py

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

src/ekore/anomalous_dimensions/polarized/time_like/__init__.py

@@ -0,0 +1,21 @@
+r"""The polarized, time-like Altarelli-Parisi splitting kernels.
+
+Normalization is given by
+
+.. math::
+    \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x)
+
+with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`.
+"""
+
+import numba as nb
+
+
+@nb.njit(cache=True)
+def gamma_ns(_order, _mode, _n, _nf):
+    raise NotImplementedError("Polarised, time-like is not yet implemented")
+
+
+@nb.njit(cache=True)
+def gamma_singlet(_order, _n, _nf):
+    raise NotImplementedError("Polarised, time-like is not yet implemented")

src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py

@@ -19,6 +19,12 @@
 import numpy as np
 
 from .... import harmonics
+
+<<<<<<< HEAD:src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+from ....harmonics import cache as c
+
+=======
+>>>>>>> master:src/eko/anomalous_dimensions/__init__.py
 from . import aem1, aem2, as1, as2, as3, as4
 
 
@@ -70,7 +76,9 @@ def gamma_ns(order, mode, n, nf):
         sx = harmonics.sx(n, max_weight=order[0] + 1)
     # now combine
     gamma_ns = np.zeros(order[0], np.complex_)
-    gamma_ns[0] = as1.gamma_ns(n, sx[0])
+
+    cache = c.reset()
+    gamma_ns[0] = as1.gamma_ns(n, cache)
     # NLO and beyond
     if order[0] >= 2:
         if mode == 10101:
@@ -145,8 +153,9 @@ def gamma_singlet(order, n, nf):
     else:
         sx = harmonics.sx(n, max_weight=order[0])
 
+    cache = c.reset()
     gamma_s = np.zeros((order[0], 2, 2), np.complex_)
-    gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
+    gamma_s[0] = as1.gamma_singlet(n, nf, cache)
     if order[0] >= 2:
         gamma_s[1] = as2.gamma_singlet(n, nf, sx)
     if order[0] >= 3:

src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py

@@ -5,9 +5,11 @@
 
 from eko import constants
 
+from ....harmonics import cache as c
+
 
 @nb.njit(cache=True)
-def gamma_ns(N, s1):
+def gamma_ns(N, cache):
     """
     Computes the leading-order non-singlet anomalous dimension.
 
@@ -17,15 +19,13 @@ def gamma_ns(N, s1):
     ----------
       N : complex
         Mellin moment
-      s1 : complex
-        harmonic sum :math:`S_{1}`
 
     Returns
     -------
       gamma_ns : complex
         Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)`
     """
-    gamma = -(3.0 - 4.0 * s1 + 2.0 / N / (N + 1.0))
+    gamma = -(3.0 - 4.0 * c.get(c.S1, cache, N) + 2.0 / N / (N + 1.0))
     result = constants.CF * gamma
     return result
 
@@ -77,7 +77,7 @@ def gamma_gq(N):
 
 
 @nb.njit(cache=True)
-def gamma_gg(N, s1, nf):
+def gamma_gg(N, nf, cache):
     """
     Computes the leading-order gluon-gluon anomalous dimension
 
@@ -87,8 +87,6 @@ def gamma_gg(N, s1, nf):
     ----------
       N : complex
         Mellin moment
-      s1 : complex
-        harmonic sum :math:`S_{1}`
       nf : int
         Number of active flavors
 
@@ -97,13 +95,13 @@ def gamma_gg(N, s1, nf):
       gamma_gg : complex
         Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)`
     """
-    gamma = s1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0)
+    gamma = c.get(c.S1, cache, N) - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0)
     result = constants.CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * constants.TR * nf
     return result
 
 
 @nb.njit(cache=True)
-def gamma_singlet(N, s1, nf):
+def gamma_singlet(N, nf, cache):
     r"""
       Computes the leading-order singlet anomalous dimension matrix
 
@@ -134,8 +132,9 @@ def gamma_singlet(N, s1, nf):
         gamma_gq : :math:`\gamma_{gq}^{(0)}`
         gamma_gg : :math:`\gamma_{gg}^{(0)}`
     """
-    gamma_qq = gamma_ns(N, s1)
+    gamma_qq = gamma_ns(N, cache)
     gamma_S_0 = np.array(
-        [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, s1, nf)]], np.complex_
+        [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, nf, cache)]],
+        np.complex_,
     )
     return gamma_S_0

src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py

@@ -5,8 +5,9 @@
 import numba as nb
 
 from eko import constants
+from eko.constants import zeta2, zeta3
+
 from .... import harmonics
-from ....harmonics.constants import zeta2, zeta3
 
 
 @nb.njit(cache=True)

src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py

@@ -10,8 +10,9 @@
 import numpy as np
 
 from eko import constants
+from eko.constants import log2, zeta2, zeta3
+
 from .... import harmonics
-from ....harmonics.constants import log2, zeta2, zeta3
 
 
 @nb.njit(cache=True)

src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py

@@ -8,7 +8,7 @@
 import numba as nb
 import numpy as np
 
-from ....harmonics.constants import zeta2, zeta3
+from eko.constants import zeta2, zeta3
 
 
 @nb.njit(cache=True)

src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py

@@ -3,8 +3,15 @@
 """
 import numba as nb
 
+<<<<<<< HEAD:src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py
+from eko.constants import CF, zeta3
+
+=======
 from eko.constants import CF
+
 from .....harmonics.constants import zeta3
+
+>>>>>>> master:src/eko/anomalous_dimensions/as4/gnsm.py
 from .....harmonics.log_functions import lm11m1, lm12m1, lm13m1