From c6336024899ab34adbb5c96ba1ff0733b886b992 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 17 Feb 2022 14:58:49 +0100 Subject: [PATCH 01/31] Implement QED lo AD --- doc/source/refs.bib | 11 +++ src/eko/anomalous_dimensions/__init__.py | 12 +-- src/eko/anomalous_dimensions/aem1.py | 88 +++++++++++++++++++ .../anomalous_dimensions/{lo.py => as1.py} | 2 +- tests/test_ad.py | 2 +- tests/test_ad_lo.py | 2 +- 6 files changed, 108 insertions(+), 9 deletions(-) create mode 100644 src/eko/anomalous_dimensions/aem1.py rename src/eko/anomalous_dimensions/{lo.py => as1.py} (99%) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index 347b62de3..fde87c366 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -419,3 +419,14 @@ @article{Bertone:2013vaa pages = {1647--1668}, year = {2014} } + +@phdthesis{Carrazza:2015dea, + author = "Carrazza, Stefano", + title = "{Parton distribution functions with QED corrections}", + eprint = "1509.00209", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + doi = "10.13130/carrazza-stefano_phd2015-07-06", + school = "Milan U.", + year = "2015" +} diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 173fe5ead..b9a191944 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -20,7 +20,7 @@ import numba as nb import numpy as np -from . import harmonics, lo, nlo, nnlo +from . import as1, harmonics, nlo, nnlo @nb.njit("Tuple((c16[:,:],c16,c16,c16[:,:],c16[:,:]))(c16[:,:])", cache=True) @@ -52,7 +52,7 @@ def exp_singlet(gamma_S): See Also -------- - eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` """ @@ -94,7 +94,7 @@ def gamma_ns(order, mode, n, nf): See Also -------- - eko.anomalous_dimensions.lo.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)` + eko.anomalous_dimensions.as1.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)` eko.anomalous_dimensions.nlo.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)` eko.anomalous_dimensions.nlo.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)` eko.anomalous_dimensions.nnlo.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)` @@ -105,7 +105,7 @@ def gamma_ns(order, mode, n, nf): sx = np.full(1, harmonics.harmonic_S1(n)) # now combine gamma_ns = np.zeros(order + 1, np.complex_) - gamma_ns[0] = lo.gamma_ns_0(n, sx[0]) + gamma_ns[0] = as1.gamma_ns_0(n, sx[0]) # NLO and beyond if order >= 1: # TODO: pass the necessary harmonics to nlo gammas @@ -150,7 +150,7 @@ def gamma_singlet(order, n, nf): See Also -------- - eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` """ @@ -161,7 +161,7 @@ def gamma_singlet(order, n, nf): sx = np.append(sx, harmonics.harmonic_S3(n)) gamma_singlet = np.zeros((order + 1, 2, 2), np.complex_) - gamma_singlet[0] = lo.gamma_singlet_0(n, sx[0], nf) + gamma_singlet[0] = as1.gamma_singlet_0(n, sx[0], nf) if order >= 1: gamma_singlet[1] = nlo.gamma_singlet_1(n, nf) if order == 2: diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py new file mode 100644 index 000000000..dd5d26636 --- /dev/null +++ b/src/eko/anomalous_dimensions/aem1.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- +import numba as nb + +from .. import constants +from . import as1 + + +@nb.njit("c16(c16)", cache=True) +def gamma_phq_0(N): + """ + Computes the leading-order photon-quark anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Parameters + ---------- + N : complex + Mellin moment + + Returns + ------- + gamma_phq_0 : complex + Leading-order photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(0)}(N)` + """ + + return as1.gamma_gq_0(N) / constants.CF + + +@nb.njit("c16(c16,u1)", cache=True) +def gamma_qph_0(N, nf: int): + """ + Computes the leading-order quark-photon anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + But adding the :math:`N_C` and the :math:`2n_f` factors from :math:`\\theta` inside the + definition of :math:`\\gamma_{q \\gamma}^{(0)}(N)`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + + Returns + ------- + gamma_qph_0 : complex + Leading-order quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(0)}(N)` + """ + return as1.gamma_qg_0(N, nf) / constants.TR * constants.NC + + +@nb.njit("c16()", cache=True) +def gamma_phph_0(): + """ + Computes the leading-order photon-photon anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Returns + ------- + gamma_phph_0 : complex + Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)` + """ + + return -4.0 / 3 + + +@nb.njit("c16(c16,c16)", cache=True) +def gamma_ns_0(N, s1): + """ + Computes the leading-order non-singlet QED anomalous dimension. + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Parameters + ---------- + N : complex + Mellin moment + s1 : complex + S1(N) + + Returns + ------- + gamma_ns_0 : complex + Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` + """ + return as1.gamma_ns_0(N, s1) / constants.CF diff --git a/src/eko/anomalous_dimensions/lo.py b/src/eko/anomalous_dimensions/as1.py similarity index 99% rename from src/eko/anomalous_dimensions/lo.py rename to src/eko/anomalous_dimensions/as1.py index 0c5e23c50..952ac4839 100644 --- a/src/eko/anomalous_dimensions/lo.py +++ b/src/eko/anomalous_dimensions/as1.py @@ -4,7 +4,7 @@ import numba as nb import numpy as np -from eko import constants +from .. import constants @nb.njit("c16(c16,c16)", cache=True) diff --git a/tests/test_ad.py b/tests/test_ad.py index 3363220c1..f60dcc4c7 100644 --- a/tests/test_ad.py +++ b/tests/test_ad.py @@ -4,8 +4,8 @@ from numpy.testing import assert_allclose, assert_almost_equal from eko import anomalous_dimensions as ad +from eko.anomalous_dimensions import as1 as ad_lo from eko.anomalous_dimensions import harmonics -from eko.anomalous_dimensions import lo as ad_lo NF = 5 diff --git a/tests/test_ad_lo.py b/tests/test_ad_lo.py index 950116295..0dc9cb078 100644 --- a/tests/test_ad_lo.py +++ b/tests/test_ad_lo.py @@ -2,7 +2,7 @@ # Test LO splitting functions import numpy as np -import eko.anomalous_dimensions.lo as ad_lo +import eko.anomalous_dimensions.as1 as ad_lo from eko.anomalous_dimensions import harmonics NF = 5 From 38e11464606742f8266d71ce42452fd84325675d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 17 Feb 2022 18:07:01 +0100 Subject: [PATCH 02/31] Change labels --- src/eko/basis_rotation.py | 4 ++-- src/eko/constants.py | 15 ++++++--------- src/eko/evolution_operator/__init__.py | 10 +++++----- 3 files changed, 13 insertions(+), 16 deletions(-) diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index e2d46fc09..887e03ccd 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -69,8 +69,8 @@ ) """|pid| representation of :data:`evol_basis`.""" -singlet_labels = ("S_qq", "S_qg", "S_gq", "S_gg") -non_singlet_labels = ("NS_m", "NS_p", "NS_v") +singlet_labels = (("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")) +non_singlet_labels = (("V3", "V3"), ("T3", "T3"), ("V", "V")) full_labels = (*singlet_labels, *non_singlet_labels) anomalous_dimensions_basis = full_labels r""" diff --git a/src/eko/constants.py b/src/eko/constants.py index 06f96ce93..3de4379b3 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -2,15 +2,6 @@ r""" This files sets the physical constants. - -The constants are: - - - :math:`N_C` the number of colors - defaults to :math:`3` - - :math:`T_R` the normalization of fundamental generators - defaults to :math:`1/2` - - :math:`C_A` second Casimir constant in the adjoint representation - defaults to - :math:`N_C = 3` - - :math:`C_F` second Casimir constant in the fundamental representation - defaults to - :math:`\frac{N_C^2-1}{2N_C} = 4/3` """ NC = 3 @@ -24,3 +15,9 @@ CF = float((NC * NC - 1.0) / (2.0 * NC)) """second Casimir constant in the fundamental representation""" + +eu2 = 4.0 / 9 +"""up quarks charge squared""" + +ed2 = 1.0 / 9 +"""down quarks charge squared""" diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index f42b972b6..e56569992 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -242,16 +242,16 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append("NS_p") + labels.append(("T3", "T3")) if order >= 1: # - becomes different starting from NLO - labels.append("NS_m") + labels.append(("V3", "V3")) if order >= 2: # v also becomes different starting from NNLO - labels.append("NS_v") + labels.append(("V", "V")) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") else: - labels.extend(singlet_labels) + labels.extend([("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")]) return labels def compute(self): @@ -267,7 +267,7 @@ def compute(self): for n in full_labels: if n in labels: # off diag singlet are zero - if n in ["S_qg", "S_gq"]: + if n in [("S", "g"), ("g", "S")]: self.op_members[n] = zero.copy() else: self.op_members[n] = eye.copy() From d0bd345f093e11529e433f975bd341e62aa748e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 18 Feb 2022 17:13:21 +0100 Subject: [PATCH 03/31] Start changing ad names --- src/eko/basis_rotation.py | 2 +- src/eko/evolution_operator/__init__.py | 22 ++++++++++++---------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index 826acf3fb..658eeee05 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -70,7 +70,7 @@ """|pid| representation of :data:`evol_basis`.""" singlet_labels = (("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")) -non_singlet_labels = (("V3", "V3"), ("T3", "T3"), ("V", "V")) +non_singlet_labels = (("ns-", None), ("ns+", None), ("nsV", None)) full_labels = (*singlet_labels, *non_singlet_labels) anomalous_dimensions_basis = full_labels r""" diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index e56569992..69470225c 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -97,11 +97,12 @@ def gamma_singlet_fact(order, n, nf, L): return gamma_singlet -@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) +@nb.njit("f8(f8,u1,string,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) def quad_ker( u, order, - mode, + mode0, + mode1, method, is_log, logx, @@ -150,7 +151,7 @@ def quad_ker( ker : float evaluated integration kernel """ - is_singlet = mode[0] == "S" + is_singlet = mode1 is not None # get transformation to N integral if logx == 0.0: return 0.0 @@ -175,11 +176,11 @@ def quad_ker( order, method, gamma_singlet, a1, a0, nf, ev_op_iterations, ev_op_max_order ) # select element of matrix - k = 0 if mode[2] == "q" else 1 - l = 0 if mode[3] == "q" else 1 + k = 0 if mode0 == "S" else 1 + l = 0 if mode1 == "S" else 1 ker = ker[k, l] else: - gamma_ns = gamma_ns_fact(order, mode, n, nf, L) + gamma_ns = gamma_ns_fact(order, mode0, n, nf, L) ker = ns.dispatcher( order, method, @@ -242,11 +243,11 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append(("T3", "T3")) + labels.append(("ns+", None)) if order >= 1: # - becomes different starting from NLO - labels.append(("V3", "V3")) + labels.append(("ns-", None)) if order >= 2: # v also becomes different starting from NNLO - labels.append(("V", "V")) + labels.append(("nsV", None)) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") @@ -318,7 +319,8 @@ def compute(self): 1.0 - self._mellin_cut, args=( self.config["order"], - label, + label[0], + label[1], self.config["method"], int_disp.log, logx, From 2c13df5e8f3d1bc3e3bde9e374b6ace14887544e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 18 Feb 2022 17:53:23 +0100 Subject: [PATCH 04/31] Start fixing tests --- src/eko/anomalous_dimensions/__init__.py | 12 +++--- src/eko/basis_rotation.py | 14 +++---- src/eko/evolution_operator/__init__.py | 18 ++++++--- src/eko/evolution_operator/physical.py | 16 ++++---- tests/test_ad.py | 8 ++-- tests/test_basis_rotation.py | 6 +-- tests/test_ev_operator.py | 48 +++++++++++++----------- 7 files changed, 68 insertions(+), 54 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index b9a191944..b708e3c68 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -80,7 +80,7 @@ def gamma_ns(order, mode, n, nf): ---------- order : int perturbative order - mode : "m" | "p" | "v" + mode : "-" | "+" | "V" sector identifier n : complex Mellin variable @@ -109,21 +109,21 @@ def gamma_ns(order, mode, n, nf): # NLO and beyond if order >= 1: # TODO: pass the necessary harmonics to nlo gammas - if mode == "p": + if mode == "+": gamma_ns_1 = nlo.gamma_nsp_1(n, nf) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here - elif mode in ["m", "v"]: + elif mode in ["-", "V"]: gamma_ns_1 = nlo.gamma_nsm_1(n, nf) gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order >= 2: sx = np.append(sx, harmonics.harmonic_S2(n)) sx = np.append(sx, harmonics.harmonic_S3(n)) - if mode == "p": + if mode == "+": gamma_ns_2 = -nnlo.gamma_nsp_2(n, nf, sx) - elif mode == "m": + elif mode == "-": gamma_ns_2 = -nnlo.gamma_nsm_2(n, nf, sx) - elif mode == "v": + elif mode == "V": gamma_ns_2 = -nnlo.gamma_nsv_2(n, nf, sx) gamma_ns[2] = gamma_ns_2 return gamma_ns diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index 658eeee05..e1e149d4b 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -101,13 +101,13 @@ """ map_ad_to_evolution = { - "S_qq": ["S.S"], - "S_qg": ["S.g"], - "S_gq": ["g.S"], - "S_gg": ["g.g"], - "NS_v": ["V.V"], - "NS_p": ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], - "NS_m": ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], + ("S", "S"): ["S.S"], + ("S", "g"): ["S.g"], + ("g", "S"): ["g.S"], + ("g", "g"): ["g.g"], + ("nsV", None): ["V.V"], + ("ns+", None): ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], + ("ns-", None): ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], } """ Map anomalous dimension sectors' names to their members diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 69470225c..cb3ffea7a 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -357,9 +357,17 @@ def copy_ns_ops(self): """Copy non-singlet kernels, if necessary""" order = self.config["order"] if order == 0: # in LO +=-=v - for label in ["NS_v", "NS_m"]: - self.op_members[label].value = self.op_members["NS_p"].value.copy() - self.op_members[label].error = self.op_members["NS_p"].error.copy() + for label in ["nsV", "ns-"]: + self.op_members[(label, None)].value = self.op_members[ + ("ns+", None) + ].value.copy() + self.op_members[(label, None)].error = self.op_members[ + ("ns+", None) + ].error.copy() elif order == 1: # in NLO -=v - self.op_members["NS_v"].value = self.op_members["NS_m"].value.copy() - self.op_members["NS_v"].error = self.op_members["NS_m"].error.copy() + self.op_members[("nsV", None)].value = self.op_members[ + ("ns-", None) + ].value.copy() + self.op_members[("nsV", None)].error = self.op_members[ + ("ns-", None) + ].error.copy() diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index 91548eac7..f7e539dc8 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -48,21 +48,21 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range): """ # constant elements m = { - "S.S": op_members["S_qq"], - "S.g": op_members["S_qg"], - "g.g": op_members["S_gg"], - "g.S": op_members["S_gq"], - "V.V": op_members["NS_v"], + "S.S": op_members[("S", "S")], + "S.g": op_members[("S", "g")], + "g.g": op_members[("g", "g")], + "g.S": op_members[("g", "S")], + "V.V": op_members[("nsV", None)], } # add elements which are already active for f in range(2, nf + 1): n = f**2 - 1 - m[f"V{n}.V{n}"] = op_members["NS_m"] - m[f"T{n}.T{n}"] = op_members["NS_p"] + m[f"V{n}.V{n}"] = op_members[("ns-", None)] + m[f"T{n}.T{n}"] = op_members[("ns+", None)] # deal with intrinsic heavy quark PDFs if intrinsic_range is not None: hqfl = "cbt" - op_id = member.OpMember.id_like(op_members["NS_v"]) + op_id = member.OpMember.id_like(op_members[("nsV", None)]) for intr_fl in intrinsic_range: if intr_fl <= nf: # light quarks are not intrinsic continue diff --git a/tests/test_ad.py b/tests/test_ad.py index f60dcc4c7..decbb2943 100644 --- a/tests/test_ad.py +++ b/tests/test_ad.py @@ -53,9 +53,9 @@ def test_eigensystem_gamma_singlet_projectors_EV(): def test_gamma_ns(): nf = 3 # LO - assert_almost_equal(ad.gamma_ns(2, "p", 1, nf)[0], 0.0) + assert_almost_equal(ad.gamma_ns(2, "+", 1, nf)[0], 0.0) # NLO - assert_allclose(ad.gamma_ns(1, "m", 1, nf), np.zeros(2), atol=2e-6) + assert_allclose(ad.gamma_ns(1, "-", 1, nf), np.zeros(2), atol=2e-6) # NNLO - assert_allclose(ad.gamma_ns(2, "m", 1, nf), np.zeros(3), atol=2e-4) - assert_allclose(ad.gamma_ns(2, "v", 1, nf), np.zeros(3), atol=8e-4) + assert_allclose(ad.gamma_ns(2, "-", 1, nf), np.zeros(3), atol=2e-4) + assert_allclose(ad.gamma_ns(2, "V", 1, nf), np.zeros(3), atol=8e-4) diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py index c55a49f78..9a04fa8b7 100644 --- a/tests/test_basis_rotation.py +++ b/tests/test_basis_rotation.py @@ -10,19 +10,19 @@ def test_ad_projector(): g = br.rotate_flavor_to_evolution[2] v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")] - s_to_s = br.ad_projector("S_qq", nf=6) + s_to_s = br.ad_projector(("S", "S"), nf=6) np.testing.assert_allclose(s @ s_to_s, s) np.testing.assert_allclose(g @ s_to_s, 0.0) np.testing.assert_allclose(v3 @ s_to_s, 0.0) - g_to_s = br.ad_projector("S_gq", nf=6) + g_to_s = br.ad_projector(("g", "S"), nf=6) np.testing.assert_allclose(s @ g_to_s, 0.0) np.testing.assert_allclose(g @ g_to_s, s) np.testing.assert_allclose(v3 @ g_to_s, 0.0) - ns_m = br.ad_projector("NS_m", nf=6) + ns_m = br.ad_projector(("ns-", None), nf=6) np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15) np.testing.assert_allclose(g @ ns_m, 0.0) diff --git a/tests/test_ev_operator.py b/tests/test_ev_operator.py index 4a657e575..4c4d2e8f6 100644 --- a/tests/test_ev_operator.py +++ b/tests/test_ev_operator.py @@ -18,13 +18,13 @@ def test_gamma_ns_fact(monkeypatch): gamma_ns = np.array([1.0, 0.5, 0.25]) monkeypatch.setattr(ad, "gamma_ns", lambda *args: gamma_ns.copy()) - gamma_ns_LO_0 = gamma_ns_fact(0, "NS_p", 1, 3, 0) + gamma_ns_LO_0 = gamma_ns_fact(0, "ns+", 1, 3, 0) np.testing.assert_allclose(gamma_ns_LO_0, gamma_ns) - gamma_ns_LO_1 = gamma_ns_fact(0, "NS_p", 1, 3, 1) + gamma_ns_LO_1 = gamma_ns_fact(0, "ns+", 1, 3, 1) np.testing.assert_allclose(gamma_ns_LO_1, gamma_ns) - gamma_ns_NLO_1 = gamma_ns_fact(1, "NS_p", 1, 3, 1) + gamma_ns_NLO_1 = gamma_ns_fact(1, "ns+", 1, 3, 1) assert gamma_ns_NLO_1[1] < gamma_ns[1] - gamma_ns_NNLO_1 = gamma_ns_fact(2, "NS_p", 1, 3, 1) + gamma_ns_NNLO_1 = gamma_ns_fact(2, "ns+", 1, 3, 1) assert gamma_ns_NNLO_1[2] - gamma_ns[2] == 8.0 @@ -56,7 +56,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode="NS_p", + mode0="ns+", + mode1=None, method="", is_log=is_log, logx=0.0, @@ -72,7 +73,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode="S_qq", + mode0="S", + mode1="S", method="", is_log=is_log, logx=1.0, @@ -88,7 +90,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode="S_qg", + mode0="S", + mode1="g", method="", is_log=is_log, logx=0.0, @@ -105,7 +108,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode="NS_p", + mode0="ns+", + mode1=None, method="", is_log=True, logx=0.0, @@ -129,9 +133,7 @@ def test_labels(self): 1, 2, ) - assert sorted(o.labels()) == sorted( - ["NS_p", "NS_m", "NS_v", "S_qq", "S_qg", "S_gq", "S_gg"] - ) + assert sorted(o.labels()) == sorted(br.full_labels) o = Operator( dict(order=1, debug_skip_non_singlet=True, debug_skip_singlet=True), {}, @@ -191,19 +193,21 @@ def test_compute(self, monkeypatch): ) # LO o.compute() - assert "NS_m" in o.op_members + assert ("ns-", None) in o.op_members np.testing.assert_allclose( - o.op_members["NS_m"].value, o.op_members["NS_p"].value + o.op_members[("ns-", None)].value, o.op_members[("ns+", None)].value ) np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_p"].value + o.op_members[("nsV", None)].value, o.op_members[("ns+", None)].value ) # NLO o.config["order"] = 1 o.compute() - assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) + assert not np.allclose( + o.op_members[("ns+", None)].value, o.op_members[("ns-", None)].value + ) np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_m"].value + o.op_members[("nsV", None)].value, o.op_members[("ns-", None)].value ) # unity operators @@ -218,13 +222,13 @@ def test_compute(self, monkeypatch): def test_pegasus_path(): def quad_ker_pegasus( - u, order, mode, method, logx, areas, a1, a0, nf, L, ev_op_iterations + u, order, mode0, method, logx, areas, a1, a0, nf, L, ev_op_iterations ): # compute the mellin inversion as done in pegasus phi = 3 / 4 * np.pi c = 1.9 n = complex(c + u * np.exp(1j * phi)) - gamma_ns = gamma_ns_fact(order, mode, n, nf, L) + gamma_ns = gamma_ns_fact(order, mode0, n, nf, L) ker = ns.dispatcher( order, method, @@ -242,7 +246,8 @@ def quad_ker_pegasus( xgrid = np.geomspace(1e-7, 1, 10) int_disp = InterpolatorDispatcher(xgrid, 1, True) order = 1 - mode = "NS_p" + mode0 = "ns+" + mode1 = None method = "" logxs = np.log(int_disp.xgrid_raw) a1 = 1 @@ -258,7 +263,8 @@ def quad_ker_pegasus( 1.0, args=( order, - mode, + mode0, + mode1, method, int_disp.log, logx, @@ -282,7 +288,7 @@ def quad_ker_pegasus( np.inf, args=( order, - mode, + mode0, method, logx, bf.areas_representation, From 7ce13e2718e5062532a2b58e9cb1e72983d2e412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 22 Feb 2022 16:36:21 +0100 Subject: [PATCH 05/31] Change non-singlet marker to empty string --- benchmarks/lha_paper_bench.py | 10 +++++++--- src/eko/basis_rotation.py | 8 ++++---- src/eko/evolution_operator/__init__.py | 24 ++++++++++++------------ src/eko/evolution_operator/physical.py | 8 ++++---- tests/test_basis_rotation.py | 2 +- tests/test_ev_operator.py | 10 +++++----- 6 files changed, 33 insertions(+), 29 deletions(-) diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py index 9ef572229..3021bf0bb 100644 --- a/benchmarks/lha_paper_bench.py +++ b/benchmarks/lha_paper_bench.py @@ -3,9 +3,13 @@ Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO) """ import numpy as np +from banana import register from ekomark.benchmark.runner import Runner +register(__file__) + + base_theory = { "ModEv": "EXA", "Q0": np.sqrt( @@ -170,8 +174,8 @@ def skip_pdfs(theory): if __name__ == "__main__": - obj = BenchmarkVFNS() - # obj = BenchmarkFFNS() + # obj = BenchmarkVFNS() + obj = BenchmarkFFNS() - obj.benchmark_plain(2) + obj.benchmark_plain(0) # obj.benchmark_sv(1) diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index e1e149d4b..fc63f36e1 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -70,7 +70,7 @@ """|pid| representation of :data:`evol_basis`.""" singlet_labels = (("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")) -non_singlet_labels = (("ns-", None), ("ns+", None), ("nsV", None)) +non_singlet_labels = (("ns-", ""), ("ns+", ""), ("nsV", "")) full_labels = (*singlet_labels, *non_singlet_labels) anomalous_dimensions_basis = full_labels r""" @@ -105,9 +105,9 @@ ("S", "g"): ["S.g"], ("g", "S"): ["g.S"], ("g", "g"): ["g.g"], - ("nsV", None): ["V.V"], - ("ns+", None): ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], - ("ns-", None): ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], + ("nsV", ""): ["V.V"], + ("ns+", ""): ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], + ("ns-", ""): ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], } """ Map anomalous dimension sectors' names to their members diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index cb3ffea7a..c12cb073b 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -151,7 +151,7 @@ def quad_ker( ker : float evaluated integration kernel """ - is_singlet = mode1 is not None + is_singlet = len(mode1) > 0 # get transformation to N integral if logx == 0.0: return 0.0 @@ -243,11 +243,11 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append(("ns+", None)) + labels.append(("ns+", "")) if order >= 1: # - becomes different starting from NLO - labels.append(("ns-", None)) + labels.append(("ns-", "")) if order >= 2: # v also becomes different starting from NNLO - labels.append(("nsV", None)) + labels.append(("nsV", "")) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") @@ -358,16 +358,16 @@ def copy_ns_ops(self): order = self.config["order"] if order == 0: # in LO +=-=v for label in ["nsV", "ns-"]: - self.op_members[(label, None)].value = self.op_members[ - ("ns+", None) + self.op_members[(label, "")].value = self.op_members[ + ("ns+", "") ].value.copy() - self.op_members[(label, None)].error = self.op_members[ - ("ns+", None) + self.op_members[(label, "")].error = self.op_members[ + ("ns+", "") ].error.copy() elif order == 1: # in NLO -=v - self.op_members[("nsV", None)].value = self.op_members[ - ("ns-", None) + self.op_members[("nsV", "")].value = self.op_members[ + ("ns-", "") ].value.copy() - self.op_members[("nsV", None)].error = self.op_members[ - ("ns-", None) + self.op_members[("nsV", "")].error = self.op_members[ + ("ns-", "") ].error.copy() diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index f7e539dc8..a8bc05709 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -52,17 +52,17 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range): "S.g": op_members[("S", "g")], "g.g": op_members[("g", "g")], "g.S": op_members[("g", "S")], - "V.V": op_members[("nsV", None)], + "V.V": op_members[("nsV", "")], } # add elements which are already active for f in range(2, nf + 1): n = f**2 - 1 - m[f"V{n}.V{n}"] = op_members[("ns-", None)] - m[f"T{n}.T{n}"] = op_members[("ns+", None)] + m[f"V{n}.V{n}"] = op_members[("ns-", "")] + m[f"T{n}.T{n}"] = op_members[("ns+", "")] # deal with intrinsic heavy quark PDFs if intrinsic_range is not None: hqfl = "cbt" - op_id = member.OpMember.id_like(op_members[("nsV", None)]) + op_id = member.OpMember.id_like(op_members[("nsV", "")]) for intr_fl in intrinsic_range: if intr_fl <= nf: # light quarks are not intrinsic continue diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py index 9a04fa8b7..c102af936 100644 --- a/tests/test_basis_rotation.py +++ b/tests/test_basis_rotation.py @@ -22,7 +22,7 @@ def test_ad_projector(): np.testing.assert_allclose(g @ g_to_s, s) np.testing.assert_allclose(v3 @ g_to_s, 0.0) - ns_m = br.ad_projector(("ns-", None), nf=6) + ns_m = br.ad_projector(("ns-", ""), nf=6) np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15) np.testing.assert_allclose(g @ ns_m, 0.0) diff --git a/tests/test_ev_operator.py b/tests/test_ev_operator.py index 4c4d2e8f6..c34d97b16 100644 --- a/tests/test_ev_operator.py +++ b/tests/test_ev_operator.py @@ -193,21 +193,21 @@ def test_compute(self, monkeypatch): ) # LO o.compute() - assert ("ns-", None) in o.op_members + assert ("ns-", "") in o.op_members np.testing.assert_allclose( - o.op_members[("ns-", None)].value, o.op_members[("ns+", None)].value + o.op_members[("ns-", "")].value, o.op_members[("ns+", "")].value ) np.testing.assert_allclose( - o.op_members[("nsV", None)].value, o.op_members[("ns+", None)].value + o.op_members[("nsV", "")].value, o.op_members[("ns+", "")].value ) # NLO o.config["order"] = 1 o.compute() assert not np.allclose( - o.op_members[("ns+", None)].value, o.op_members[("ns-", None)].value + o.op_members[("ns+", "")].value, o.op_members[("ns-", "")].value ) np.testing.assert_allclose( - o.op_members[("nsV", None)].value, o.op_members[("ns-", None)].value + o.op_members[("nsV", "")].value, o.op_members[("ns-", "")].value ) # unity operators From f2c703b283791f35307abf479ddd71c3c42894e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 22 Feb 2022 17:43:07 +0100 Subject: [PATCH 06/31] Fic matching conditions --- src/eko/matching_conditions/__init__.py | 28 +++--- .../operator_matrix_element.py | 43 ++++++--- tests/test_ev_operator.py | 6 +- tests/test_matching.py | 27 +++--- tests/test_ome.py | 90 +++++++++++-------- 5 files changed, 114 insertions(+), 80 deletions(-) diff --git a/src/eko/matching_conditions/__init__.py b/src/eko/matching_conditions/__init__.py index 3b7ae3c72..ce15e7a75 100644 --- a/src/eko/matching_conditions/__init__.py +++ b/src/eko/matching_conditions/__init__.py @@ -44,11 +44,11 @@ def split_ad_to_evol_map( """ m = { - "S.S": ome_members["S_qq"], - "S.g": ome_members["S_qg"], # This is always zero for the time being - "g.S": ome_members["S_gq"], - "g.g": ome_members["S_gg"], - "V.V": ome_members["NS_qq"], + "S.S": ome_members[("S", "S")], + "S.g": ome_members[("S", "g")], # This is always zero for the time being + "g.S": ome_members[("g", "S")], + "g.g": ome_members[("g", "g")], + "V.V": ome_members[("V", "V")], } # add elements which are already active @@ -61,15 +61,15 @@ def split_ad_to_evol_map( hq = br.quark_names[nf] m.update( { - f"{hq}-.V": ome_members["NS_Hq"], - f"{hq}+.S": ome_members["S_Hq"], - f"{hq}+.g": ome_members["S_Hg"], + f"{hq}-.V": ome_members[("H-", "V")], + f"{hq}+.S": ome_members[("H+", "S")], + f"{hq}+.g": ome_members[("H+", "g")], } ) # intrinsic matching if len(intrinsic_range) != 0: - op_id = member.OpMember.id_like(ome_members["NS_qq"]) + op_id = member.OpMember.id_like(ome_members[("V", "V")]) for intr_fl in intrinsic_range: ihq = br.quark_names[intr_fl - 1] # find name if intr_fl > nf + 1: @@ -80,11 +80,11 @@ def split_ad_to_evol_map( # match the missing contribution from h+ and h- m.update( { - f"{ihq}+.{ihq}+": ome_members["S_HH"], - # f"S.{ihq}+": ome_members["S_qH"], - f"g.{ihq}+": ome_members["S_gH"], - f"{ihq}-.{ihq}-": ome_members["NS_HH"], - # f"V.{ihq}-": ome_members["NS_qH"], + f"{ihq}+.{ihq}+": ome_members[("H+", "H+")], + # f"S.{ihq}+": ome_members[("S", "H+")], + f"g.{ihq}+": ome_members[("g", "H+")], + f"{ihq}-.{ihq}-": ome_members[("H-", "H-")], + # f"V.{ihq}-": ome_members[("V", "H-")], } ) return cls.promote_names(m, q2_thr) diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index 708ca218e..aa068bfdf 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -145,8 +145,10 @@ def build_ome(A, order, a_s, backward_method): return ome -@nb.njit("f8(f8,u1,string,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) -def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_msbar): +@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) +def quad_ker( + u, order, mode0, mode1, is_log, logx, areas, a_s, L, backward_method, is_msbar +): """ Raw kernel inside quad @@ -177,15 +179,15 @@ def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_ms ker : float evaluated integration kernel """ - is_singlet = mode[0] == "S" + is_singlet = mode0 in ["S", "g", "H+"] # get transformation to N integral r = 0.4 * 16.0 / (1.0 - logx) if is_singlet: o = 1.0 - indeces = {"g": 0, "q": 1, "H": 2} + indices = {"g": 0, "S": 1, "H+": 2} else: o = 0.0 - indeces = {"q": 0, "H": 1} + indices = {"V": 0, "H-": 1} n = mellin.Talbot_path(u, r, o) jac = mellin.Talbot_jac(u, r, o) @@ -215,7 +217,7 @@ def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_ms ker = build_ome(A, order, a_s, backward_method) # select the need matrix element - ker = ker[indeces[mode[-2]], indeces[mode[-1]]] + ker = ker[indices[mode0], indices[mode1]] if ker == 0.0: return 0.0 @@ -267,10 +269,10 @@ def labels(self): if self.config["debug_skip_non_singlet"]: logger.warning("Matching: skipping non-singlet sector") else: - labels.extend(["NS_qq", "NS_Hq"]) + labels.extend([("V", "V"), ("H-", "V")]) if self.is_intrinsic or self.backward_method != "": # intrisic labels, which are not zero at NLO - labels.append("NS_HH") + labels.append(("H-", "H-")) # if self.backward_method == "exact": # # this contribution starts at NNLO, we don't have it for the moment # labels.append("NS_qH") @@ -279,9 +281,9 @@ def labels(self): if self.config["debug_skip_singlet"]: logger.warning("Matching: skipping singlet sector") else: - labels.extend([*singlet_labels, "S_Hg", "S_Hq"]) + labels.extend([*singlet_labels, ("H+", "g"), ("H+", "S")]) if self.is_intrinsic or self.backward_method != "": - labels.extend(["S_gH", "S_HH"]) + labels.extend([("g", "H+"), ("H+", "H+")]) # if self.backward_method == "exact": # labels.extend(["S_qH"]) return labels @@ -304,7 +306,7 @@ def compute(self, q2, L, is_msbar): grid_size = len(self.int_disp.xgrid) labels = self.labels() for n in labels: - if n[-1] == n[-2]: + if n[0] == n[1]: self.ome_members[n] = OpMember( np.eye(grid_size), np.zeros((grid_size, grid_size)) ) @@ -339,7 +341,8 @@ def compute(self, q2, L, is_msbar): 1.0 - self._mellin_cut, args=( self.config["order"], - label, + label[0], + label[1], self.int_disp.log, logx, bf.areas_representation, @@ -372,7 +375,13 @@ def copy_ome(self): """Add the missing |OME|, if necessary""" grid_size = len(self.int_disp.xgrid) # basic labels skipped with skip debug - for label in ["NS_qq", "S_Hg", "S_Hq", "NS_Hq", *singlet_labels]: + for label in [ + ("V", "V"), + ("H+", "g"), + ("H+", "S"), + ("H-", "V"), + *singlet_labels, + ]: if label not in self.ome_members: self.ome_members[label] = OpMember( np.zeros((grid_size, grid_size)), np.zeros((grid_size, grid_size)) @@ -380,7 +389,13 @@ def copy_ome(self): # intrinsic labels not computed yet if self.is_intrinsic: - for label in ["S_qH", "NS_qH", "NS_HH", "S_HH", "S_gH"]: + for label in [ + ("S", "H+"), + ("V", "H-"), + ("H-", "H-"), + ("H+", "H+"), + ("g", "H+"), + ]: if label not in self.ome_members: self.ome_members[label] = OpMember( np.zeros((grid_size, grid_size)), diff --git a/tests/test_ev_operator.py b/tests/test_ev_operator.py index c34d97b16..6e5df0b18 100644 --- a/tests/test_ev_operator.py +++ b/tests/test_ev_operator.py @@ -57,7 +57,7 @@ def test_quad_ker(monkeypatch): u=0, order=0, mode0="ns+", - mode1=None, + mode1="", method="", is_log=is_log, logx=0.0, @@ -109,7 +109,7 @@ def test_quad_ker(monkeypatch): u=0, order=0, mode0="ns+", - mode1=None, + mode1="", method="", is_log=True, logx=0.0, @@ -247,7 +247,7 @@ def quad_ker_pegasus( int_disp = InterpolatorDispatcher(xgrid, 1, True) order = 1 mode0 = "ns+" - mode1 = None + mode1 = "" method = "" logxs = np.log(int_disp.xgrid_raw) a1 = 1 diff --git a/tests/test_matching.py b/tests/test_matching.py index 02f1d0d04..dacab976b 100644 --- a/tests/test_matching.py +++ b/tests/test_matching.py @@ -3,6 +3,7 @@ import numpy as np from numpy.testing import assert_almost_equal +from eko import basis_rotation as br from eko import member from eko.matching_conditions import MatchingCondition @@ -17,17 +18,19 @@ class TestMatchingCondition: def mkOME(self): ome = {} - for key in ["qq", "qg", "gq", "gg", "Hq", "Hg"]: - ome.update({f"S_{key}": mkOM(self.shape)}) - if "g" not in key: - ome.update({f"NS_{key}": mkOM(self.shape)}) + for key in [ + *br.singlet_labels, + ("H+", "S"), + ("H+", "g"), + ("V", "V"), + ("H-", "V"), + ]: + ome.update({key: mkOM(self.shape)}) return ome def update_intrinsic_OME(self, ome): - for key in ["HH", "qH", "gH"]: - ome.update({f"S_{key}": mkOM(self.shape)}) - if "g" not in key: - ome.update({f"NS_{key}": mkOM(self.shape)}) + for key in [("H+", "H+"), ("H-", "H-"), ("V", "H-"), ("S", "H+"), ("g", "H+")]: + ome.update({key: mkOM(self.shape)}) def test_split_ad_to_evol_map(self): ome = self.mkOME() @@ -54,7 +57,7 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( a.op_members[member.MemberName("V.V")].value, - ome["NS_qq"].value, + ome[("V", "V")].value, ) # # if alpha is zero, nothing non-trivial should happen b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, []) @@ -107,13 +110,13 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( d.op_members[member.MemberName("b-.V")].value, - ome["NS_Hq"].value, + ome[("H-", "V")].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.S")].value, - ome["S_Hq"].value, + ome[("H+", "S")].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.g")].value, - ome["S_Hg"].value, + ome[("H+", "g")].value, ) diff --git a/tests/test_ome.py b/tests/test_ome.py index 9cd9c2105..109d04c2a 100644 --- a/tests/test_ome.py +++ b/tests/test_ome.py @@ -93,7 +93,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode="NS_qq", + mode0="V", + mode1="V", is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -106,7 +107,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode="S_qq", + mode0="S", + mode1="S", is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -119,7 +121,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode="S_qg", + mode0="S", + mode1="g", is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -131,12 +134,13 @@ def test_quad_ker(monkeypatch): np.testing.assert_allclose(res_s, 0.0) # test expanded intrisic inverse kernels - labels = ["NS_qq", *singlet_labels] + labels = [("V", "V"), *singlet_labels] for label in labels: res_ns = quad_ker( u=0, order=2, - mode=label, + mode0=label[0], + mode1=label[1], is_log=True, logx=0.0, areas=np.zeros(3), @@ -153,21 +157,22 @@ def test_quad_ker(monkeypatch): # test exact intrisic inverse kernel labels.extend( [ - "S_Hq", - "S_Hg", - "S_HH", - "S_qH", - "S_gH", - "NS_qH", - "NS_HH", - "NS_Hq", + ("H+", "S"), + ("H+", "g"), + ("H+", "H+"), + ("S", "H+"), + ("g", "H+"), + ("V", "H-"), + ("H-", "H-"), + ("H-", "V"), ] ) for label in labels: res_ns = quad_ker( u=0, order=2, - mode=label, + mode0=label[0], + mode1=label[1], is_log=True, logx=0.0, areas=np.zeros(3), @@ -185,7 +190,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode="NS_qq", + mode0="V", + mode1="V", is_log=True, logx=0.0, areas=np.array([0.01, 0.1, 1.0]), @@ -245,13 +251,19 @@ def test_labels(self): ) o = OperatorMatrixElement(g.config, g.managers, is_backward=False) labels = o.labels() - test_labels = ["NS_qq", "NS_Hq"] + test_labels = [("V", "V"), ("H-", "V")] for l in test_labels: if skip_ns: assert l not in labels else: assert l in labels - test_labels = ["S_qq", "S_Hq", "S_gg", "S_Hg", "S_gH"] + test_labels = [ + ("S", "S"), + ("H+", "S"), + ("g", "g"), + ("H+", "g"), + ("g", "H+"), + ] for l in test_labels: if skip_singlet: assert l not in labels @@ -280,24 +292,28 @@ def test_compute_lo(self): o = OperatorMatrixElement(g.config, g.managers, is_backward=False) o.compute(self.theory_card["mb"] ** 2, L=0, is_msbar=False) - dim = o.ome_members["NS_qq"].value.shape - for idx in ["S", "NS"]: + dim = o.ome_members[("V", "V")].value.shape + for indices in [("S", "H+"), ("V", "H-")]: np.testing.assert_allclose( - o.ome_members[f"{idx}_qq"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(indices[0], indices[0])].value, np.eye(dim[0]), atol=1e-8 ) np.testing.assert_allclose( - o.ome_members[f"{idx}_HH"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(indices[1], indices[1])].value, np.eye(dim[0]), atol=1e-8 + ) + np.testing.assert_allclose( + o.ome_members[(indices[0], indices[1])].value, np.zeros(dim) + ) + np.testing.assert_allclose( + o.ome_members[(indices[1], indices[0])].value, np.zeros(dim) ) - np.testing.assert_allclose(o.ome_members[f"{idx}_qH"].value, np.zeros(dim)) - np.testing.assert_allclose(o.ome_members[f"{idx}_Hq"].value, np.zeros(dim)) np.testing.assert_allclose( - o.ome_members["S_gg"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[("g", "g")].value, np.eye(dim[0]), atol=1e-8 ) np.testing.assert_allclose( - o.ome_members["S_qg"].value, o.ome_members["S_gq"].value + o.ome_members["S", "g"].value, o.ome_members[("g", "S")].value ) np.testing.assert_allclose( - o.ome_members["S_Hg"].value, o.ome_members["S_gH"].value + o.ome_members[("H+", "g")].value, o.ome_members[("g", "H+")].value ) def test_compute_nlo(self): @@ -326,20 +342,20 @@ def test_compute_nlo(self): dim = len(operators_card["interpolation_xgrid"]) shape = (dim, dim) - for idx in ["S", "NS"]: - assert o.ome_members[f"{idx}_qq"].value.shape == shape - assert o.ome_members[f"{idx}_HH"].value.shape == shape - assert o.ome_members[f"{idx}_qH"].value.shape == shape - assert o.ome_members[f"{idx}_Hq"].value.shape == shape + for indices in [("S", "H+"), ("V", "H-")]: + assert o.ome_members[(indices[0], indices[0])].value.shape == shape + assert o.ome_members[(indices[1], indices[1])].value.shape == shape + assert o.ome_members[(indices[0], indices[1])].value.shape == shape + assert o.ome_members[(indices[1], indices[0])].value.shape == shape np.testing.assert_allclose( - o.ome_members[f"{idx}_qH"].value, np.zeros(shape) + o.ome_members[(indices[0], indices[1])].value, np.zeros(shape) ) np.testing.assert_allclose( - o.ome_members[f"{idx}_Hq"].value, np.zeros(shape) + o.ome_members[(indices[1], indices[0])].value, np.zeros(shape) ) - assert o.ome_members["S_gg"].value.shape == shape + assert o.ome_members[("g", "g")].value.shape == shape np.testing.assert_allclose( - o.ome_members["S_qg"].value, o.ome_members["S_gq"].value + o.ome_members[("S", "g")].value, o.ome_members[("g", "S")].value ) - assert o.ome_members["S_Hg"].value.shape == shape - assert o.ome_members["S_gH"].value.shape == shape + assert o.ome_members[("H+", "g")].value.shape == shape + assert o.ome_members[("g", "H+")].value.shape == shape From 615c9a48a10ed2196562c42615a33ea7a71abda5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 24 Feb 2022 17:21:57 +0100 Subject: [PATCH 07/31] Fix heading --- doc/source/theory/pQCD.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 3a36b84cd..0f0e44793 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -48,7 +48,7 @@ is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, whil same expression for POLE masses is reported in Appendix A. QCD Splitting Functions -------------------- +----------------------- In the case in which only the |QCD| corrections are considered, the Altarelli-Parisi splitting kernels can be expanded in powers of the strong coupling :math:`a_s(\mu^2)` and are given by :cite:`Moch:2004pa,Vogt:2004mw` From e9bd255efb20897f66b9ddf92998774e22a15d6b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 25 Feb 2022 15:56:40 +0100 Subject: [PATCH 08/31] Replace labels with integers --- src/eko/anomalous_dimensions/__init__.py | 15 ++++--- src/eko/basis_rotation.py | 38 +++++++++++++----- src/eko/evolution_operator/__init__.py | 50 +++++++++++++----------- src/eko/evolution_operator/physical.py | 18 +++++---- src/eko/matching_conditions/__init__.py | 32 ++++++++------- tests/test_ad.py | 15 +++++-- tests/test_basis_rotation.py | 6 +-- tests/test_ev_operator.py | 42 +++++++++++--------- 8 files changed, 130 insertions(+), 86 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index b708e3c68..811d7cc85 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -20,6 +20,7 @@ import numba as nb import numpy as np +from .. import basis_rotation as br from . import as1, harmonics, nlo, nnlo @@ -71,7 +72,7 @@ def exp_singlet(gamma_S): return exp, lambda_p, lambda_m, e_p, e_m -@nb.njit("c16[:](u1,string,c16,u1)", cache=True) +@nb.njit("c16[:](u1,u1,c16,u1)", cache=True) def gamma_ns(order, mode, n, nf): r""" Computes the tower of the non-singlet anomalous dimensions @@ -109,21 +110,23 @@ def gamma_ns(order, mode, n, nf): # NLO and beyond if order >= 1: # TODO: pass the necessary harmonics to nlo gammas - if mode == "+": + if mode == 10101: gamma_ns_1 = nlo.gamma_nsp_1(n, nf) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here - elif mode in ["-", "V"]: + elif mode in [10201, 10200]: gamma_ns_1 = nlo.gamma_nsm_1(n, nf) + else: + raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order >= 2: sx = np.append(sx, harmonics.harmonic_S2(n)) sx = np.append(sx, harmonics.harmonic_S3(n)) - if mode == "+": + if mode == 10101: gamma_ns_2 = -nnlo.gamma_nsp_2(n, nf, sx) - elif mode == "-": + elif mode == 10201: gamma_ns_2 = -nnlo.gamma_nsm_2(n, nf, sx) - elif mode == "V": + elif mode == 10200: gamma_ns_2 = -nnlo.gamma_nsv_2(n, nf, sx) gamma_ns[2] = gamma_ns_2 return gamma_ns diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index fc63f36e1..d3c89a769 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -69,13 +69,21 @@ ) """|pid| representation of :data:`evol_basis`.""" -singlet_labels = (("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")) -non_singlet_labels = (("ns-", ""), ("ns+", ""), ("nsV", "")) +non_singlet_pids_map = {"ns-": 10201, "ns+": 10101, "nsV": 10200} + +singlet_labels = ((100, 100), (100, 21), (21, 100), (21, 21)) +non_singlet_labels = ( + (non_singlet_pids_map["ns-"], 0), + (non_singlet_pids_map["ns+"], 0), + (non_singlet_pids_map["nsV"], 0), +) full_labels = (*singlet_labels, *non_singlet_labels) anomalous_dimensions_basis = full_labels r""" Sorted elements in Anomalous Dimensions Basis as :obj:`str`. """ +matching_hplus_pid = 90 +matching_hminus_pid = 91 # Tranformation from physical basis to QCD evolution basis rotate_flavor_to_evolution = np.array( @@ -101,13 +109,25 @@ """ map_ad_to_evolution = { - ("S", "S"): ["S.S"], - ("S", "g"): ["S.g"], - ("g", "S"): ["g.S"], - ("g", "g"): ["g.g"], - ("nsV", ""): ["V.V"], - ("ns+", ""): ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], - ("ns-", ""): ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], + (100, 100): ["S.S"], + (100, 21): ["S.g"], + (21, 100): ["g.S"], + (21, 21): ["g.g"], + (non_singlet_pids_map["nsV"], 0): ["V.V"], + (non_singlet_pids_map["ns+"], 0): [ + "T3.T3", + "T8.T8", + "T15.T15", + "T24.T24", + "T35.T35", + ], + (non_singlet_pids_map["ns-"], 0): [ + "V3.V3", + "V8.V8", + "V15.V15", + "V24.V24", + "V35.V35", + ], } """ Map anomalous dimension sectors' names to their members diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index c12cb073b..d5f10a46a 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -13,8 +13,8 @@ from scipy import integrate from .. import anomalous_dimensions as ad +from .. import basis_rotation as br from .. import beta, interpolation, mellin -from ..basis_rotation import full_labels, singlet_labels from ..kernels import non_singlet as ns from ..kernels import singlet as s from ..member import OpMember @@ -22,7 +22,7 @@ logger = logging.getLogger(__name__) -@nb.njit("c16[:](u1,string,c16,u1,f8)", cache=True) +@nb.njit("c16[:](u1,u1,c16,u1,f8)", cache=True) def gamma_ns_fact(order, mode, n, nf, L): """ Adjust the anomalous dimensions with the scale variations. @@ -45,7 +45,7 @@ def gamma_ns_fact(order, mode, n, nf, L): gamma_ns : numpy.ndarray adjusted non-singlet anomalous dimensions """ - gamma_ns = ad.gamma_ns(order, mode[-1], n, nf) + gamma_ns = ad.gamma_ns(order, mode, n, nf) # since we are modifying *in-place* be carefull, that the order matters! # and indeed, we need to adjust the high elements first if order >= 2: @@ -97,7 +97,7 @@ def gamma_singlet_fact(order, n, nf, L): return gamma_singlet -@nb.njit("f8(f8,u1,string,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) +@nb.njit("f8(f8,u1,u1,u1,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) def quad_ker( u, order, @@ -151,7 +151,7 @@ def quad_ker( ker : float evaluated integration kernel """ - is_singlet = len(mode1) > 0 + is_singlet = mode1 != 0 # get transformation to N integral if logx == 0.0: return 0.0 @@ -176,8 +176,8 @@ def quad_ker( order, method, gamma_singlet, a1, a0, nf, ev_op_iterations, ev_op_max_order ) # select element of matrix - k = 0 if mode0 == "S" else 1 - l = 0 if mode1 == "S" else 1 + k = 0 if mode0 == 100 else 1 + l = 0 if mode1 == 100 else 1 ker = ker[k, l] else: gamma_ns = gamma_ns_fact(order, mode0, n, nf, L) @@ -243,16 +243,16 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append(("ns+", "")) + labels.append((br.non_singlet_pids_map["ns+"], 0)) if order >= 1: # - becomes different starting from NLO - labels.append(("ns-", "")) + labels.append((br.non_singlet_pids_map["ns-"], 0)) if order >= 2: # v also becomes different starting from NNLO - labels.append(("nsV", "")) + labels.append((br.non_singlet_pids_map["nsV"], 0)) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") else: - labels.extend([("S", "S"), ("S", "g"), ("g", "S"), ("g", "g")]) + labels.extend(br.singlet_labels) return labels def compute(self): @@ -265,10 +265,10 @@ def compute(self): labels = self.labels() eye = OpMember(np.eye(grid_size), np.zeros((grid_size, grid_size))) zero = OpMember(*[np.zeros((grid_size, grid_size))] * 2) - for n in full_labels: + for n in br.full_labels: if n in labels: # off diag singlet are zero - if n in [("S", "g"), ("g", "S")]: + if n in br.singlet_labels and n[0] != n[1]: self.op_members[n] = zero.copy() else: self.op_members[n] = eye.copy() @@ -358,16 +358,20 @@ def copy_ns_ops(self): order = self.config["order"] if order == 0: # in LO +=-=v for label in ["nsV", "ns-"]: - self.op_members[(label, "")].value = self.op_members[ - ("ns+", "") + self.op_members[ + (br.non_singlet_pids_map[label], 0) + ].value = self.op_members[ + (br.non_singlet_pids_map["ns+"], 0) ].value.copy() - self.op_members[(label, "")].error = self.op_members[ - ("ns+", "") + self.op_members[ + (br.non_singlet_pids_map[label], 0) + ].error = self.op_members[ + (br.non_singlet_pids_map["ns+"], 0) ].error.copy() elif order == 1: # in NLO -=v - self.op_members[("nsV", "")].value = self.op_members[ - ("ns-", "") - ].value.copy() - self.op_members[("nsV", "")].error = self.op_members[ - ("ns-", "") - ].error.copy() + self.op_members[ + (br.non_singlet_pids_map["nsV"], 0) + ].value = self.op_members[(br.non_singlet_pids_map["ns-"], 0)].value.copy() + self.op_members[ + (br.non_singlet_pids_map["nsV"], 0) + ].error = self.op_members[(br.non_singlet_pids_map["ns-"], 0)].error.copy() diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index a8bc05709..aed503520 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -48,21 +48,23 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range): """ # constant elements m = { - "S.S": op_members[("S", "S")], - "S.g": op_members[("S", "g")], - "g.g": op_members[("g", "g")], - "g.S": op_members[("g", "S")], - "V.V": op_members[("nsV", "")], + "S.S": op_members[(100, 100)], + "S.g": op_members[(100, 21)], + "g.g": op_members[(21, 21)], + "g.S": op_members[(21, 100)], + "V.V": op_members[(br.non_singlet_pids_map["nsV"], 0)], } # add elements which are already active for f in range(2, nf + 1): n = f**2 - 1 - m[f"V{n}.V{n}"] = op_members[("ns-", "")] - m[f"T{n}.T{n}"] = op_members[("ns+", "")] + m[f"V{n}.V{n}"] = op_members[(br.non_singlet_pids_map["ns-"], 0)] + m[f"T{n}.T{n}"] = op_members[(br.non_singlet_pids_map["ns+"], 0)] # deal with intrinsic heavy quark PDFs if intrinsic_range is not None: hqfl = "cbt" - op_id = member.OpMember.id_like(op_members[("nsV", "")]) + op_id = member.OpMember.id_like( + op_members[(br.non_singlet_pids_map["nsV"], 0)] + ) for intr_fl in intrinsic_range: if intr_fl <= nf: # light quarks are not intrinsic continue diff --git a/src/eko/matching_conditions/__init__.py b/src/eko/matching_conditions/__init__.py index ce15e7a75..b76d7a1aa 100644 --- a/src/eko/matching_conditions/__init__.py +++ b/src/eko/matching_conditions/__init__.py @@ -44,11 +44,11 @@ def split_ad_to_evol_map( """ m = { - "S.S": ome_members[("S", "S")], - "S.g": ome_members[("S", "g")], # This is always zero for the time being - "g.S": ome_members[("g", "S")], - "g.g": ome_members[("g", "g")], - "V.V": ome_members[("V", "V")], + "S.S": ome_members[(100, 100)], + "S.g": ome_members[(100, 21)], # This is always zero for the time being + "g.S": ome_members[(21, 100)], + "g.g": ome_members[(21, 21)], + "V.V": ome_members[(200, 200)], } # add elements which are already active @@ -61,15 +61,15 @@ def split_ad_to_evol_map( hq = br.quark_names[nf] m.update( { - f"{hq}-.V": ome_members[("H-", "V")], - f"{hq}+.S": ome_members[("H+", "S")], - f"{hq}+.g": ome_members[("H+", "g")], + f"{hq}-.V": ome_members[(br.matching_hminus_pid, 200)], + f"{hq}+.S": ome_members[(br.matching_hplus_pid, 100)], + f"{hq}+.g": ome_members[(br.matching_hplus_pid, 21)], } ) # intrinsic matching if len(intrinsic_range) != 0: - op_id = member.OpMember.id_like(ome_members[("V", "V")]) + op_id = member.OpMember.id_like(ome_members[(200, 200)]) for intr_fl in intrinsic_range: ihq = br.quark_names[intr_fl - 1] # find name if intr_fl > nf + 1: @@ -80,11 +80,15 @@ def split_ad_to_evol_map( # match the missing contribution from h+ and h- m.update( { - f"{ihq}+.{ihq}+": ome_members[("H+", "H+")], - # f"S.{ihq}+": ome_members[("S", "H+")], - f"g.{ihq}+": ome_members[("g", "H+")], - f"{ihq}-.{ihq}-": ome_members[("H-", "H-")], - # f"V.{ihq}-": ome_members[("V", "H-")], + f"{ihq}+.{ihq}+": ome_members[ + (br.matching_hplus_pid, br.matching_hplus_pid) + ], + # f"S.{ihq}+": ome_members[(100, br.matching_hplus_pid)], + f"g.{ihq}+": ome_members[(21, br.matching_hplus_pid)], + f"{ihq}-.{ihq}-": ome_members[ + (br.matching_hminus_pid, br.matching_hminus_pid) + ], + # f"V.{ihq}-": ome_members[(200, br.matching_hminus_pid)], } ) return cls.promote_names(m, q2_thr) diff --git a/tests/test_ad.py b/tests/test_ad.py index decbb2943..3255e1751 100644 --- a/tests/test_ad.py +++ b/tests/test_ad.py @@ -4,6 +4,7 @@ from numpy.testing import assert_allclose, assert_almost_equal from eko import anomalous_dimensions as ad +from eko import basis_rotation as br from eko.anomalous_dimensions import as1 as ad_lo from eko.anomalous_dimensions import harmonics @@ -53,9 +54,15 @@ def test_eigensystem_gamma_singlet_projectors_EV(): def test_gamma_ns(): nf = 3 # LO - assert_almost_equal(ad.gamma_ns(2, "+", 1, nf)[0], 0.0) + assert_almost_equal(ad.gamma_ns(2, br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0) # NLO - assert_allclose(ad.gamma_ns(1, "-", 1, nf), np.zeros(2), atol=2e-6) + assert_allclose( + ad.gamma_ns(1, br.non_singlet_pids_map["ns-"], 1, nf), np.zeros(2), atol=2e-6 + ) # NNLO - assert_allclose(ad.gamma_ns(2, "-", 1, nf), np.zeros(3), atol=2e-4) - assert_allclose(ad.gamma_ns(2, "V", 1, nf), np.zeros(3), atol=8e-4) + assert_allclose( + ad.gamma_ns(2, br.non_singlet_pids_map["ns-"], 1, nf), np.zeros(3), atol=2e-4 + ) + assert_allclose( + ad.gamma_ns(2, br.non_singlet_pids_map["nsV"], 1, nf), np.zeros(3), atol=8e-4 + ) diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py index c102af936..faf320b73 100644 --- a/tests/test_basis_rotation.py +++ b/tests/test_basis_rotation.py @@ -10,19 +10,19 @@ def test_ad_projector(): g = br.rotate_flavor_to_evolution[2] v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")] - s_to_s = br.ad_projector(("S", "S"), nf=6) + s_to_s = br.ad_projector((100, 100), nf=6) np.testing.assert_allclose(s @ s_to_s, s) np.testing.assert_allclose(g @ s_to_s, 0.0) np.testing.assert_allclose(v3 @ s_to_s, 0.0) - g_to_s = br.ad_projector(("g", "S"), nf=6) + g_to_s = br.ad_projector((21, 100), nf=6) np.testing.assert_allclose(s @ g_to_s, 0.0) np.testing.assert_allclose(g @ g_to_s, s) np.testing.assert_allclose(v3 @ g_to_s, 0.0) - ns_m = br.ad_projector(("ns-", ""), nf=6) + ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6) np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15) np.testing.assert_allclose(g @ ns_m, 0.0) diff --git a/tests/test_ev_operator.py b/tests/test_ev_operator.py index 6e5df0b18..055e267d3 100644 --- a/tests/test_ev_operator.py +++ b/tests/test_ev_operator.py @@ -18,13 +18,13 @@ def test_gamma_ns_fact(monkeypatch): gamma_ns = np.array([1.0, 0.5, 0.25]) monkeypatch.setattr(ad, "gamma_ns", lambda *args: gamma_ns.copy()) - gamma_ns_LO_0 = gamma_ns_fact(0, "ns+", 1, 3, 0) + gamma_ns_LO_0 = gamma_ns_fact(0, br.non_singlet_pids_map["ns+"], 1, 3, 0) np.testing.assert_allclose(gamma_ns_LO_0, gamma_ns) - gamma_ns_LO_1 = gamma_ns_fact(0, "ns+", 1, 3, 1) + gamma_ns_LO_1 = gamma_ns_fact(0, br.non_singlet_pids_map["ns+"], 1, 3, 1) np.testing.assert_allclose(gamma_ns_LO_1, gamma_ns) - gamma_ns_NLO_1 = gamma_ns_fact(1, "ns+", 1, 3, 1) + gamma_ns_NLO_1 = gamma_ns_fact(1, br.non_singlet_pids_map["ns+"], 1, 3, 1) assert gamma_ns_NLO_1[1] < gamma_ns[1] - gamma_ns_NNLO_1 = gamma_ns_fact(2, "ns+", 1, 3, 1) + gamma_ns_NNLO_1 = gamma_ns_fact(2, br.non_singlet_pids_map["ns+"], 1, 3, 1) assert gamma_ns_NNLO_1[2] - gamma_ns[2] == 8.0 @@ -56,8 +56,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode0="ns+", - mode1="", + mode0=br.non_singlet_pids_map["ns+"], + mode1=0, method="", is_log=is_log, logx=0.0, @@ -73,8 +73,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode0="S", - mode1="S", + mode0=100, + mode1=100, method="", is_log=is_log, logx=1.0, @@ -90,8 +90,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode0="S", - mode1="g", + mode0=100, + mode1=21, method="", is_log=is_log, logx=0.0, @@ -108,8 +108,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode0="ns+", - mode1="", + mode0=br.non_singlet_pids_map["ns+"], + mode1=0, method="", is_log=True, logx=0.0, @@ -193,21 +193,25 @@ def test_compute(self, monkeypatch): ) # LO o.compute() - assert ("ns-", "") in o.op_members + assert (br.non_singlet_pids_map["ns-"], 0) in o.op_members np.testing.assert_allclose( - o.op_members[("ns-", "")].value, o.op_members[("ns+", "")].value + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, ) np.testing.assert_allclose( - o.op_members[("nsV", "")].value, o.op_members[("ns+", "")].value + o.op_members[(br.non_singlet_pids_map["nsV"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, ) # NLO o.config["order"] = 1 o.compute() assert not np.allclose( - o.op_members[("ns+", "")].value, o.op_members[("ns-", "")].value + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, ) np.testing.assert_allclose( - o.op_members[("nsV", "")].value, o.op_members[("ns-", "")].value + o.op_members[(br.non_singlet_pids_map["nsV"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, ) # unity operators @@ -246,8 +250,8 @@ def quad_ker_pegasus( xgrid = np.geomspace(1e-7, 1, 10) int_disp = InterpolatorDispatcher(xgrid, 1, True) order = 1 - mode0 = "ns+" - mode1 = "" + mode0 = br.non_singlet_pids_map["ns+"] + mode1 = 0 method = "" logxs = np.log(int_disp.xgrid_raw) a1 = 1 From 24326745293ba492528d2dbab06c7179947e7bcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 25 Feb 2022 16:11:54 +0100 Subject: [PATCH 09/31] Replace string with pids in matching --- .../operator_matrix_element.py | 49 ++++++++----- tests/test_matching.py | 24 ++++--- tests/test_ome.py | 69 ++++++++++--------- 3 files changed, 80 insertions(+), 62 deletions(-) diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index aa068bfdf..56c647915 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -11,9 +11,9 @@ import numpy as np from scipy import integrate +from .. import basis_rotation as br from .. import interpolation, mellin from ..anomalous_dimensions import harmonics -from ..basis_rotation import singlet_labels from ..member import OpMember from . import nlo, nnlo @@ -145,7 +145,7 @@ def build_ome(A, order, a_s, backward_method): return ome -@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) +@nb.njit("f8(f8,u1,u1,u1,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) def quad_ker( u, order, mode0, mode1, is_log, logx, areas, a_s, L, backward_method, is_msbar ): @@ -179,15 +179,15 @@ def quad_ker( ker : float evaluated integration kernel """ - is_singlet = mode0 in ["S", "g", "H+"] + is_singlet = mode0 in [100, 21, 90] # get transformation to N integral r = 0.4 * 16.0 / (1.0 - logx) if is_singlet: o = 1.0 - indices = {"g": 0, "S": 1, "H+": 2} + indices = {21: 0, 100: 1, 90: 2} else: o = 0.0 - indices = {"V": 0, "H-": 1} + indices = {200: 0, 91: 1} n = mellin.Talbot_path(u, r, o) jac = mellin.Talbot_jac(u, r, o) @@ -269,10 +269,10 @@ def labels(self): if self.config["debug_skip_non_singlet"]: logger.warning("Matching: skipping non-singlet sector") else: - labels.extend([("V", "V"), ("H-", "V")]) + labels.extend([(200, 200), (br.matching_hminus_pid, 200)]) if self.is_intrinsic or self.backward_method != "": # intrisic labels, which are not zero at NLO - labels.append(("H-", "H-")) + labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) # if self.backward_method == "exact": # # this contribution starts at NNLO, we don't have it for the moment # labels.append("NS_qH") @@ -281,9 +281,20 @@ def labels(self): if self.config["debug_skip_singlet"]: logger.warning("Matching: skipping singlet sector") else: - labels.extend([*singlet_labels, ("H+", "g"), ("H+", "S")]) + labels.extend( + [ + *br.singlet_labels, + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, 100), + ] + ) if self.is_intrinsic or self.backward_method != "": - labels.extend([("g", "H+"), ("H+", "H+")]) + labels.extend( + [ + (21, br.matching_hplus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + ] + ) # if self.backward_method == "exact": # labels.extend(["S_qH"]) return labels @@ -376,11 +387,11 @@ def copy_ome(self): grid_size = len(self.int_disp.xgrid) # basic labels skipped with skip debug for label in [ - ("V", "V"), - ("H+", "g"), - ("H+", "S"), - ("H-", "V"), - *singlet_labels, + (200, 200), + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, 100), + (br.matching_hminus_pid, 200), + *br.singlet_labels, ]: if label not in self.ome_members: self.ome_members[label] = OpMember( @@ -390,11 +401,11 @@ def copy_ome(self): # intrinsic labels not computed yet if self.is_intrinsic: for label in [ - ("S", "H+"), - ("V", "H-"), - ("H-", "H-"), - ("H+", "H+"), - ("g", "H+"), + (100, br.matching_hplus_pid), + (200, br.matching_hminus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + (21, br.matching_hplus_pid), ]: if label not in self.ome_members: self.ome_members[label] = OpMember( diff --git a/tests/test_matching.py b/tests/test_matching.py index dacab976b..6aecbc3e8 100644 --- a/tests/test_matching.py +++ b/tests/test_matching.py @@ -20,16 +20,22 @@ def mkOME(self): ome = {} for key in [ *br.singlet_labels, - ("H+", "S"), - ("H+", "g"), - ("V", "V"), - ("H-", "V"), + (br.matching_hplus_pid, 100), + (br.matching_hplus_pid, 21), + (200, 200), + (br.matching_hminus_pid, 200), ]: ome.update({key: mkOM(self.shape)}) return ome def update_intrinsic_OME(self, ome): - for key in [("H+", "H+"), ("H-", "H-"), ("V", "H-"), ("S", "H+"), ("g", "H+")]: + for key in [ + (br.matching_hplus_pid, br.matching_hplus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (200, br.matching_hminus_pid), + (100, br.matching_hplus_pid), + (21, br.matching_hplus_pid), + ]: ome.update({key: mkOM(self.shape)}) def test_split_ad_to_evol_map(self): @@ -57,7 +63,7 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( a.op_members[member.MemberName("V.V")].value, - ome[("V", "V")].value, + ome[(200, 200)].value, ) # # if alpha is zero, nothing non-trivial should happen b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, []) @@ -110,13 +116,13 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( d.op_members[member.MemberName("b-.V")].value, - ome[("H-", "V")].value, + ome[(br.matching_hminus_pid, 200)].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.S")].value, - ome[("H+", "S")].value, + ome[(br.matching_hplus_pid, 100)].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.g")].value, - ome[("H+", "g")].value, + ome[(br.matching_hplus_pid, 21)].value, ) diff --git a/tests/test_ome.py b/tests/test_ome.py index 109d04c2a..bae1c155c 100644 --- a/tests/test_ome.py +++ b/tests/test_ome.py @@ -4,8 +4,8 @@ import numpy as np +from eko import basis_rotation as br from eko import interpolation, mellin -from eko.basis_rotation import singlet_labels from eko.evolution_operator.grid import OperatorGrid from eko.interpolation import InterpolatorDispatcher from eko.matching_conditions.operator_matrix_element import ( @@ -93,8 +93,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode0="V", - mode1="V", + mode0=200, + mode1=200, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -107,8 +107,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode0="S", - mode1="S", + mode0=100, + mode1=100, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -121,8 +121,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode0="S", - mode1="g", + mode0=100, + mode1=21, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -134,7 +134,7 @@ def test_quad_ker(monkeypatch): np.testing.assert_allclose(res_s, 0.0) # test expanded intrisic inverse kernels - labels = [("V", "V"), *singlet_labels] + labels = [(200, 200), *br.singlet_labels] for label in labels: res_ns = quad_ker( u=0, @@ -157,14 +157,14 @@ def test_quad_ker(monkeypatch): # test exact intrisic inverse kernel labels.extend( [ - ("H+", "S"), - ("H+", "g"), - ("H+", "H+"), - ("S", "H+"), - ("g", "H+"), - ("V", "H-"), - ("H-", "H-"), - ("H-", "V"), + (br.matching_hplus_pid, 100), + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, br.matching_hplus_pid), + (100, br.matching_hplus_pid), + (21, br.matching_hplus_pid), + (200, br.matching_hminus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (br.matching_hminus_pid, 200), ] ) for label in labels: @@ -190,8 +190,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode0="V", - mode1="V", + mode0=200, + mode1=200, is_log=True, logx=0.0, areas=np.array([0.01, 0.1, 1.0]), @@ -251,18 +251,18 @@ def test_labels(self): ) o = OperatorMatrixElement(g.config, g.managers, is_backward=False) labels = o.labels() - test_labels = [("V", "V"), ("H-", "V")] + test_labels = [(200, 200), (br.matching_hminus_pid, 200)] for l in test_labels: if skip_ns: assert l not in labels else: assert l in labels test_labels = [ - ("S", "S"), - ("H+", "S"), - ("g", "g"), - ("H+", "g"), - ("g", "H+"), + (100, 100), + (br.matching_hplus_pid, 100), + (21, 21), + (br.matching_hplus_pid, 21), + (21, br.matching_hplus_pid), ] for l in test_labels: if skip_singlet: @@ -292,8 +292,8 @@ def test_compute_lo(self): o = OperatorMatrixElement(g.config, g.managers, is_backward=False) o.compute(self.theory_card["mb"] ** 2, L=0, is_msbar=False) - dim = o.ome_members[("V", "V")].value.shape - for indices in [("S", "H+"), ("V", "H-")]: + dim = o.ome_members[(200, 200)].value.shape + for indices in [(100, br.matching_hplus_pid), (200, br.matching_hminus_pid)]: np.testing.assert_allclose( o.ome_members[(indices[0], indices[0])].value, np.eye(dim[0]), atol=1e-8 ) @@ -307,13 +307,14 @@ def test_compute_lo(self): o.ome_members[(indices[1], indices[0])].value, np.zeros(dim) ) np.testing.assert_allclose( - o.ome_members[("g", "g")].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(21, 21)].value, np.eye(dim[0]), atol=1e-8 ) np.testing.assert_allclose( - o.ome_members["S", "g"].value, o.ome_members[("g", "S")].value + o.ome_members[100, 21].value, o.ome_members[(21, 100)].value ) np.testing.assert_allclose( - o.ome_members[("H+", "g")].value, o.ome_members[("g", "H+")].value + o.ome_members[(br.matching_hplus_pid, 21)].value, + o.ome_members[(21, br.matching_hplus_pid)].value, ) def test_compute_nlo(self): @@ -342,7 +343,7 @@ def test_compute_nlo(self): dim = len(operators_card["interpolation_xgrid"]) shape = (dim, dim) - for indices in [("S", "H+"), ("V", "H-")]: + for indices in [(100, br.matching_hplus_pid), (200, br.matching_hminus_pid)]: assert o.ome_members[(indices[0], indices[0])].value.shape == shape assert o.ome_members[(indices[1], indices[1])].value.shape == shape assert o.ome_members[(indices[0], indices[1])].value.shape == shape @@ -353,9 +354,9 @@ def test_compute_nlo(self): np.testing.assert_allclose( o.ome_members[(indices[1], indices[0])].value, np.zeros(shape) ) - assert o.ome_members[("g", "g")].value.shape == shape + assert o.ome_members[(21, 21)].value.shape == shape np.testing.assert_allclose( - o.ome_members[("S", "g")].value, o.ome_members[("g", "S")].value + o.ome_members[(100, 21)].value, o.ome_members[(21, 100)].value ) - assert o.ome_members[("H+", "g")].value.shape == shape - assert o.ome_members[("g", "H+")].value.shape == shape + assert o.ome_members[(br.matching_hplus_pid, 21)].value.shape == shape + assert o.ome_members[(21, br.matching_hplus_pid)].value.shape == shape From 4addb86cb07c68a7eec75e0db4fad46c31785794 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 25 Feb 2022 16:53:10 +0100 Subject: [PATCH 10/31] Change mode to u2 --- benchmarks/lha_paper_bench.py | 6 +++--- src/eko/anomalous_dimensions/__init__.py | 3 ++- src/eko/evolution_operator/__init__.py | 4 ++-- .../operator_matrix_element.py | 2 +- tests/test_runner.py | 19 +++++++++++++++++++ 5 files changed, 27 insertions(+), 7 deletions(-) diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py index 3021bf0bb..b9a40ba4d 100644 --- a/benchmarks/lha_paper_bench.py +++ b/benchmarks/lha_paper_bench.py @@ -174,8 +174,8 @@ def skip_pdfs(theory): if __name__ == "__main__": - # obj = BenchmarkVFNS() - obj = BenchmarkFFNS() + obj = BenchmarkVFNS() + # obj = BenchmarkFFNS() - obj.benchmark_plain(0) + obj.benchmark_plain(2) # obj.benchmark_sv(1) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 811d7cc85..16468566c 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -72,7 +72,7 @@ def exp_singlet(gamma_S): return exp, lambda_p, lambda_m, e_p, e_m -@nb.njit("c16[:](u1,u1,c16,u1)", cache=True) +@nb.njit("c16[:](u1,u2,c16,u1)", cache=True) def gamma_ns(order, mode, n, nf): r""" Computes the tower of the non-singlet anomalous dimensions @@ -116,6 +116,7 @@ def gamma_ns(order, mode, n, nf): elif mode in [10201, 10200]: gamma_ns_1 = nlo.gamma_nsm_1(n, nf) else: + print(mode) raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index d5f10a46a..20adf6f93 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -22,7 +22,7 @@ logger = logging.getLogger(__name__) -@nb.njit("c16[:](u1,u1,c16,u1,f8)", cache=True) +@nb.njit("c16[:](u1,u2,c16,u1,f8)", cache=True) def gamma_ns_fact(order, mode, n, nf, L): """ Adjust the anomalous dimensions with the scale variations. @@ -97,7 +97,7 @@ def gamma_singlet_fact(order, n, nf, L): return gamma_singlet -@nb.njit("f8(f8,u1,u1,u1,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) +@nb.njit("f8(f8,u1,u2,u2,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1)", cache=True) def quad_ker( u, order, diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index 56c647915..c7fb1340a 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -145,7 +145,7 @@ def build_ome(A, order, a_s, backward_method): return ome -@nb.njit("f8(f8,u1,u1,u1,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) +@nb.njit("f8(f8,u1,u2,u2,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) def quad_ker( u, order, mode0, mode1, is_log, logx, areas, a_s, L, backward_method, is_msbar ): diff --git a/tests/test_runner.py b/tests/test_runner.py index b31f88331..ef33273ef 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -115,3 +115,22 @@ def check_shapes(o, txs, ixs, theory_card, operators_card): assert "operator_errors" in ops assert ops["operators"].shape == op_shape assert ops["operator_errors"].shape == op_shape + + +def test_vfns(): + # change targetbasis + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operators_card) + tc["kcThr"] = 1.0 + tc["kbThr"] = 1.0 + tc["PTO"] = 2 + oc["debug_skip_non_singlet"] = False + r = eko.runner.Runner(tc, oc) + o = r.get_output() + check_shapes( + o, + o["interpolation_xgrid"], + o["interpolation_xgrid"], + tc, + oc, + ) From 9b8c1310eab032542e6c42b11f0b7b64cfacdcb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 25 Feb 2022 17:05:35 +0100 Subject: [PATCH 11/31] Remove previous debug print --- src/eko/anomalous_dimensions/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 16468566c..78dd48fa6 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -116,7 +116,6 @@ def gamma_ns(order, mode, n, nf): elif mode in [10201, 10200]: gamma_ns_1 = nlo.gamma_nsm_1(n, nf) else: - print(mode) raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond From de4d4eb1630357d617fe80d4283261527585cda7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 25 Feb 2022 17:16:01 +0100 Subject: [PATCH 12/31] Fix benchmark_evol_to_unity.py --- tests/benchmark_evol_to_unity.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/tests/benchmark_evol_to_unity.py b/tests/benchmark_evol_to_unity.py index e9f33fbe7..6f0eacdf7 100644 --- a/tests/benchmark_evol_to_unity.py +++ b/tests/benchmark_evol_to_unity.py @@ -2,6 +2,7 @@ import numpy as np +from eko import basis_rotation as br from eko.evolution_operator import Operator from eko.evolution_operator.grid import OperatorGrid from eko.interpolation import InterpolatorDispatcher @@ -74,38 +75,39 @@ def test_operator_grid( o.compute() o_back.compute() - dim = o_back.op_members["NS_v"].value.shape - for k in ["NS_v", "NS_m", "NS_p"]: + dim = o_back.op_members[(br.non_singlet_pids_map["nsV"], 0)].value.shape + for k in ["nsV", "ns-", "ns+"]: np.testing.assert_allclose( - o.op_members[k].value @ o_back.op_members[k].value, + o.op_members[(br.non_singlet_pids_map[k], 0)].value + @ o_back.op_members[(br.non_singlet_pids_map[k], 0)].value, np.eye(dim[0]), atol=7e-2, ) # qq np.testing.assert_allclose( - o_back.op_members["S_qq"].value @ o.op_members["S_qq"].value - + o_back.op_members["S_qg"].value @ o.op_members["S_gq"].value, + o_back.op_members[(100, 100)].value @ o.op_members[(100, 100)].value + + o_back.op_members[(100, 21)].value @ o.op_members[(21, 100)].value, np.eye(dim[0]), atol=7e-2, ) # qg np.testing.assert_allclose( - o_back.op_members["S_qq"].value @ o.op_members["S_qg"].value - + o_back.op_members["S_qg"].value @ o.op_members["S_gg"].value, + o_back.op_members[(100, 100)].value @ o.op_members[(100, 21)].value + + o_back.op_members[(100, 21)].value @ o.op_members[(21, 21)].value, np.zeros(dim), atol=7e-4, ) # gg np.testing.assert_allclose( - o_back.op_members["S_gg"].value @ o.op_members["S_gg"].value - + o_back.op_members["S_gq"].value @ o.op_members["S_qg"].value, + o_back.op_members[(21, 21)].value @ o.op_members[(21, 21)].value + + o_back.op_members[(21, 100)].value @ o.op_members[(100, 21)].value, np.eye(dim[0]), atol=9e-2, ) # gq np.testing.assert_allclose( - o_back.op_members["S_gg"].value @ o.op_members["S_gq"].value - + o_back.op_members["S_gq"].value @ o.op_members["S_qq"].value, + o_back.op_members[(21, 21)].value @ o.op_members[(21, 100)].value + + o_back.op_members[(21, 100)].value @ o.op_members[(100, 100)].value, np.zeros(dim), atol=2e-3, ) @@ -122,7 +124,7 @@ def test_operator_grid( # ome.compute( q2, L) # ome_back.compute(q2, L) - # dim = ome.ome_members["S_qq"].value.shape + # dim = ome.ome_members[(100, 100)].value.shape # ome_tensor = np.zeros((3,3,dim[0],dim[0])) # ome_tensor_back = ome_tensor # idx_dict = dict(zip(["g", "q", "H"],[0,1,2])) From 0ff06cb9f0c1642b8a2cf1d6e3692c7aa646910b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 1 Mar 2022 16:24:51 +0100 Subject: [PATCH 13/31] Test aem1.py and NotImplementedError in anomalous_dimensions/__init__.py --- src/eko/constants.py | 4 ++-- tests/test_ad_lo.py | 36 +++++++++++++++++++++++++++++------- tests/test_kernels_ns.py | 3 +++ 3 files changed, 34 insertions(+), 9 deletions(-) diff --git a/src/eko/constants.py b/src/eko/constants.py index 3de4379b3..470168305 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -11,10 +11,10 @@ """the normalization of fundamental generators""" CA = float(NC) -"""second Casimir constant in the adjoint representation""" +"""second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3""" CF = float((NC * NC - 1.0) / (2.0 * NC)) -"""second Casimir constant in the fundamental representation""" +"""second Casimir constant in the fundamental representation - defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3""" eu2 = 4.0 / 9 """up quarks charge squared""" diff --git a/tests/test_ad_lo.py b/tests/test_ad_lo.py index 0dc9cb078..98b60e9e8 100644 --- a/tests/test_ad_lo.py +++ b/tests/test_ad_lo.py @@ -2,7 +2,8 @@ # Test LO splitting functions import numpy as np -import eko.anomalous_dimensions.as1 as ad_lo +import eko.anomalous_dimensions.aem1 as ad_lo_aem1 +import eko.anomalous_dimensions.as1 as ad_lo_as1 from eko.anomalous_dimensions import harmonics NF = 5 @@ -12,14 +13,18 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) s1 = harmonics.harmonic_S1(N) - np.testing.assert_almost_equal(ad_lo.gamma_ns_0(N, s1), 0) + np.testing.assert_almost_equal(ad_lo_as1.gamma_ns_0(N, s1), 0) def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) - np.testing.assert_almost_equal(ad_lo.gamma_ns_0(N, s1) + ad_lo.gamma_gq_0(N), 0) + np.testing.assert_almost_equal( + ad_lo_as1.gamma_ns_0(N, s1) + ad_lo_as1.gamma_gq_0(N), + +ad_lo_aem1.gamma_ns_0(N, s1) + ad_lo_aem1.gamma_phq_0(N), + 0, + ) def test_gluon_momentum_conservation(): @@ -27,24 +32,41 @@ def test_gluon_momentum_conservation(): N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) np.testing.assert_almost_equal( - ad_lo.gamma_qg_0(N, NF) + ad_lo.gamma_gg_0(N, s1, NF), 0 + ad_lo_as1.gamma_qg_0(N, NF) + ad_lo_as1.gamma_gg_0(N, s1, NF), 0 ) def test_gamma_qg_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) - np.testing.assert_almost_equal(ad_lo.gamma_qg_0(N, NF), res) + np.testing.assert_almost_equal(ad_lo_as1.gamma_qg_0(N, NF), res) def test_gamma_gq_0(): N = complex(0.0, 1.0) res = complex(4.0, -4.0) / 3.0 - np.testing.assert_almost_equal(ad_lo.gamma_gq_0(N), res) + np.testing.assert_almost_equal(ad_lo_as1.gamma_gq_0(N), res) def test_gamma_gg_0(): N = complex(0.0, 1.0) s1 = harmonics.harmonic_S1(N) res = complex(5.195725159621, 10.52008856962) - np.testing.assert_almost_equal(ad_lo.gamma_gg_0(N, s1, NF), res) + np.testing.assert_almost_equal(ad_lo_as1.gamma_gg_0(N, s1, NF), res) + + +def test_gamma_phq_0(): + N = complex(0.0, 1.0) + res = complex(4.0, -4.0) / 3.0 / 4 * 3 + np.testing.assert_almost_equal(ad_lo_aem1.gamma_phq_0(N), res) + + +def test_gamma_qph_0(): + N = complex(1.0, 0.0) + res = complex(-20.0 / 3.0, 0.0) * 3 / 0.5 + np.testing.assert_almost_equal(ad_lo_aem1.gamma_qph_0(N, NF), res) + + +def test_gamma_phph_0(): + res = complex(-4.0 / 3, 0.0) + np.testing.assert_almost_equal(ad_lo_aem1.gamma_phph_0(), res) diff --git a/tests/test_kernels_ns.py b/tests/test_kernels_ns.py index 3f2a7dd9e..f97d8e7eb 100644 --- a/tests/test_kernels_ns.py +++ b/tests/test_kernels_ns.py @@ -2,6 +2,7 @@ import numpy as np import pytest +from eko import anomalous_dimensions as ad from eko import beta from eko.kernels import non_singlet as ns @@ -114,6 +115,8 @@ def test_ode_nlo(): def test_error(): with pytest.raises(NotImplementedError): ns.dispatcher(3, "iterate-exact", np.random.rand(3) + 0j, 0.2, 0.1, 3, 10) + with pytest.raises(NotImplementedError): + ad.gamma_ns(1, 10202, 1, 3) def test_gamma_usage(): From b88ab720122f964408823306d0911603cd13e8cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 4 Mar 2022 12:37:04 +0100 Subject: [PATCH 14/31] Add the description of mode0 and mode1 --- src/eko/anomalous_dimensions/aem1.py | 2 +- src/eko/evolution_operator/__init__.py | 6 ++++-- src/eko/matching_conditions/operator_matrix_element.py | 6 ++++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py index dd5d26636..240c419df 100644 --- a/src/eko/anomalous_dimensions/aem1.py +++ b/src/eko/anomalous_dimensions/aem1.py @@ -27,7 +27,7 @@ def gamma_phq_0(N): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qph_0(N, nf: int): +def gamma_qph_0(N, nf): """ Computes the leading-order quark-photon anomalous dimension diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 20adf6f93..12a2e5796 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -125,8 +125,10 @@ def quad_ker( perturbation order method : str method - mode : str - sector element + mode0: int + pid for first sector element + mode1 : int + pid for second sector element is_log : boolean is a logarithmic interpolation logx : float diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index c7fb1340a..7623580d0 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -158,8 +158,10 @@ def quad_ker( quad argument order : int perturbation order - mode : str - element in the singlet sector + mode0 : int + pid for first element in the singlet sector + mode1 : int + pid for second element in the singlet sector is_log : boolean logarithmic interpolation logx : float From 72914498b94bf057ac3b1f48067c1655ebbc950d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 4 Mar 2022 15:37:08 +0100 Subject: [PATCH 15/31] Change non_singlet_pids_map into non_singlet_labels --- src/eko/evolution_operator/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 12a2e5796..4f047ea35 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -245,11 +245,11 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append((br.non_singlet_pids_map["ns+"], 0)) + labels.append(br.non_singlet_labels[1]) if order >= 1: # - becomes different starting from NLO - labels.append((br.non_singlet_pids_map["ns-"], 0)) + labels.append(br.non_singlet_labels[0]) if order >= 2: # v also becomes different starting from NNLO - labels.append((br.non_singlet_pids_map["nsV"], 0)) + labels.append(br.non_singlet_labels[2]) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") From 4cce3251222aaf9be3bff43adf83a27191822861 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 4 Mar 2022 16:26:49 +0100 Subject: [PATCH 16/31] Change nlo and nnlo into as2 and as3 in anomalous_dimensions --- src/eko/anomalous_dimensions/__init__.py | 34 +++++++++---------- .../anomalous_dimensions/{nlo.py => as2.py} | 0 .../anomalous_dimensions/{nnlo.py => as3.py} | 0 tests/benchmark_ad.py | 2 +- tests/test_ad_nlo.py | 2 +- tests/test_ad_nnlo.py | 2 +- 6 files changed, 20 insertions(+), 20 deletions(-) rename src/eko/anomalous_dimensions/{nlo.py => as2.py} (100%) rename src/eko/anomalous_dimensions/{nnlo.py => as3.py} (100%) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 78dd48fa6..7993bdd9c 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -21,7 +21,7 @@ import numpy as np from .. import basis_rotation as br -from . import as1, harmonics, nlo, nnlo +from . import as1, as2, as3, harmonics @nb.njit("Tuple((c16[:,:],c16,c16,c16[:,:],c16[:,:]))(c16[:,:])", cache=True) @@ -54,8 +54,8 @@ def exp_singlet(gamma_S): See Also -------- eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` """ # compute eigenvalues det = np.sqrt( @@ -96,11 +96,11 @@ def gamma_ns(order, mode, n, nf): See Also -------- eko.anomalous_dimensions.as1.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)` - eko.anomalous_dimensions.nlo.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsm_2 : :math:`\gamma_{ns,-}^{(2)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsv_2 : :math:`\gamma_{ns,v}^{(2)}(N)` + eko.anomalous_dimensions.as2.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)` + eko.anomalous_dimensions.as2.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsm_2 : :math:`\gamma_{ns,-}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsv_2 : :math:`\gamma_{ns,v}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) @@ -111,10 +111,10 @@ def gamma_ns(order, mode, n, nf): if order >= 1: # TODO: pass the necessary harmonics to nlo gammas if mode == 10101: - gamma_ns_1 = nlo.gamma_nsp_1(n, nf) + gamma_ns_1 = as2.gamma_nsp_1(n, nf) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here elif mode in [10201, 10200]: - gamma_ns_1 = nlo.gamma_nsm_1(n, nf) + gamma_ns_1 = as2.gamma_nsm_1(n, nf) else: raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 @@ -123,11 +123,11 @@ def gamma_ns(order, mode, n, nf): sx = np.append(sx, harmonics.harmonic_S2(n)) sx = np.append(sx, harmonics.harmonic_S3(n)) if mode == 10101: - gamma_ns_2 = -nnlo.gamma_nsp_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsp_2(n, nf, sx) elif mode == 10201: - gamma_ns_2 = -nnlo.gamma_nsm_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsm_2(n, nf, sx) elif mode == 10200: - gamma_ns_2 = -nnlo.gamma_nsv_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsv_2(n, nf, sx) gamma_ns[2] = gamma_ns_2 return gamma_ns @@ -154,8 +154,8 @@ def gamma_singlet(order, n, nf): See Also -------- eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) @@ -166,8 +166,8 @@ def gamma_singlet(order, n, nf): gamma_singlet = np.zeros((order + 1, 2, 2), np.complex_) gamma_singlet[0] = as1.gamma_singlet_0(n, sx[0], nf) if order >= 1: - gamma_singlet[1] = nlo.gamma_singlet_1(n, nf) + gamma_singlet[1] = as2.gamma_singlet_1(n, nf) if order == 2: sx = np.append(sx, harmonics.harmonic_S4(n)) - gamma_singlet[2] = -nnlo.gamma_singlet_2(n, nf, sx) + gamma_singlet[2] = -as3.gamma_singlet_2(n, nf, sx) return gamma_singlet diff --git a/src/eko/anomalous_dimensions/nlo.py b/src/eko/anomalous_dimensions/as2.py similarity index 100% rename from src/eko/anomalous_dimensions/nlo.py rename to src/eko/anomalous_dimensions/as2.py diff --git a/src/eko/anomalous_dimensions/nnlo.py b/src/eko/anomalous_dimensions/as3.py similarity index 100% rename from src/eko/anomalous_dimensions/nnlo.py rename to src/eko/anomalous_dimensions/as3.py diff --git a/tests/benchmark_ad.py b/tests/benchmark_ad.py index 465b75b3a..e37109cb0 100644 --- a/tests/benchmark_ad.py +++ b/tests/benchmark_ad.py @@ -2,8 +2,8 @@ """Benchmark the NLO anomalous dimensions against PEGASUS""" import numpy as np +import eko.anomalous_dimensions.as2 as ad_nlo import eko.anomalous_dimensions.harmonics as h -import eko.anomalous_dimensions.nlo as ad_nlo from eko.constants import CA, CF, TR diff --git a/tests/test_ad_nlo.py b/tests/test_ad_nlo.py index 82033d3a7..d3c69b6fb 100644 --- a/tests/test_ad_nlo.py +++ b/tests/test_ad_nlo.py @@ -2,8 +2,8 @@ # Test NLO anomalous dims import numpy as np +import eko.anomalous_dimensions.as2 as ad_nlo import eko.anomalous_dimensions.harmonics as h -import eko.anomalous_dimensions.nlo as ad_nlo from eko.constants import CA, CF NF = 5 diff --git a/tests/test_ad_nnlo.py b/tests/test_ad_nnlo.py index 923d4a3c7..c37bed82a 100644 --- a/tests/test_ad_nnlo.py +++ b/tests/test_ad_nnlo.py @@ -2,7 +2,7 @@ # Test NNLO anomalous dims import numpy as np -import eko.anomalous_dimensions.nnlo as ad_nnlo +import eko.anomalous_dimensions.as3 as ad_nnlo from eko.anomalous_dimensions import harmonics NF = 5 From 3d24bc21a5be2d456b715a1131e5fa95f85bae12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Wed, 9 Mar 2022 15:18:35 +0100 Subject: [PATCH 17/31] Apply some pylint suggestions --- src/eko/kernels/singlet.py | 10 +++++----- src/eko/strong_coupling.py | 2 +- src/ekomark/benchmark/runner.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index aa5d19f33..40ac12e0b 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -518,22 +518,22 @@ def dispatcher( # pylint: disable=too-many-return-statements # Common method for NLO and NNLO if method in ["iterate-exact", "iterate-expanded"]: return eko_iterate(gamma_singlet, a1, a0, nf, order, ev_op_iterations) - elif method == "perturbative-exact": + if method == "perturbative-exact": return eko_perturbative( gamma_singlet, a1, a0, nf, order, ev_op_iterations, ev_op_max_order, True ) - elif method == "perturbative-expanded": + if method == "perturbative-expanded": return eko_perturbative( gamma_singlet, a1, a0, nf, order, ev_op_iterations, ev_op_max_order, False ) - elif method in ["truncated", "ordered-truncated"]: + if method in ["truncated", "ordered-truncated"]: return eko_truncated(gamma_singlet, a1, a0, nf, order, ev_op_iterations) # These methods are scattered for nlo and nnlo - elif method == "decompose-exact": + if method == "decompose-exact": if order == 1: return nlo_decompose_exact(gamma_singlet, a1, a0, nf) return nnlo_decompose_exact(gamma_singlet, a1, a0, nf) - elif method == "decompose-expanded": + if method == "decompose-expanded": if order == 1: return nlo_decompose_expanded(gamma_singlet, a1, a0, nf) return nnlo_decompose_expanded(gamma_singlet, a1, a0, nf) diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 0642ecdef..455a48b5d 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -23,7 +23,7 @@ def strong_coupling_mod_ev(mod_ev): """Map ModEv key to the available strong coupling evolution methods""" if mod_ev in ["EXA", "iterate-exact", "decompose-exact", "perturbative-exact"]: return "exact" - elif mod_ev in [ + if mod_ev in [ "TRN", "truncated", "ordered-truncated", diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index 40ccab8a7..1eabf6fe4 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -129,7 +129,7 @@ def run_external(self, theory, ocard, pdf): ocard, rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "lhapdf": + if self.external.lower() == "lhapdf": from .external import lhapdf_utils # here theory card is not needed @@ -139,7 +139,7 @@ def run_external(self, theory, ocard, pdf): self.skip_pdfs(theory), rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "pegasus": + if self.external.lower() == "pegasus": from .external import pegasus_utils return pegasus_utils.compute_pegasus_data( @@ -149,7 +149,7 @@ def run_external(self, theory, ocard, pdf): rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "apfel": + if self.external.lower() == "apfel": from .external import apfel_utils return apfel_utils.compute_apfel_data( From b28395af5ee5378bff0bc816b588c68b14c8e1fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Wed, 9 Mar 2022 15:25:01 +0100 Subject: [PATCH 18/31] Fix quad_ker signature --- 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 078aadec8..20018066d 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -122,7 +122,7 @@ def integrand( return self.path.prefactor * pj * self.path.jac -@nb.njit("f8(f8,u1,u2,u2,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1,u1)", cache=True) +@nb.njit("f8(f8,u1,u2,u2,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1,u1)", cache=True) def quad_ker( u, order, From 846e550b36b137692b4f2136f4451ed428abd85d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Wed, 9 Mar 2022 16:38:46 +0100 Subject: [PATCH 19/31] Fix documentation of gamma_ns in anomalous_dimensions/__init__.py --- src/eko/anomalous_dimensions/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 7993bdd9c..bf3305957 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -81,7 +81,7 @@ def gamma_ns(order, mode, n, nf): ---------- order : int perturbative order - mode : "-" | "+" | "V" + mode : 10201 | 10101 | 10200 sector identifier n : complex Mellin variable From 3552d01faa0c6d04ecd191ece6f9aeb83b5edd78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Wed, 9 Mar 2022 16:49:50 +0100 Subject: [PATCH 20/31] Add function update_colors in constants.py --- src/eko/constants.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/eko/constants.py b/src/eko/constants.py index 470168305..44b96da63 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -8,16 +8,23 @@ """the number of colors""" TR = float(1.0 / 2.0) -"""the normalization of fundamental generators""" +"""the normalization of fundamental generators - defaults to :math:`T_F = 1/2`""" CA = float(NC) -"""second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3""" +"""second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3`""" CF = float((NC * NC - 1.0) / (2.0 * NC)) -"""second Casimir constant in the fundamental representation - defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3""" +"""second Casimir constant in the fundamental representation - defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`""" eu2 = 4.0 / 9 """up quarks charge squared""" ed2 = 1.0 / 9 """down quarks charge squared""" + + +def update_colors(nc): + global NC, CA, CF + NC = int(nc) + CA = float(NC) + CF = float((NC * NC - 1.0) / (2.0 * NC)) From d7397f74d7fb414cf3c4d745fb8ccaa4fd67954f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Wed, 9 Mar 2022 17:50:09 +0100 Subject: [PATCH 21/31] Remove index 0 from QED LO splitting functions --- src/eko/anomalous_dimensions/aem1.py | 16 ++++++++-------- tests/eko/test_ad_lo.py | 8 ++++---- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py index 240c419df..22ce20914 100644 --- a/src/eko/anomalous_dimensions/aem1.py +++ b/src/eko/anomalous_dimensions/aem1.py @@ -6,7 +6,7 @@ @nb.njit("c16(c16)", cache=True) -def gamma_phq_0(N): +def gamma_phq(N): """ Computes the leading-order photon-quark anomalous dimension @@ -19,7 +19,7 @@ def gamma_phq_0(N): Returns ------- - gamma_phq_0 : complex + gamma_phq : complex Leading-order photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(0)}(N)` """ @@ -27,7 +27,7 @@ def gamma_phq_0(N): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qph_0(N, nf): +def gamma_qph(N, nf): """ Computes the leading-order quark-photon anomalous dimension @@ -44,14 +44,14 @@ def gamma_qph_0(N, nf): Returns ------- - gamma_qph_0 : complex + gamma_qph : complex Leading-order quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(0)}(N)` """ return as1.gamma_qg_0(N, nf) / constants.TR * constants.NC @nb.njit("c16()", cache=True) -def gamma_phph_0(): +def gamma_phph(): """ Computes the leading-order photon-photon anomalous dimension @@ -59,7 +59,7 @@ def gamma_phph_0(): Returns ------- - gamma_phph_0 : complex + gamma_phph : complex Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)` """ @@ -67,7 +67,7 @@ def gamma_phph_0(): @nb.njit("c16(c16,c16)", cache=True) -def gamma_ns_0(N, s1): +def gamma_ns(N, s1): """ Computes the leading-order non-singlet QED anomalous dimension. @@ -82,7 +82,7 @@ def gamma_ns_0(N, s1): Returns ------- - gamma_ns_0 : complex + gamma_ns : complex Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` """ return as1.gamma_ns_0(N, s1) / constants.CF diff --git a/tests/eko/test_ad_lo.py b/tests/eko/test_ad_lo.py index 98b60e9e8..2e794b421 100644 --- a/tests/eko/test_ad_lo.py +++ b/tests/eko/test_ad_lo.py @@ -22,7 +22,7 @@ def test_quark_momentum_conservation(): s1 = harmonics.harmonic_S1(N) np.testing.assert_almost_equal( ad_lo_as1.gamma_ns_0(N, s1) + ad_lo_as1.gamma_gq_0(N), - +ad_lo_aem1.gamma_ns_0(N, s1) + ad_lo_aem1.gamma_phq_0(N), + +ad_lo_aem1.gamma_ns(N, s1) + ad_lo_aem1.gamma_phq(N), 0, ) @@ -58,15 +58,15 @@ def test_gamma_gg_0(): def test_gamma_phq_0(): N = complex(0.0, 1.0) res = complex(4.0, -4.0) / 3.0 / 4 * 3 - np.testing.assert_almost_equal(ad_lo_aem1.gamma_phq_0(N), res) + np.testing.assert_almost_equal(ad_lo_aem1.gamma_phq(N), res) def test_gamma_qph_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) * 3 / 0.5 - np.testing.assert_almost_equal(ad_lo_aem1.gamma_qph_0(N, NF), res) + np.testing.assert_almost_equal(ad_lo_aem1.gamma_qph(N, NF), res) def test_gamma_phph_0(): res = complex(-4.0 / 3, 0.0) - np.testing.assert_almost_equal(ad_lo_aem1.gamma_phph_0(), res) + np.testing.assert_almost_equal(ad_lo_aem1.gamma_phph(), res) From f535c916c04c4de5a201d1ffafaf3dc3837902ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 10 Mar 2022 10:21:37 +0100 Subject: [PATCH 22/31] Remove indices 0,1,2 from lo,nlo,nnlo QCD AD --- benchmarks/eko/benchmark_ad.py | 8 ++-- src/eko/anomalous_dimensions/__init__.py | 42 +++++++++---------- src/eko/anomalous_dimensions/aem1.py | 6 +-- src/eko/anomalous_dimensions/as1.py | 35 ++++++++-------- src/eko/anomalous_dimensions/as2.py | 45 ++++++++++---------- src/eko/anomalous_dimensions/as3.py | 52 ++++++++++++------------ tests/eko/test_ad.py | 4 +- tests/eko/test_ad_lo.py | 27 ++++++------ tests/eko/test_ad_nlo.py | 20 ++++----- tests/eko/test_ad_nnlo.py | 10 ++--- 10 files changed, 125 insertions(+), 124 deletions(-) diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py index 66bbe2e3d..1087db0e0 100644 --- a/benchmarks/eko/benchmark_ad.py +++ b/benchmarks/eko/benchmark_ad.py @@ -3,7 +3,7 @@ import numpy as np import pytest -import eko.anomalous_dimensions.as2 as ad_nlo +import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h from eko.constants import CA, CF, TR @@ -136,8 +136,8 @@ def check_gamma_1_pegasus(N, NF): P1NSP = CF * ((CF - CA / 2.0) * PNPA + CA * PNSB + TR * NF * PNSC) P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC) - np.testing.assert_allclose(ad_nlo.gamma_nsp_1(N, NF), -P1NSP) - np.testing.assert_allclose(ad_nlo.gamma_nsm_1(N, NF), -P1NSM) + np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF), -P1NSP) + np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF), -P1NSM) NS = N * N NT = NS * N @@ -254,7 +254,7 @@ def check_gamma_1_pegasus(N, NF): P1Sgq = (CF * CF * PGQA + CF * CA * PGQB + TR * NF * CF * PGQC) * 4.0 P1Sgg = (CA * CA * PGGA + TR * NF * (CA * PGGB + CF * PGGC)) * 4.0 - gS1 = ad_nlo.gamma_singlet_1(N, NF) + gS1 = ad_as2.gamma_singlet(N, NF) np.testing.assert_allclose(gS1[0, 0], -P1Sqq) np.testing.assert_allclose(gS1[0, 1], -P1Sqg) np.testing.assert_allclose(gS1[1, 0], -P1Sgq) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index bf3305957..bdf7d67ef 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -53,9 +53,9 @@ def exp_singlet(gamma_S): See Also -------- - eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.as2.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.as3.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` """ # compute eigenvalues det = np.sqrt( @@ -95,26 +95,26 @@ def gamma_ns(order, mode, n, nf): See Also -------- - eko.anomalous_dimensions.as1.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)` - eko.anomalous_dimensions.as2.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)` - eko.anomalous_dimensions.as2.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)` - eko.anomalous_dimensions.as3.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)` - eko.anomalous_dimensions.as3.gamma_nsm_2 : :math:`\gamma_{ns,-}^{(2)}(N)` - eko.anomalous_dimensions.as3.gamma_nsv_2 : :math:`\gamma_{ns,v}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)` + eko.anomalous_dimensions.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) # now combine gamma_ns = np.zeros(order + 1, np.complex_) - gamma_ns[0] = as1.gamma_ns_0(n, sx[0]) + gamma_ns[0] = as1.gamma_ns(n, sx[0]) # NLO and beyond if order >= 1: # TODO: pass the necessary harmonics to nlo gammas if mode == 10101: - gamma_ns_1 = as2.gamma_nsp_1(n, nf) + gamma_ns_1 = as2.gamma_nsp(n, nf) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here elif mode in [10201, 10200]: - gamma_ns_1 = as2.gamma_nsm_1(n, nf) + gamma_ns_1 = as2.gamma_nsm(n, nf) else: raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 @@ -123,11 +123,11 @@ def gamma_ns(order, mode, n, nf): sx = np.append(sx, harmonics.harmonic_S2(n)) sx = np.append(sx, harmonics.harmonic_S3(n)) if mode == 10101: - gamma_ns_2 = -as3.gamma_nsp_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsp(n, nf, sx) elif mode == 10201: - gamma_ns_2 = -as3.gamma_nsm_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsm(n, nf, sx) elif mode == 10200: - gamma_ns_2 = -as3.gamma_nsv_2(n, nf, sx) + gamma_ns_2 = -as3.gamma_nsv(n, nf, sx) gamma_ns[2] = gamma_ns_2 return gamma_ns @@ -153,9 +153,9 @@ def gamma_singlet(order, n, nf): See Also -------- - eko.anomalous_dimensions.as1.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.as2.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.as3.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) @@ -164,10 +164,10 @@ def gamma_singlet(order, n, nf): sx = np.append(sx, harmonics.harmonic_S3(n)) gamma_singlet = np.zeros((order + 1, 2, 2), np.complex_) - gamma_singlet[0] = as1.gamma_singlet_0(n, sx[0], nf) + gamma_singlet[0] = as1.gamma_singlet(n, sx[0], nf) if order >= 1: - gamma_singlet[1] = as2.gamma_singlet_1(n, nf) + gamma_singlet[1] = as2.gamma_singlet(n, nf) if order == 2: sx = np.append(sx, harmonics.harmonic_S4(n)) - gamma_singlet[2] = -as3.gamma_singlet_2(n, nf, sx) + gamma_singlet[2] = -as3.gamma_singlet(n, nf, sx) return gamma_singlet diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py index 22ce20914..2151f384b 100644 --- a/src/eko/anomalous_dimensions/aem1.py +++ b/src/eko/anomalous_dimensions/aem1.py @@ -23,7 +23,7 @@ def gamma_phq(N): Leading-order photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(0)}(N)` """ - return as1.gamma_gq_0(N) / constants.CF + return as1.gamma_gq(N) / constants.CF @nb.njit("c16(c16,u1)", cache=True) @@ -47,7 +47,7 @@ def gamma_qph(N, nf): gamma_qph : complex Leading-order quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(0)}(N)` """ - return as1.gamma_qg_0(N, nf) / constants.TR * constants.NC + return as1.gamma_qg(N, nf) / constants.TR * constants.NC @nb.njit("c16()", cache=True) @@ -85,4 +85,4 @@ def gamma_ns(N, s1): gamma_ns : complex Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` """ - return as1.gamma_ns_0(N, s1) / constants.CF + return as1.gamma_ns(N, s1) / constants.CF diff --git a/src/eko/anomalous_dimensions/as1.py b/src/eko/anomalous_dimensions/as1.py index 847fb4ba5..2efc0cae5 100644 --- a/src/eko/anomalous_dimensions/as1.py +++ b/src/eko/anomalous_dimensions/as1.py @@ -8,7 +8,7 @@ @nb.njit("c16(c16,c16)", cache=True) -def gamma_ns_0(N, s1): +def gamma_ns(N, s1): """ Computes the leading-order non-singlet anomalous dimension. @@ -23,7 +23,7 @@ def gamma_ns_0(N, s1): Returns ------- - gamma_ns_0 : complex + 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)) @@ -32,7 +32,7 @@ def gamma_ns_0(N, s1): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_0(N, nf): +def gamma_qg(N, nf): """ Computes the leading-order quark-gluon anomalous dimension @@ -47,7 +47,7 @@ def gamma_qg_0(N, nf): Returns ------- - gamma_qg_0 : complex + gamma_qg : complex Leading-order quark-gluon anomalous dimension :math:`\\gamma_{qg}^{(0)}(N)` """ gamma = -(N**2 + N + 2.0) / (N * (N + 1.0) * (N + 2.0)) @@ -56,7 +56,7 @@ def gamma_qg_0(N, nf): @nb.njit("c16(c16)", cache=True) -def gamma_gq_0(N): +def gamma_gq(N): """ Computes the leading-order gluon-quark anomalous dimension @@ -69,7 +69,7 @@ def gamma_gq_0(N): Returns ------- - gamma_gq_0 : complex + gamma_gq : complex Leading-order gluon-quark anomalous dimension :math:`\\gamma_{gq}^{(0)}(N)` """ gamma = -(N**2 + N + 2.0) / (N * (N + 1.0) * (N - 1.0)) @@ -78,7 +78,7 @@ def gamma_gq_0(N): @nb.njit("c16(c16,c16,u1)", cache=True) -def gamma_gg_0(N, s1, nf): +def gamma_gg(N, s1, nf): """ Computes the leading-order gluon-gluon anomalous dimension @@ -95,7 +95,7 @@ def gamma_gg_0(N, s1, nf): Returns ------- - gamma_gg_0 : complex + 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) @@ -104,7 +104,7 @@ def gamma_gg_0(N, s1, nf): @nb.njit("c16[:,:](c16,c16,u1)", cache=True) -def gamma_singlet_0(N, s1, nf): +def gamma_singlet(N, s1, nf): r""" Computes the leading-order singlet anomalous dimension matrix @@ -130,14 +130,13 @@ def gamma_singlet_0(N, s1, nf): See Also -------- - gamma_ns_0 : :math:`\gamma_{qq}^{(0)}` - gamma_qg_0 : :math:`\gamma_{qg}^{(0)}` - gamma_gq_0 : :math:`\gamma_{gq}^{(0)}` - gamma_gg_0 : :math:`\gamma_{gg}^{(0)}` + gamma_ns : :math:`\gamma_{qq}^{(0)}` + gamma_qg : :math:`\gamma_{qg}^{(0)}` + gamma_gq : :math:`\gamma_{gq}^{(0)}` + gamma_gg : :math:`\gamma_{gg}^{(0)}` """ - gamma_qq = gamma_ns_0(N, s1) - gamma_qg = gamma_qg_0(N, nf) - gamma_gq = gamma_gq_0(N) - gamma_gg = gamma_gg_0(N, s1, nf) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_ns(N, s1) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, s1, nf)]], np.complex_ + ) return gamma_S_0 diff --git a/src/eko/anomalous_dimensions/as2.py b/src/eko/anomalous_dimensions/as2.py index 3466be105..7c477c591 100644 --- a/src/eko/anomalous_dimensions/as2.py +++ b/src/eko/anomalous_dimensions/as2.py @@ -15,7 +15,7 @@ @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsm_1(n, nf): +def gamma_nsm(n, nf): """ Computes the |NLO| valence-like non-singlet anomalous dimension. @@ -30,7 +30,7 @@ def gamma_nsm_1(n, nf): Returns ------- - gamma_nsm_1 : complex + gamma_nsm : complex |NLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(1)}(N)` """ @@ -58,7 +58,7 @@ def gamma_nsm_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsp_1(n, nf): +def gamma_nsp(n, nf): """ Computes the |NLO| singlet-like non-singlet anomalous dimension. @@ -73,7 +73,7 @@ def gamma_nsp_1(n, nf): Returns ------- - gamma_nsp_1 : complex + gamma_nsp : complex |NLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(1)}(N)` """ @@ -99,7 +99,7 @@ def gamma_nsp_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_ps_1(n, nf): +def gamma_ps(n, nf): """ Computes the |NLO| pure-singlet quark-quark anomalous dimension. @@ -114,7 +114,7 @@ def gamma_ps_1(n, nf): Returns ------- - gamma_ps_1 : complex + gamma_ps : complex |NLO| pure-singlet quark-quark anomalous dimension :math:`\\gamma_{ps}^{(1)}(N)` """ @@ -126,7 +126,7 @@ def gamma_ps_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_1(n, nf): +def gamma_qg(n, nf): """ Computes the |NLO| quark-gluon singlet anomalous dimension. @@ -141,7 +141,7 @@ def gamma_qg_1(n, nf): Returns ------- - gamma_qg_1 : complex + gamma_qg : complex |NLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(1)}(N)` """ @@ -159,7 +159,7 @@ def gamma_qg_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gq_1(n, nf): +def gamma_gq(n, nf): """ Computes the |NLO| gluon-quark singlet anomalous dimension. @@ -174,7 +174,7 @@ def gamma_gq_1(n, nf): Returns ------- - gamma_gq_1 : complex + gamma_gq : complex |NLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(1)}(N)` """ @@ -195,7 +195,7 @@ def gamma_gq_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gg_1(n, nf): +def gamma_gg(n, nf): """ Computes the |NLO| gluon-gluon singlet anomalous dimension. @@ -210,7 +210,7 @@ def gamma_gg_1(n, nf): Returns ------- - gamma_gg_1 : complex + gamma_gg : complex |NLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(1)}(N)` """ @@ -233,7 +233,7 @@ def gamma_gg_1(n, nf): @nb.njit("c16[:,:](c16,u1)", cache=True) -def gamma_singlet_1(N, nf): +def gamma_singlet(N, nf): r""" Computes the next-leading-order singlet anomalous dimension matrix @@ -257,15 +257,14 @@ def gamma_singlet_1(N, nf): See Also -------- - gamma_nsp_1 : :math:`\gamma_{qq}^{(1)}` - gamma_ps_1 : :math:`\gamma_{qq}^{(1)}` - gamma_qg_1 : :math:`\gamma_{qg}^{(1)}` - gamma_gq_1 : :math:`\gamma_{gq}^{(1)}` - gamma_gg_1 : :math:`\gamma_{gg}^{(1)}` + gamma_nsp : :math:`\gamma_{qq}^{(1)}` + gamma_ps : :math:`\gamma_{qq}^{(1)}` + gamma_qg : :math:`\gamma_{qg}^{(1)}` + gamma_gq : :math:`\gamma_{gq}^{(1)}` + gamma_gg : :math:`\gamma_{gg}^{(1)}` """ - gamma_qq = gamma_nsp_1(N, nf) + gamma_ps_1(N, nf) - gamma_qg = gamma_qg_1(N, nf) - gamma_gq = gamma_gq_1(N, nf) - gamma_gg = gamma_gg_1(N, nf) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_nsp(N, nf) + gamma_ps(N, nf) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N, nf), gamma_gg(N, nf)]], np.complex_ + ) return gamma_S_0 diff --git a/src/eko/anomalous_dimensions/as3.py b/src/eko/anomalous_dimensions/as3.py index 46a63eb5c..c0864d2ff 100644 --- a/src/eko/anomalous_dimensions/as3.py +++ b/src/eko/anomalous_dimensions/as3.py @@ -17,7 +17,7 @@ @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsm_2(n, nf, sx): +def gamma_nsm(n, nf, sx): """ Computes the |NNLO| valence-like non-singlet anomalous dimension. @@ -34,7 +34,7 @@ def gamma_nsm_2(n, nf, sx): Returns ------- - gamma_nsm_2 : complex + gamma_nsm : complex |NNLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(2)}(N)` """ @@ -94,7 +94,7 @@ def gamma_nsm_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsp_2(n, nf, sx): +def gamma_nsp(n, nf, sx): """ Computes the |NNLO| singlet-like non-singlet anomalous dimension. @@ -111,7 +111,7 @@ def gamma_nsp_2(n, nf, sx): Returns ------- - gamma_nsp_2 : complex + gamma_nsp : complex |NNLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(2)}(N)` """ @@ -171,7 +171,7 @@ def gamma_nsp_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsv_2(n, nf, sx): +def gamma_nsv(n, nf, sx): """ Computes the |NNLO| valence non-singlet anomalous dimension. @@ -188,7 +188,7 @@ def gamma_nsv_2(n, nf, sx): Returns ------- - gamma_nsv_2 : complex + gamma_nsv : complex |NNLO| valence non-singlet anomalous dimension :math:`\\gamma_{ns,v}^{(2)}(N)` """ @@ -221,12 +221,12 @@ def gamma_nsv_2(n, nf, sx): + 46.18 * E2 ) - result = gamma_nsm_2(n, nf, sx) + nf * ps2 + result = gamma_nsm(n, nf, sx) + nf * ps2 return result @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_ps_2(n, nf, sx): +def gamma_ps(n, nf, sx): """ Computes the |NNLO| pure-singlet quark-quark anomalous dimension. @@ -243,7 +243,7 @@ def gamma_ps_2(n, nf, sx): Returns ------- - gamma_ps_2 : complex + gamma_ps : complex |NNLO| pure-singlet quark-quark anomalous dimension :math:`\\gamma_{ps}^{(2)}(N)` """ @@ -298,7 +298,7 @@ def gamma_ps_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_qg_2(n, nf, sx): +def gamma_qg(n, nf, sx): """ Computes the |NNLO| quark-gluon singlet anomalous dimension. @@ -315,7 +315,7 @@ def gamma_qg_2(n, nf, sx): Returns ------- - gamma_qg_2 : complex + gamma_qg : complex |NNLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(2)}(N)` """ @@ -372,7 +372,7 @@ def gamma_qg_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gq_2(n, nf, sx): +def gamma_gq(n, nf, sx): """ Computes the |NNLO| gluon-quark singlet anomalous dimension. @@ -389,7 +389,7 @@ def gamma_gq_2(n, nf, sx): Returns ------- - gamma_gq_2 : complex + gamma_gq : complex |NNLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(2)}(N)` """ @@ -462,7 +462,7 @@ def gamma_gq_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gg_2(n, nf, sx): +def gamma_gg(n, nf, sx): """ Computes the |NNLO| gluon-gluon singlet anomalous dimension. @@ -479,7 +479,7 @@ def gamma_gg_2(n, nf, sx): Returns ------- - gamma_gg_2 : complex + gamma_gg : complex |NNLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(2)}(N)` """ @@ -550,7 +550,7 @@ def gamma_gg_2(n, nf, sx): @nb.njit("c16[:,:](c16,u1,c16[:])", cache=True) -def gamma_singlet_2(N, nf, sx): +def gamma_singlet(N, nf, sx): r""" Computes the |NNLO| singlet anomalous dimension matrix @@ -578,15 +578,15 @@ def gamma_singlet_2(N, nf, sx): See Also -------- - gamma_nsp_2 : :math:`\gamma_{qq}^{(2)}` - gamma_ps_2 : :math:`\gamma_{qq}^{(2)}` - gamma_qg_2 : :math:`\gamma_{qg}^{(2)}` - gamma_gq_2 : :math:`\gamma_{gq}^{(2)}` - gamma_gg_2 : :math:`\gamma_{gg}^{(2)}` + gamma_nsp : :math:`\gamma_{qq}^{(2)}` + gamma_ps : :math:`\gamma_{qq}^{(2)}` + gamma_qg : :math:`\gamma_{qg}^{(2)}` + gamma_gq : :math:`\gamma_{gq}^{(2)}` + gamma_gg : :math:`\gamma_{gg}^{(2)}` """ - gamma_qq = gamma_nsp_2(N, nf, sx) + gamma_ps_2(N, nf, sx) - gamma_qg = gamma_qg_2(N, nf, sx) - gamma_gq = gamma_gq_2(N, nf, sx) - gamma_gg = gamma_gg_2(N, nf, sx) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_nsp(N, nf, sx) + gamma_ps(N, nf, sx) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf, sx)], [gamma_gq(N, nf, sx), gamma_gg(N, nf, sx)]], + np.complex_, + ) return gamma_S_0 diff --git a/tests/eko/test_ad.py b/tests/eko/test_ad.py index 3255e1751..2627decd5 100644 --- a/tests/eko/test_ad.py +++ b/tests/eko/test_ad.py @@ -5,7 +5,7 @@ from eko import anomalous_dimensions as ad from eko import basis_rotation as br -from eko.anomalous_dimensions import as1 as ad_lo +from eko.anomalous_dimensions import as1 as ad_as1 from eko.anomalous_dimensions import harmonics NF = 5 @@ -14,7 +14,7 @@ def test_eigensystem_gamma_singlet_0_values(): n = 3 s1 = harmonics.harmonic_S1(n) - gamma_S_0 = ad_lo.gamma_singlet_0(3, s1, NF) + gamma_S_0 = ad_as1.gamma_singlet(3, s1, NF) res = ad.exp_singlet(gamma_S_0) lambda_p = complex(12.273612971466964, 0) lambda_m = complex(5.015275917421917, 0) diff --git a/tests/eko/test_ad_lo.py b/tests/eko/test_ad_lo.py index 2e794b421..4dfafe877 100644 --- a/tests/eko/test_ad_lo.py +++ b/tests/eko/test_ad_lo.py @@ -2,8 +2,8 @@ # Test LO splitting functions import numpy as np -import eko.anomalous_dimensions.aem1 as ad_lo_aem1 -import eko.anomalous_dimensions.as1 as ad_lo_as1 +import eko.anomalous_dimensions.aem1 as ad_aem1 +import eko.anomalous_dimensions.as1 as ad_as1 from eko.anomalous_dimensions import harmonics NF = 5 @@ -13,7 +13,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) s1 = harmonics.harmonic_S1(N) - np.testing.assert_almost_equal(ad_lo_as1.gamma_ns_0(N, s1), 0) + np.testing.assert_almost_equal(ad_as1.gamma_ns(N, s1), 0) def test_quark_momentum_conservation(): @@ -21,8 +21,11 @@ def test_quark_momentum_conservation(): N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) np.testing.assert_almost_equal( - ad_lo_as1.gamma_ns_0(N, s1) + ad_lo_as1.gamma_gq_0(N), - +ad_lo_aem1.gamma_ns(N, s1) + ad_lo_aem1.gamma_phq(N), + ad_as1.gamma_ns(N, s1) + ad_as1.gamma_gq(N), + 0, + ) + np.testing.assert_almost_equal( + ad_aem1.gamma_ns(N, s1) + ad_aem1.gamma_phq(N), 0, ) @@ -32,41 +35,41 @@ def test_gluon_momentum_conservation(): N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) np.testing.assert_almost_equal( - ad_lo_as1.gamma_qg_0(N, NF) + ad_lo_as1.gamma_gg_0(N, s1, NF), 0 + ad_as1.gamma_qg(N, NF) + ad_as1.gamma_gg(N, s1, NF), 0 ) def test_gamma_qg_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) - np.testing.assert_almost_equal(ad_lo_as1.gamma_qg_0(N, NF), res) + np.testing.assert_almost_equal(ad_as1.gamma_qg(N, NF), res) def test_gamma_gq_0(): N = complex(0.0, 1.0) res = complex(4.0, -4.0) / 3.0 - np.testing.assert_almost_equal(ad_lo_as1.gamma_gq_0(N), res) + np.testing.assert_almost_equal(ad_as1.gamma_gq(N), res) def test_gamma_gg_0(): N = complex(0.0, 1.0) s1 = harmonics.harmonic_S1(N) res = complex(5.195725159621, 10.52008856962) - np.testing.assert_almost_equal(ad_lo_as1.gamma_gg_0(N, s1, NF), res) + np.testing.assert_almost_equal(ad_as1.gamma_gg(N, s1, NF), res) def test_gamma_phq_0(): N = complex(0.0, 1.0) res = complex(4.0, -4.0) / 3.0 / 4 * 3 - np.testing.assert_almost_equal(ad_lo_aem1.gamma_phq(N), res) + np.testing.assert_almost_equal(ad_aem1.gamma_phq(N), res) def test_gamma_qph_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) * 3 / 0.5 - np.testing.assert_almost_equal(ad_lo_aem1.gamma_qph(N, NF), res) + np.testing.assert_almost_equal(ad_aem1.gamma_qph(N, NF), res) def test_gamma_phph_0(): res = complex(-4.0 / 3, 0.0) - np.testing.assert_almost_equal(ad_lo_aem1.gamma_phph(), res) + np.testing.assert_almost_equal(ad_aem1.gamma_phph(), res) diff --git a/tests/eko/test_ad_nlo.py b/tests/eko/test_ad_nlo.py index d3c69b6fb..1f7e45fa1 100644 --- a/tests/eko/test_ad_nlo.py +++ b/tests/eko/test_ad_nlo.py @@ -2,7 +2,7 @@ # Test NLO anomalous dims import numpy as np -import eko.anomalous_dimensions.as2 as ad_nlo +import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h from eko.constants import CA, CF @@ -11,9 +11,9 @@ def test_gamma_1(): # number conservation - np.testing.assert_allclose(ad_nlo.gamma_nsm_1(1, NF), 0.0, atol=2e-6) + np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF), 0.0, atol=2e-6) - gS1 = ad_nlo.gamma_singlet_1(2, NF) + gS1 = ad_as2.gamma_singlet(2, NF) # gluon momentum conservation # the CA*NF term seems to be tough to compute, so raise the constraint ... np.testing.assert_allclose(gS1[0, 1] + gS1[1, 1], 0.0, atol=4e-5) @@ -25,10 +25,10 @@ def test_gamma_1(): # reference values are obtained from MMa # Non singlet sector np.testing.assert_allclose( - ad_nlo.gamma_nsp_1(2, NF), (-112.0 * CF + 376.0 * CA - 64.0 * NF) * CF / 27.0 + ad_as2.gamma_nsp(2, NF), (-112.0 * CF + 376.0 * CA - 64.0 * NF) * CF / 27.0 ) # singlet sector - np.testing.assert_allclose(ad_nlo.gamma_ps_1(2, NF), -40.0 * CF * NF / 27.0) + np.testing.assert_allclose(ad_as2.gamma_ps(2, NF), -40.0 * CF * NF / 27.0) np.testing.assert_allclose(gS1[0, 1], (-74.0 * CF - 35.0 * CA) * NF / 27.0) # qg np.testing.assert_allclose( gS1[1, 0], (112.0 * CF - 376.0 * CA + 104.0 * NF) * CF / 27.0 @@ -36,7 +36,7 @@ def test_gamma_1(): # add additional point at (analytical) continuation point np.testing.assert_allclose( - ad_nlo.gamma_nsm_1(2, NF), + ad_as2.gamma_nsm(2, NF), ( (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.zeta3) * CF + (373.0 / 9.0 - 34.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA @@ -45,7 +45,7 @@ def test_gamma_1(): * CF, ) np.testing.assert_allclose( - ad_nlo.gamma_nsp_1(3, NF), + ad_as2.gamma_nsp(3, NF), ( (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.zeta3) * CF + (459.0 / 8.0 - 43.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA @@ -53,8 +53,8 @@ def test_gamma_1(): ) * CF, ) - np.testing.assert_allclose(ad_nlo.gamma_ps_1(3, NF), -1391.0 * CF * NF / 5400.0) - gS1 = ad_nlo.gamma_singlet_1(3, NF) + np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * CF * NF / 5400.0) + gS1 = ad_as2.gamma_singlet(3, NF) np.testing.assert_allclose( gS1[1, 0], ( @@ -73,7 +73,7 @@ def test_gamma_1(): ), rtol=6e-7, ) # gg - gS1 = ad_nlo.gamma_singlet_1(4, NF) + gS1 = ad_as2.gamma_singlet(4, NF) np.testing.assert_allclose( gS1[0, 1], (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * NF ) # qg diff --git a/tests/eko/test_ad_nnlo.py b/tests/eko/test_ad_nnlo.py index c37bed82a..58381478a 100644 --- a/tests/eko/test_ad_nnlo.py +++ b/tests/eko/test_ad_nnlo.py @@ -2,7 +2,7 @@ # Test NNLO anomalous dims import numpy as np -import eko.anomalous_dimensions.as3 as ad_nnlo +import eko.anomalous_dimensions.as3 as ad_as3 from eko.anomalous_dimensions import harmonics NF = 5 @@ -26,13 +26,13 @@ def test_gamma_2(): # number conservation - each is 0 on its own, see :cite:`Moch:2004pa` N = 1 sx = get_sx(N) - np.testing.assert_allclose(ad_nnlo.gamma_nsv_2(N, NF, sx), 0.000960586, rtol=3e-7) - np.testing.assert_allclose(ad_nnlo.gamma_nsm_2(N, NF, sx), -0.000594225, rtol=6e-7) + np.testing.assert_allclose(ad_as3.gamma_nsv(N, NF, sx), 0.000960586, rtol=3e-7) + np.testing.assert_allclose(ad_as3.gamma_nsm(N, NF, sx), -0.000594225, rtol=6e-7) # get singlet sector N = 2 sx = get_sx(N) - gS2 = ad_nnlo.gamma_singlet_2(N, NF, sx) + gS2 = ad_as3.gamma_singlet(N, NF, sx) # gluon momentum conservation np.testing.assert_allclose(gS2[0, 1] + gS2[1, 1], 0.00388726, rtol=2e-6) @@ -42,4 +42,4 @@ def test_gamma_2(): assert gS2.shape == (2, 2) # test nsv_2 equal to referece value - np.testing.assert_allclose(ad_nnlo.gamma_nsv_2(N, NF, sx), -188.325593, rtol=3e-7) + np.testing.assert_allclose(ad_as3.gamma_nsv(N, NF, sx), -188.325593, rtol=3e-7) From 5c5286639cb59bbde2c92ed7e2611f19a4f8a42f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 10 Mar 2022 11:47:15 +0100 Subject: [PATCH 23/31] Test function update_colors --- tests/eko/test_ad_nlo.py | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/tests/eko/test_ad_nlo.py b/tests/eko/test_ad_nlo.py index 1f7e45fa1..e098c113e 100644 --- a/tests/eko/test_ad_nlo.py +++ b/tests/eko/test_ad_nlo.py @@ -4,7 +4,7 @@ import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h -from eko.constants import CA, CF +from eko.constants import CA, CF, update_colors NF = 5 @@ -77,3 +77,33 @@ def test_gamma_1(): np.testing.assert_allclose( gS1[0, 1], (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * NF ) # qg + + update_colors(4) + from eko.constants import CA as ca + from eko.constants import CF as cf + + np.testing.assert_allclose(ca, 4.0) + gS1 = ad_as2.gamma_singlet(3, NF) + np.testing.assert_allclose( + gS1[1, 0], + ( + 973.0 / 432.0 * cf + + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * ca + + 61.0 / 54.0 * NF + ) + * cf, + ) # gq + np.testing.assert_allclose( + gS1[1, 1], + ( + (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) * ca**2 + - 967.0 / 270.0 * ca * NF + + 541.0 / 216.0 * cf * NF + ), + rtol=6e-7, + ) # gg + gS1 = ad_as2.gamma_singlet(4, NF) + np.testing.assert_allclose( + gS1[0, 1], (-56317.0 / 18000.0 * cf + 16387.0 / 9000.0 * ca) * NF + ) # qg + update_colors(3) From e2aebc0a98fcb61296decc7fd8f181ae78423a99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 10 Mar 2022 12:27:34 +0100 Subject: [PATCH 24/31] Add doc string in update_colors --- src/eko/constants.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/eko/constants.py b/src/eko/constants.py index 44b96da63..1305abc92 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -24,6 +24,14 @@ def update_colors(nc): + """ + Updates the number of colors to :math:`NC = nc` and the Casimirs for a generic value of :math:`NC` + + Parameters + ---------- + nc : int + Number of colors + """ global NC, CA, CF NC = int(nc) CA = float(NC) From 50ac8584feefbcc48f4453a4e1e5175871638979 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 10:43:11 +0100 Subject: [PATCH 25/31] Fix test_ad_nlo.py --- src/eko/anomalous_dimensions/__init__.py | 4 +- src/eko/evolution_operator/__init__.py | 4 +- tests/eko/test_ad_nlo.py | 65 +++++++++++++----------- 3 files changed, 39 insertions(+), 34 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index bdf7d67ef..cde46ad76 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -116,7 +116,9 @@ def gamma_ns(order, mode, n, nf): elif mode in [10201, 10200]: gamma_ns_1 = as2.gamma_nsm(n, nf) else: - raise NotImplementedError("Non-singlet sector is not implemented") + raise NotImplementedError( + "Non-singlet sector " + str(mode) + " is not implemented" + ) gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order >= 2: diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 20018066d..3ed1db0a5 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -42,9 +42,9 @@ def select_singlet_element(ker, mode0, mode1): ker : numpy.ndarray singlet integration kernel mode0 : int - first sector element + id for first sector element mode1 : int - second sector element + id for second sector element Returns ------- ker : complex diff --git a/tests/eko/test_ad_nlo.py b/tests/eko/test_ad_nlo.py index e098c113e..23a541fae 100644 --- a/tests/eko/test_ad_nlo.py +++ b/tests/eko/test_ad_nlo.py @@ -4,7 +4,7 @@ import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h -from eko.constants import CA, CF, update_colors +from eko import constants as const NF = 5 @@ -25,85 +25,88 @@ def test_gamma_1(): # reference values are obtained from MMa # Non singlet sector np.testing.assert_allclose( - ad_as2.gamma_nsp(2, NF), (-112.0 * CF + 376.0 * CA - 64.0 * NF) * CF / 27.0 + ad_as2.gamma_nsp(2, NF), + (-112.0 * const.CF + 376.0 * const.CA - 64.0 * NF) * const.CF / 27.0, ) # singlet sector - np.testing.assert_allclose(ad_as2.gamma_ps(2, NF), -40.0 * CF * NF / 27.0) - np.testing.assert_allclose(gS1[0, 1], (-74.0 * CF - 35.0 * CA) * NF / 27.0) # qg + np.testing.assert_allclose(ad_as2.gamma_ps(2, NF), -40.0 * const.CF * NF / 27.0) np.testing.assert_allclose( - gS1[1, 0], (112.0 * CF - 376.0 * CA + 104.0 * NF) * CF / 27.0 + gS1[0, 1], (-74.0 * const.CF - 35.0 * const.CA) * NF / 27.0 + ) # qg + np.testing.assert_allclose( + gS1[1, 0], (112.0 * const.CF - 376.0 * const.CA + 104.0 * NF) * const.CF / 27.0 ) # gq # add additional point at (analytical) continuation point np.testing.assert_allclose( ad_as2.gamma_nsm(2, NF), ( - (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.zeta3) * CF - + (373.0 / 9.0 - 34.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA + (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.zeta3) * const.CF + + (373.0 / 9.0 - 34.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * const.CA - 64.0 * NF / 27.0 ) - * CF, + * const.CF, ) np.testing.assert_allclose( ad_as2.gamma_nsp(3, NF), ( - (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.zeta3) * CF - + (459.0 / 8.0 - 43.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA + (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.zeta3) * const.CF + + (459.0 / 8.0 - 43.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * const.CA - 415.0 * NF / 108.0 ) - * CF, + * const.CF, ) - np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * CF * NF / 5400.0) + np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * const.CF * NF / 5400.0) gS1 = ad_as2.gamma_singlet(3, NF) np.testing.assert_allclose( gS1[1, 0], ( - 973.0 / 432.0 * CF - + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * CA + 973.0 / 432.0 * const.CF + + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * const.CA + 61.0 / 54.0 * NF ) - * CF, + * const.CF, ) # gq np.testing.assert_allclose( gS1[1, 1], ( - (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) * CA**2 - - 967.0 / 270.0 * CA * NF - + 541.0 / 216.0 * CF * NF + (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) + * const.CA**2 + - 967.0 / 270.0 * const.CA * NF + + 541.0 / 216.0 * const.CF * NF ), rtol=6e-7, ) # gg gS1 = ad_as2.gamma_singlet(4, NF) np.testing.assert_allclose( - gS1[0, 1], (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * NF + gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF ) # qg - update_colors(4) - from eko.constants import CA as ca - from eko.constants import CF as cf + const.update_colors(4) - np.testing.assert_allclose(ca, 4.0) + np.testing.assert_allclose(const.CA, 4.0) gS1 = ad_as2.gamma_singlet(3, NF) np.testing.assert_allclose( gS1[1, 0], ( - 973.0 / 432.0 * cf - + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * ca + 973.0 / 432.0 * const.CF + + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * const.CA + 61.0 / 54.0 * NF ) - * cf, + * const.CF, ) # gq np.testing.assert_allclose( gS1[1, 1], ( - (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) * ca**2 - - 967.0 / 270.0 * ca * NF - + 541.0 / 216.0 * cf * NF + (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) + * const.CA**2 + - 967.0 / 270.0 * const.CA * NF + + 541.0 / 216.0 * const.CF * NF ), rtol=6e-7, ) # gg gS1 = ad_as2.gamma_singlet(4, NF) np.testing.assert_allclose( - gS1[0, 1], (-56317.0 / 18000.0 * cf + 16387.0 / 9000.0 * ca) * NF + gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF ) # qg - update_colors(3) + const.update_colors(3) From 26ae80cd889b503bcb9dfb6fb72e9c0770a2e17c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 10:56:39 +0100 Subject: [PATCH 26/31] Revert anomalous_dimensions/__init__.py --- src/eko/anomalous_dimensions/__init__.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index cde46ad76..bdf7d67ef 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -116,9 +116,7 @@ def gamma_ns(order, mode, n, nf): elif mode in [10201, 10200]: gamma_ns_1 = as2.gamma_nsm(n, nf) else: - raise NotImplementedError( - "Non-singlet sector " + str(mode) + " is not implemented" - ) + raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order >= 2: From 5283db71c8f98abd3c4798da7c00aa2c54d99f33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 12:26:47 +0100 Subject: [PATCH 27/31] Fix TR doc in constants.py --- src/eko/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eko/constants.py b/src/eko/constants.py index 1305abc92..5ac8e580a 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -8,7 +8,7 @@ """the number of colors""" TR = float(1.0 / 2.0) -"""the normalization of fundamental generators - defaults to :math:`T_F = 1/2`""" +"""the normalization of fundamental generators - defaults to :math:`T_R = 1/2`""" CA = float(NC) """second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3`""" From ad574e67f091112e2c6ea595ef2b8690cadf1eca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 14:56:14 +0100 Subject: [PATCH 28/31] Fix warning global-statement --- src/eko/constants.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/eko/constants.py b/src/eko/constants.py index 5ac8e580a..a9bb89027 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -32,7 +32,8 @@ def update_colors(nc): nc : int Number of colors """ - global NC, CA, CF + global NC, CA, CF # pylint: disable=global-statement + NC = int(nc) CA = float(NC) CF = float((NC * NC - 1.0) / (2.0 * NC)) From 5861c6aabf1548a813f476a191eda9d84c301af1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 15:03:33 +0100 Subject: [PATCH 29/31] Change elif with if in kernels/non_singlet.py --- src/eko/kernels/non_singlet.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py index a0750f9e2..0229e75df 100644 --- a/src/eko/kernels/non_singlet.py +++ b/src/eko/kernels/non_singlet.py @@ -338,30 +338,30 @@ def dispatcher( if order == 0: return lo_exact(gamma_ns, a1, a0, nf) # NLO - elif order == 1: + if order == 1: if method in [ "iterate-expanded", "decompose-expanded", "perturbative-expanded", ]: return nlo_expanded(gamma_ns, a1, a0, nf) - elif method == "truncated": + if method == "truncated": return nlo_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) - elif method == "ordered-truncated": + if method == "ordered-truncated": return nlo_ordered_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) # if method in ["iterate-exact", "decompose-exact", "perturbative-exact"]: return nlo_exact(gamma_ns, a1, a0, nf) # NNLO - elif order == 2: + if order == 2: if method in [ "iterate-expanded", "decompose-expanded", "perturbative-expanded", ]: return nnlo_expanded(gamma_ns, a1, a0, nf) - elif method == "truncated": + if method == "truncated": return nnlo_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) - elif method == "ordered-truncated": + if method == "ordered-truncated": return nnlo_ordered_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) # if method in ["iterate-exact", "decompose-exact", "perturbative-exact"]: return nnlo_exact(gamma_ns, a1, a0, nf) From f0c123bd0941435f810b1dded0b6a6c7e7da9971 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 17:19:57 +0100 Subject: [PATCH 30/31] Fix def of gamma_phph --- src/eko/anomalous_dimensions/aem1.py | 11 ++++++++--- tests/eko/test_ad_lo.py | 4 ++-- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py index 2151f384b..953c143c8 100644 --- a/src/eko/anomalous_dimensions/aem1.py +++ b/src/eko/anomalous_dimensions/aem1.py @@ -50,20 +50,25 @@ def gamma_qph(N, nf): return as1.gamma_qg(N, nf) / constants.TR * constants.NC -@nb.njit("c16()", cache=True) -def gamma_phph(): +@nb.njit("c16(u1)", cache=True) +def gamma_phph(nf): """ Computes the leading-order photon-photon anomalous dimension Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + Parameters + ---------- + nf : int + Number of active flavors + Returns ------- gamma_phph : complex Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)` """ - return -4.0 / 3 + return -2 / 3 * constants.NC * 2 * nf @nb.njit("c16(c16,c16)", cache=True) diff --git a/tests/eko/test_ad_lo.py b/tests/eko/test_ad_lo.py index 4dfafe877..c0d352676 100644 --- a/tests/eko/test_ad_lo.py +++ b/tests/eko/test_ad_lo.py @@ -71,5 +71,5 @@ def test_gamma_qph_0(): def test_gamma_phph_0(): - res = complex(-4.0 / 3, 0.0) - np.testing.assert_almost_equal(ad_aem1.gamma_phph(), res) + res = complex(-2.0 / 3 * 3 * 2 * NF, 0.0) + np.testing.assert_almost_equal(ad_aem1.gamma_phph(NF), res) From acd65aa9df2bdce017a013f7017533d705ac719c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 11 Mar 2022 17:42:56 +0100 Subject: [PATCH 31/31] Fix again def of gamma_phph --- src/eko/anomalous_dimensions/aem1.py | 2 +- tests/eko/test_ad_lo.py | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py index 953c143c8..de620655a 100644 --- a/src/eko/anomalous_dimensions/aem1.py +++ b/src/eko/anomalous_dimensions/aem1.py @@ -68,7 +68,7 @@ def gamma_phph(nf): Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)` """ - return -2 / 3 * constants.NC * 2 * nf + return 2 / 3 * constants.NC * 2 * nf @nb.njit("c16(c16,c16)", cache=True) diff --git a/tests/eko/test_ad_lo.py b/tests/eko/test_ad_lo.py index c0d352676..d08e79e32 100644 --- a/tests/eko/test_ad_lo.py +++ b/tests/eko/test_ad_lo.py @@ -39,6 +39,12 @@ def test_gluon_momentum_conservation(): ) +def test_photon_momentum_conservation(): + # gluon momentum + N = complex(2.0, 0.0) + np.testing.assert_almost_equal(ad_aem1.gamma_qph(N, NF) + ad_aem1.gamma_phph(NF), 0) + + def test_gamma_qg_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) @@ -71,5 +77,5 @@ def test_gamma_qph_0(): def test_gamma_phph_0(): - res = complex(-2.0 / 3 * 3 * 2 * NF, 0.0) + res = complex(2.0 / 3 * 3 * 2 * NF, 0.0) np.testing.assert_almost_equal(ad_aem1.gamma_phph(NF), res)