Python 3.13 & numpy 2.0 support

.github/workflows/unittests.yml

@@ -6,7 +6,7 @@ jobs:
   test:
     strategy:
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
       fail-fast: false
 
     uses: ./.github/workflows/python-poetry-tests-with-data.yml

CHANGELOG.md

@@ -18,6 +18,7 @@ Items without prefix refer to a global change.
 ### Changed
 - data: Change version naming for unreleased eko versions from `0.0.0` to `0.0.0-post.{distance}+{commit hash}` ([#448](https://github.com/NNPDF/eko/pull/448)) ([#465](https://github.com/NNPDF/eko/pull/465))
 - py: `dump_blocks` and `dump_info` in `ekobox.genpdf.export` now accept paths pointing to the actual output file. The actual output path is now also returned from these functions (soft breaking change) ([#469](https://github.com/NNPDF/eko/pull/469))
+- py: Add support for Python 3.13, lift minimum NumPy version to 2
 
 ### Fixed
 - py: Remove usage of `ev_op_iterations` from truncated evolution as this is an inconsistent choice ([#459](https://github.com/NNPDF/eko/pull/459))

benchmarks/ekobox/genpdf/benchmark_init.py

@@ -97,7 +97,7 @@ def benchmark_genpdf_dict(tmp_path, cd):
 @pytest.mark.isolated
 def benchmark_genpdf_custom(tmp_path, cd):
     with cd(tmp_path):
-        c = np.zeros_like(br.flavor_basis_pids, dtype=np.float_)
+        c = np.zeros_like(br.flavor_basis_pids, dtype=np.float64)
         c[br.flavor_basis_pids.index(1)] = 1.0
         c[br.flavor_basis_pids.index(2)] = 0.5
         genpdf.generate_pdf(

pyproject.toml

@@ -40,15 +40,16 @@ packages = [
 ]
 
 [tool.poetry.dependencies]
-python = "^3.9,<3.13"
-numpy = "^1.24"
+python = "^3.9,<3.14"
+numpy = "^2"
 scipy = "^1.10.1"
 PyYAML = "^6.0"
 lz4 = "^4.0.2"
-numba = "^0.59.0"
+# 0.60.0 is the last version compatible with 3.9 and 0.61 the first with 3.13
+numba = "~0.60.0 || ~0.61.0"
 packaging = "^24.1"
 # ekomark
-banana-hep = { version = "^0.6.12", optional = true }
+banana-hep = { version = "^0.6.16", optional = true }
 sqlalchemy = { version = "^1.4.21", optional = true }
 pandas = { version = "^2.1.4", optional = true }
 matplotlib = { version = "^3.5.1", optional = true }

src/eko/evolution_operator/operator_matrix_element.py

@@ -78,7 +78,7 @@ def build_ome(A, matching_order, a_s, backward_method):
     # Local bi = 1+x+x**2+x**3;
     # Print;
     # .end
-    ome = np.eye(len(A[0]), dtype=np.complex_)
+    ome = np.eye(len(A[0]), dtype=np.complex128)
     A = A[:, :, :]
     A = np.ascontiguousarray(A)
     if backward_method is MatchingMethods.BACKWARD_EXPANDED:

src/eko/evolution_operator/quad_ker.py

@@ -206,7 +206,7 @@ def build_ome(A, matching_order, a_s, backward_method):
     # Local bi = 1+x+x**2+x**3;
     # Print;
     # .end
-    ome = np.eye(len(A[0]), dtype=np.complex_)
+    ome = np.eye(len(A[0]), dtype=np.complex128)
     A = A[:, :, :]
     A = np.ascontiguousarray(A)
     if backward_method is InversionMethod.EXPANDED:

src/eko/interpolation.py

@@ -413,7 +413,7 @@ class XGrid:
     """
 
     def __init__(self, xgrid: Union[Sequence, npt.NDArray], log: bool = True):
-        ugrid = np.array(np.unique(xgrid), np.float_)
+        ugrid = np.array(np.unique(xgrid), np.float64)
         if len(xgrid) != len(ugrid):
             raise ValueError(f"xgrid is not unique: {xgrid}")
         if len(xgrid) < 2:

src/eko/kernels/non_singlet.py

@@ -27,7 +27,7 @@ def U_vec(gamma_ns, beta, order):
     np.ndarray
         U vector
     """
-    U = np.zeros(order[0], dtype=np.complex_)
+    U = np.zeros(order[0], dtype=np.complex128)
     beta0 = beta[0]
     R0 = gamma_ns[0] / beta0
     U[0] = 1.0

src/eko/kernels/non_singlet_qed.py

@@ -24,7 +24,7 @@ def contract_gammas(gamma_ns, aem):
         non-singlet anomalous dimensions
     """
     vec_alphaem = np.array(
-        [aem**i for i in range(gamma_ns.shape[1])], dtype=np.complex_
+        [aem**i for i in range(gamma_ns.shape[1])], dtype=np.complex128
     )
     return gamma_ns @ vec_alphaem
 

src/eko/kernels/singlet.py

@@ -324,12 +324,12 @@ def eko_iterate(gamma_singlet, a1, a0, beta_vec, order, ev_op_iterations):
         singlet iterated (exact) EKO
     """
     a_steps = np.geomspace(a0, a1, 1 + ev_op_iterations)
-    e = np.identity(2, np.complex_)
+    e = np.identity(2, np.complex128)
     al = a_steps[0]
     for ah in a_steps[1:]:
         a_half = (ah + al) / 2.0
         delta_a = ah - al
-        gamma_summed = np.zeros((2, 2), dtype=np.complex_)
+        gamma_summed = np.zeros((2, 2), dtype=np.complex128)
         beta_summed = 0
         for i in range(order[0]):
             gamma_summed += gamma_singlet[i] * a_half**i
@@ -368,7 +368,7 @@ def r_vec(gamma_singlet, beta, ev_op_max_order, order, is_exact):
         R vector
     """
     r = np.zeros(
-        (ev_op_max_order[0] + 1, 2, 2), dtype=np.complex_
+        (ev_op_max_order[0] + 1, 2, 2), dtype=np.complex128
     )  # k = 0 .. max_order
     beta0 = beta[0]
     # fill explicit elements
@@ -416,15 +416,15 @@ def u_vec(r, ev_op_max_order):
     numpy.ndarray
         U vector
     """
-    u = np.zeros((ev_op_max_order[0], 2, 2), np.complex_)  # k = 0 .. max_order
+    u = np.zeros((ev_op_max_order[0], 2, 2), np.complex128)  # k = 0 .. max_order
     # init
-    u[0] = np.identity(2, np.complex_)
+    u[0] = np.identity(2, np.complex128)
     _, r_p, r_m, e_p, e_m = ad.exp_matrix_2D(r[0])
     e_p = np.ascontiguousarray(e_p)
     e_m = np.ascontiguousarray(e_m)
     for kk in range(1, ev_op_max_order[0]):
         # compute R'
-        rp = np.zeros((2, 2), dtype=np.complex_)
+        rp = np.zeros((2, 2), dtype=np.complex128)
         for jj in range(kk):
             rp += np.ascontiguousarray(r[kk - jj]) @ u[jj]
         # now compose U
@@ -456,7 +456,7 @@ def sum_u(uvec, a):
         sum
     """
     p = 1.0
-    res = np.zeros((2, 2), dtype=np.complex_)
+    res = np.zeros((2, 2), dtype=np.complex128)
     for uk in uvec:
         res += p * uk
         p *= a
@@ -499,7 +499,7 @@ def eko_perturbative(
     """
     r = r_vec(gamma_singlet, beta, ev_op_max_order, order, is_exact)
     uk = u_vec(r, ev_op_max_order)
-    e = np.identity(2, np.complex_)
+    e = np.identity(2, np.complex128)
     # iterate elements
     a_steps = np.geomspace(a0, a1, 1 + ev_op_iterations)
     al = a_steps[0]
@@ -592,7 +592,7 @@ def dispatcher(  # pylint: disable=too-many-return-statements
     betalist = [beta.beta_qcd((2 + i, 0), nf) for i in range(order[0])]
     # for SV expanded we still fall in here, but we don't need to do anything
     if a1 == a0:
-        return np.eye(len(gamma_singlet[0]), dtype=np.complex_)
+        return np.eye(len(gamma_singlet[0]), dtype=np.complex128)
 
     # use always exact in LO
     if order[0] == 1:

src/eko/kernels/singlet_qed.py

@@ -35,7 +35,7 @@ def eko_iterate(gamma_singlet, as_list, a_half, nf, order, ev_op_iterations, dim
     e_s^{order} : numpy.ndarray
         singlet QEDxQCD iterated (exact) EKO
     """
-    e = np.identity(dim, np.complex_)
+    e = np.identity(dim, np.complex128)
     betaQCD = np.zeros((order[0] + 1, order[1] + 1))
     for i in range(1, order[0] + 1):
         betaQCD[i, 0] = beta.beta_qcd((i + 1, 0), nf)
@@ -46,7 +46,7 @@ def eko_iterate(gamma_singlet, as_list, a_half, nf, order, ev_op_iterations, dim
         as_half = a_half[step - 1, 0]
         aem = a_half[step - 1, 1]
         delta_a = ah - al
-        gamma = np.zeros((dim, dim), np.complex_)
+        gamma = np.zeros((dim, dim), np.complex128)
         betatot = 0
         for i in range(0, order[0] + 1):
             for j in range(0, order[1] + 1):

src/eko/scale_variations/expanded.py

@@ -153,7 +153,7 @@ def singlet_variation(gamma, a_s, order, nf, L, dim):
     numpy.ndarray
         scale variation kernel
     """
-    sv_ker = np.eye(dim, dtype=np.complex_)
+    sv_ker = np.eye(dim, dtype=np.complex128)
     gamma = np.ascontiguousarray(gamma)
     if order[0] >= 2:
         sv_ker += a_s * variation_as1(gamma, L)

src/ekobox/genpdf/__init__.py

@@ -150,7 +150,7 @@ def generate_pdf(
     >>> from eko import basis_rotation as br
     >>> from ekobox import genpdf
     >>> import numpy as np
-    >>> anti_qed_singlet = np.zeros_like(br.flavor_basis_pids, dtype=np.float_)
+    >>> anti_qed_singlet = np.zeros_like(br.flavor_basis_pids, dtype=np.float64)
     >>> anti_qed_singlet[br.flavor_basis_pids.index(1)] = -4
     >>> anti_qed_singlet[br.flavor_basis_pids.index(-1)] = -4
     >>> anti_qed_singlet[br.flavor_basis_pids.index(2)] = 1

src/ekobox/genpdf/load.py

@@ -66,8 +66,8 @@ def load_blocks_from_file(pdfset_name, member):
         # section head
         next_head_section = cnt.index("---\n", head_section + 1)
         # determine participating pids
-        xgrid = np.array(cnt[head_section + 1].strip().split(" "), dtype=np.float_)
-        Qgrid = np.array(cnt[head_section + 2].strip().split(" "), dtype=np.float_)
+        xgrid = np.array(cnt[head_section + 1].strip().split(" "), dtype=np.float64)
+        Qgrid = np.array(cnt[head_section + 2].strip().split(" "), dtype=np.float64)
         pids = np.array(cnt[head_section + 3].strip().split(" "), dtype=np.int_)
         # data
         data = []

src/ekore/anomalous_dimensions/__init__.py

@@ -53,7 +53,7 @@ def exp_matrix_2D(gamma_S):
     lambda_p = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] + det)
     lambda_m = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] - det)
     # compute projectors
-    identity = np.identity(2, np.complex_)
+    identity = np.identity(2, np.complex128)
     c = 1.0 / det
     e_p = +c * (gamma_S - lambda_m * identity)
     e_m = -c * (gamma_S - lambda_p * identity)
@@ -80,15 +80,15 @@ def exp_matrix(gamma):
         projectors on the eigenspaces of the matrix gamma :math:`\gamma(N)`
     """
     dim = gamma.shape[0]
-    e = np.zeros((dim, dim, dim), np.complex_)
+    e = np.zeros((dim, dim, dim), np.complex128)
     # if dim == 2:
     #     exp, lambda_p, lambda_m, e_p, e_m = exp_matrix_2D(gamma)
     #     e[0] = e_p
     #     e[1] = e_m
     #     return exp, np.array([lambda_p, lambda_m]), e
     w, v = np.linalg.eig(gamma)
     v_inv = np.linalg.inv(v)
-    exp = np.zeros((dim, dim), np.complex_)
+    exp = np.zeros((dim, dim), np.complex128)
     for i in range(dim):
         e[i] = np.outer(v[:, i], v_inv[i])
         exp += e[i] * np.exp(w[i])

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

@@ -36,7 +36,7 @@ def gamma_ns(order, mode, n, nf):
         non-singlet anomalous dimensions
     """
     cache = c.reset()
-    gamma_ns = np.zeros(order[0], np.complex_)
+    gamma_ns = np.zeros(order[0], np.complex128)
     gamma_ns[0] = as1.gamma_ns(n, cache)
     # NLO and beyond
     if order[0] >= 2:
@@ -81,7 +81,7 @@ def gamma_singlet(order, n, nf):
         singlet anomalous dimensions matrices
     """
     cache = c.reset()
-    gamma_s = np.zeros((order[0], 2, 2), np.complex_)
+    gamma_s = np.zeros((order[0], 2, 2), np.complex128)
     gamma_s[0] = as1.gamma_singlet(n, cache, nf)
     if order[0] >= 2:
         gamma_s[1] = as2.gamma_singlet(n, nf, cache)

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

@@ -106,6 +106,6 @@ def gamma_singlet(N, cache, nf):
     gamma_qq = gamma_ns(N, cache)
     gamma_S_0 = np.array(
         [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, cache, nf)]],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_0

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

@@ -193,6 +193,6 @@ def gamma_singlet(n, nf, cache):
             [gamma_qq, gamma_qg(n, nf, cache)],
             [gamma_gq(n, nf, cache), gamma_gg(n, nf, cache)],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_0

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

@@ -437,6 +437,6 @@ def gamma_singlet(N, nf, cache):
             [gamma_qq, gamma_qg(N, nf, cache)],
             [gamma_gq(N, nf, cache), gamma_gg(N, nf, cache)],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_0

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

@@ -50,7 +50,7 @@ def gamma_ns(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     """
     cache = c.reset()
     # now combine
-    gamma_ns = np.zeros(order[0], np.complex_)
+    gamma_ns = np.zeros(order[0], np.complex128)
     gamma_ns[0] = as1.gamma_ns(n, cache)
     # NLO and beyond
     if order[0] >= 2:
@@ -122,7 +122,7 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
         singlet anomalous dimensions matrices
     """
     cache = c.reset()
-    gamma_s = np.zeros((order[0], 2, 2), np.complex_)
+    gamma_s = np.zeros((order[0], 2, 2), np.complex128)
     gamma_s[0] = as1.gamma_singlet(n, cache, nf)
     if order[0] >= 2:
         gamma_s[1] = as2.gamma_singlet(n, nf, cache)
@@ -162,7 +162,7 @@ def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     """
     cache = c.reset()
     # now combine
-    gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex_)
+    gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex128)
     gamma_ns[1, 0] = as1.gamma_ns(n, cache)
     gamma_ns[0, 1] = choose_ns_ad_aem1(mode, n, cache)
     gamma_ns[1, 1] = choose_ns_ad_as1aem1(mode, n, cache)
@@ -310,7 +310,7 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
         singlet anomalous dimensions matrices
     """
     cache = c.reset()
-    gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex_)
+    gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex128)
     gamma_s[1, 0] = as1.gamma_singlet_qed(n, cache, nf)
     gamma_s[0, 1] = aem1.gamma_singlet(n, nf, cache)
     gamma_s[1, 1] = as1aem1.gamma_singlet(n, nf, cache)
@@ -354,7 +354,7 @@ def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
         valence anomalous dimensions matrices
     """
     cache = c.reset()
-    gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex_)
+    gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex128)
     gamma_v[1, 0] = as1.gamma_valence_qed(n, cache)
     gamma_v[0, 1] = aem1.gamma_valence(n, nf, cache)
     gamma_v[1, 1] = as1aem1.gamma_valence(n, nf, cache)

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

@@ -144,7 +144,7 @@ def gamma_singlet(N, nf, cache):
                 e2delta * gamma_nonsinglet,
             ],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_01
 
@@ -179,6 +179,6 @@ def gamma_valence(N, nf, cache):
             [e2avg, vue2m],
             [vde2m, e2delta],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_V_01 * gamma_ns(N, cache)

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

@@ -397,7 +397,7 @@ def gamma_singlet(N, nf, cache):
                 + vu * vd * e2m**2 * gamma_pure_singlet,
             ],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_02
 
@@ -437,6 +437,6 @@ def gamma_valence(N, nf, cache):
                 vd * constants.eu2 * gamma_ns_m_u + vu * constants.ed2 * gamma_ns_m_d,
             ],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_V_02

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

@@ -130,7 +130,7 @@ def gamma_singlet(N, cache, nf):
     gamma_qq = gamma_ns(N, cache)
     gamma_S_0 = np.array(
         [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, cache, nf)]],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_0
 
@@ -170,7 +170,7 @@ def gamma_singlet_qed(N, cache, nf):
             [gamma_qg(N, nf), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j],
             [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_qq],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S
 
@@ -203,6 +203,6 @@ def gamma_valence_qed(N, cache):
             [1.0, 0.0],
             [0.0, 1.0],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_V * gamma_ns(N, cache)

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

@@ -426,7 +426,7 @@ def gamma_singlet(N, nf, cache):
                 e2delta * gamma_ns_p,
             ],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_S_11
 
@@ -455,6 +455,6 @@ def gamma_valence(N, nf, cache):
             [e2avg, vue2m],
             [vde2m, e2delta],
         ],
-        np.complex_,
+        np.complex128,
     )
     return gamma_V_11 * gamma_nsm(N, cache)