From 603ca1b188b96f4ef25d4f9ba6cf67649e9f71d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 29 Nov 2022 15:19:31 +0100 Subject: [PATCH 01/14] Init cache of different types of harmonic sums --- src/eko/harmonics/__init__.py | 93 +++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index 02cbe78c2..b284c9f39 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -265,3 +265,96 @@ def compute_qed_ns_cache(n, s1): g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) sx_qed_ns.append(g3Np2) return sx_qed_ns + + +@nb.njit(cache=True) +def compute_cache_S_n_half(n): + r"""Get the harmonics sums cache S_i(n/2). + + Parameters + ---------- + n : complex + Mellin moment + + Returns + ------- + list + harmonic sums cache. It contains: + + .. math :: + [S_1(n/2), S_2(n/2), S_3(n/2)] + """ + return [S1(n / 2.0), S2(n / 2.0), S3(n / 2.0)] + + +@nb.njit(cache=True) +def compute_cache_S_n_half_minus_1halph(n): + r"""Get the harmonics sums cache S_i((n-1)/2). + + Parameters + ---------- + n : complex + Mellin moment + + Returns + ------- + list + harmonic sums cache. It contains: + + .. math :: + [S_1((n-1)/2), S_2((n-1)/2), S_3((n-1)/2)] + """ + return [S1((n - 1) / 2), S2((n - 1) / 2), S3((n - 1) / 2)] + + +@nb.njit(cache=True) +def compute_cache_S_n_half_plus_1halph(n, sx): + r"""Get the harmonics sums cache S_i((n+1)/2). + + Parameters + ---------- + n : complex + Mellin moment + sx : list + harmonic sums cache S_i((n-1)/2) + + Returns + ------- + list + harmonic sums cache. It contains: + + .. math :: + [S_1((n+1)/2), S_2((n+1)/2), S_3((n+1)/2)] + """ + S1_nh = polygamma.recursive_harmonic_sum(sx[0], n, 1, 1) + S2_nh = polygamma.recursive_harmonic_sum(sx[1], n, 1, 2) + S3_nh = polygamma.recursive_harmonic_sum(sx[2], n, 1, 3) + return [S1_nh, S2_nh, S3_nh] + + +@nb.njit(cache=True) +def compute_cache_g3(n, s1, qed): + r"""Get the functionis g3(n) and g3(n+2). + + Parameters + ---------- + n : complex + Mellin moment + s1 : list + harmonic sum :math:`S_1(n)` + + + Returns + ------- + list + g3 functions cache. It contains: + + .. math :: + [g3(n), g3(n+2)] + """ + g_list = [g_functions.mellin_g3(n, s1)] + if qed: + S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) + g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) + g_list.append(g3Np2) + return g_list From aa4ffba71cc7133365dab8109066573ab3d934c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 2 Dec 2022 14:49:44 +0100 Subject: [PATCH 02/14] Change compute_new_sx_cache --- src/eko/anomalous_dimensions/__init__.py | 34 +++++++------ src/eko/anomalous_dimensions/aem2.py | 36 ++++++------- src/eko/anomalous_dimensions/as1aem1.py | 44 ++++++++-------- src/eko/anomalous_dimensions/as2.py | 64 +++++++++++++----------- src/eko/harmonics/__init__.py | 36 +++++++------ tests/eko/test_ad.py | 6 ++- tests/eko/test_ad_aem2.py | 4 +- tests/eko/test_ad_as1aem1.py | 4 +- tests/eko/test_ad_as2.py | 23 +++++---- 9 files changed, 134 insertions(+), 117 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index cd39b6652..d0fd75d46 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -152,11 +152,12 @@ def gamma_ns(order, mode, n, nf): gamma_ns[0] = as1.gamma_ns(n, sx[0]) # NLO and beyond if order[0] >= 2: + sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, False) if mode == 10101: - gamma_ns_1 = as2.gamma_nsp(n, nf, sx) + gamma_ns_1 = as2.gamma_nsp(n, nf, sx, sx_new) # 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(n, nf, sx) + gamma_ns_1 = as2.gamma_nsm(n, nf, sx, sx_new) else: raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 @@ -227,7 +228,8 @@ def gamma_singlet(order, n, nf): gamma_s = np.zeros((order[0], 2, 2), np.complex_) gamma_s[0] = as1.gamma_singlet(n, sx[0], nf) if order[0] >= 2: - gamma_s[1] = as2.gamma_singlet(n, nf, sx) + sx_new = harmonics.compute_new_sx_cache(n, sx[0], True, False) + gamma_s[1] = as2.gamma_singlet(n, nf, sx, sx_new) if order[0] >= 3: gamma_s[2] = as3.gamma_singlet(n, nf, sx) if order[0] >= 4: @@ -280,7 +282,7 @@ def gamma_ns_qed(order, mode, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0]) + sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) # now combine gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex_) if order[0] >= 1: @@ -288,18 +290,18 @@ def gamma_ns_qed(order, mode, n, nf): if order[1] >= 1: gamma_ns[0, 1] = choose_ns_ad_aem1(mode, n, sx) if order[0] >= 1 and order[1] >= 1: - gamma_ns[1, 1] = choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed) + gamma_ns[1, 1] = choose_ns_ad_as1aem1(mode, n, sx, sx_new) # NLO and beyond if order[0] >= 2: if mode in [10102, 10103]: - gamma_ns[2, 0] = as2.gamma_nsp(n, nf, sx) + gamma_ns[2, 0] = as2.gamma_nsp(n, nf, sx, sx_new) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here elif mode in [10202, 10203]: - gamma_ns[2, 0] = as2.gamma_nsm(n, nf, sx) + gamma_ns[2, 0] = as2.gamma_nsm(n, nf, sx, sx_new) else: raise NotImplementedError("Non-singlet sector is not implemented") if order[1] >= 2: - gamma_ns[0, 2] = choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed) + gamma_ns[0, 2] = choose_ns_ad_aem2(mode, n, nf, sx, sx_new) # NNLO and beyond if order[0] >= 3: if mode in [10102, 10103]: @@ -427,18 +429,18 @@ def gamma_singlet_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0]) + sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex_) if order[0] >= 1: gamma_s[1, 0] = as1.gamma_singlet_qed(n, sx[0], nf) if order[1] >= 1: gamma_s[0, 1] = aem1.gamma_singlet(n, nf, sx) if order[0] >= 1 and order[1] >= 1: - gamma_s[1, 1] = as1aem1.gamma_singlet(n, nf, sx, sx_ns_qed) + gamma_s[1, 1] = as1aem1.gamma_singlet(n, nf, sx, sx_new) if order[0] >= 2: - gamma_s[2, 0] = as2.gamma_singlet_qed(n, nf, sx) + gamma_s[2, 0] = as2.gamma_singlet_qed(n, nf, sx, sx_new) if order[1] >= 2: - gamma_s[0, 2] = aem2.gamma_singlet(n, nf, sx, sx_ns_qed) + gamma_s[0, 2] = aem2.gamma_singlet(n, nf, sx, sx_new) if order[0] >= 3: gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, sx) return gamma_s @@ -479,18 +481,18 @@ def gamma_valence_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0]) + sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex_) if order[0] >= 1: gamma_v[1, 0] = as1.gamma_valence_qed(n, sx[0]) if order[1] >= 1: gamma_v[0, 1] = aem1.gamma_valence(n, nf, sx) if order[0] >= 1 and order[1] >= 1: - gamma_v[1, 1] = as1aem1.gamma_valence(n, nf, sx, sx_ns_qed) + gamma_v[1, 1] = as1aem1.gamma_valence(n, nf, sx, sx_new) if order[0] >= 2: - gamma_v[2, 0] = as2.gamma_valence_qed(n, nf, sx) + gamma_v[2, 0] = as2.gamma_valence_qed(n, nf, sx, sx_new) if order[1] >= 2: - gamma_v[0, 2] = aem2.gamma_valence(n, nf, sx, sx_ns_qed) + gamma_v[0, 2] = aem2.gamma_valence(n, nf, sx, sx_new) if order[0] >= 3: gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, sx) return gamma_v diff --git a/src/eko/anomalous_dimensions/aem2.py b/src/eko/anomalous_dimensions/aem2.py index bb4cbc8f7..09a053829 100644 --- a/src/eko/anomalous_dimensions/aem2.py +++ b/src/eko/anomalous_dimensions/aem2.py @@ -151,7 +151,7 @@ def gamma_phd(N, nf, sx): @nb.njit(cache=True) -def gamma_nspu(N, nf, sx, sx_ns_qed): +def gamma_nspu(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` singlet-like non-singlet anomalous dimension for up quarks. Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u. @@ -184,13 +184,11 @@ def gamma_nspu(N, nf, sx, sx_ns_qed): - 80.0 / 9.0 * S1 + 16.0 / 3.0 * S2 ) * eSigma2 - return ( - constants.eu2 * as1aem1.gamma_nsp(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp - ) + return constants.eu2 * as1aem1.gamma_nsp(N, sx, sx_new) / constants.CF / 2.0 + tmp @nb.njit(cache=True) -def gamma_nspd(N, nf, sx, sx_ns_qed): +def gamma_nspd(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` singlet-like non-singlet anomalous dimension for down quarks. Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d. @@ -223,13 +221,11 @@ def gamma_nspd(N, nf, sx, sx_ns_qed): - 80.0 / 9.0 * S1 + 16.0 / 3.0 * S2 ) * eSigma2 - return ( - constants.ed2 * as1aem1.gamma_nsp(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp - ) + return constants.ed2 * as1aem1.gamma_nsp(N, sx, sx_new) / constants.CF / 2.0 + tmp @nb.njit(cache=True) -def gamma_nsmu(N, nf, sx, sx_ns_qed): +def gamma_nsmu(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` valence-like non-singlet anomalous dimension for up quarks. Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u. @@ -262,13 +258,11 @@ def gamma_nsmu(N, nf, sx, sx_ns_qed): - 80.0 / 9.0 * S1 + 16.0 / 3.0 * S2 ) * eSigma2 - return ( - constants.eu2 * as1aem1.gamma_nsm(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp - ) + return constants.eu2 * as1aem1.gamma_nsm(N, sx, sx_new) / constants.CF / 2.0 + tmp @nb.njit(cache=True) -def gamma_nsmd(N, nf, sx, sx_ns_qed): +def gamma_nsmd(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` valence-like non-singlet anomalous dimension for down quarks. Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d. @@ -301,9 +295,7 @@ def gamma_nsmd(N, nf, sx, sx_ns_qed): - 80.0 / 9.0 * S1 + 16.0 / 3.0 * S2 ) * eSigma2 - return ( - constants.ed2 * as1aem1.gamma_nsm(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp - ) + return constants.ed2 * as1aem1.gamma_nsm(N, sx, sx_new) / constants.CF / 2.0 + tmp @nb.njit(cache=True) @@ -336,7 +328,7 @@ def gamma_ps(N, nf): @nb.njit(cache=True) -def gamma_singlet(N, nf, sx, sx_ns_qed): +def gamma_singlet(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` singlet sector. Parameters @@ -363,8 +355,8 @@ def gamma_singlet(N, nf, sx, sx_ns_qed): gamma_ph_d = gamma_phd(N, nf, sx) gamma_u_ph = gamma_uph(N, nf, sx) gamma_d_ph = gamma_dph(N, nf, sx) - gamma_ns_p_u = gamma_nspu(N, nf, sx, sx_ns_qed) - gamma_ns_p_d = gamma_nspd(N, nf, sx, sx_ns_qed) + gamma_ns_p_u = gamma_nspu(N, nf, sx, sx_new) + gamma_ns_p_d = gamma_nspd(N, nf, sx, sx_new) gamma_pure_singlet = gamma_ps(N, nf) gamma_S_02 = np.array( [ @@ -408,7 +400,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed): @nb.njit(cache=True) -def gamma_valence(N, nf, sx, sx_ns_qed): +def gamma_valence(N, nf, sx, sx_new): r"""Compute the :math:`O(a_{em}^2)` valence sector. Parameters @@ -429,8 +421,8 @@ def gamma_valence(N, nf, sx, sx_ns_qed): nd = nf - nu vu = nu / nf vd = nd / nf - gamma_ns_m_u = gamma_nsmu(N, nf, sx, sx_ns_qed) - gamma_ns_m_d = gamma_nsmd(N, nf, sx, sx_ns_qed) + gamma_ns_m_u = gamma_nsmu(N, nf, sx, sx_new) + gamma_ns_m_d = gamma_nsmd(N, nf, sx, sx_new) gamma_V_02 = np.array( [ [ diff --git a/src/eko/anomalous_dimensions/as1aem1.py b/src/eko/anomalous_dimensions/as1aem1.py index 7b8eaa442..8b4b9a740 100644 --- a/src/eko/anomalous_dimensions/as1aem1.py +++ b/src/eko/anomalous_dimensions/as1aem1.py @@ -241,7 +241,7 @@ def gamma_gg(): @nb.njit(cache=True) -def gamma_nsp(N, sx, sx_ns_qed): +def gamma_nsp(N, sx, sx_new): r"""Compute the :math:`O(a_s^1a_{em}^1)` singlet-like non-singlet anomalous dimension. Implements sum of Eqs. (33-34) of :cite:`deFlorian:2015ujt`. @@ -263,14 +263,14 @@ def gamma_nsp(N, sx, sx_ns_qed): S1 = sx[0] S2 = sx[1] S3 = sx[2] - S1h = sx_ns_qed[0] - S2h = sx_ns_qed[1] - S3h = sx_ns_qed[2] - S1p1h = sx_ns_qed[3] - S2p1h = sx_ns_qed[4] - S3p1h = sx_ns_qed[5] - g3N = sx_ns_qed[6] - g3Np2 = sx_ns_qed[7] + S1h = sx_new[0] + S2h = sx_new[1] + S3h = sx_new[2] + g3N = sx_new[3] + g3Np2 = sx_new[7] + S1p1h = sx_new[8] + S2p1h = sx_new[9] + S3p1h = sx_new[10] result = ( +32.0 * zeta2 * S1h - 32.0 * zeta2 * S1p1h @@ -308,7 +308,7 @@ def gamma_nsp(N, sx, sx_ns_qed): @nb.njit(cache=True) -def gamma_nsm(N, sx, sx_ns_qed): +def gamma_nsm(N, sx, sx_new): r"""Compute the :math:`O(a_s^1a_{em}^1)` valence-like non-singlet anomalous dimension. Implements difference between Eqs. (33-34) of :cite:`deFlorian:2015ujt`. @@ -330,14 +330,14 @@ def gamma_nsm(N, sx, sx_ns_qed): S1 = sx[0] S2 = sx[1] S3 = sx[2] - S1h = sx_ns_qed[0] - S2h = sx_ns_qed[1] - S3h = sx_ns_qed[2] - S1p1h = sx_ns_qed[3] - S2p1h = sx_ns_qed[4] - S3p1h = sx_ns_qed[5] - g3N = sx_ns_qed[6] - g3Np2 = sx_ns_qed[7] + S1h = sx_new[0] + S2h = sx_new[1] + S3h = sx_new[2] + g3N = sx_new[3] + g3Np2 = sx_new[7] + S1p1h = sx_new[8] + S2p1h = sx_new[9] + S3p1h = sx_new[10] result = ( -32.0 * zeta2 * S1h - 8.0 / (N + N**2) * S2h @@ -371,7 +371,7 @@ def gamma_nsm(N, sx, sx_ns_qed): @nb.njit(cache=True) -def gamma_singlet(N, nf, sx, sx_ns_qed): +def gamma_singlet(N, nf, sx, sx_new): r"""Compute the :math:`O(a_s^1a_{em}^1)` singlet sector. Parameters @@ -397,7 +397,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed): gamma_ph_q = gamma_phq(N, sx) gamma_q_g = gamma_qg(N, nf, sx) gamma_q_ph = gamma_qph(N, nf, sx) - gamma_ns_p = gamma_nsp(N, sx, sx_ns_qed) + gamma_ns_p = gamma_nsp(N, sx, sx_new) gamma_S_11 = np.array( [ [ @@ -431,7 +431,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed): @nb.njit(cache=True) -def gamma_valence(N, nf, sx, sx_ns_qed): +def gamma_valence(N, nf, sx, sx_new): r"""Compute the :math:`O(a_s^1a_{em}^1)` valence sector. Parameters @@ -459,4 +459,4 @@ def gamma_valence(N, nf, sx, sx_ns_qed): ], np.complex_, ) - return gamma_V_11 * gamma_nsm(N, sx, sx_ns_qed) + return gamma_V_11 * gamma_nsm(N, sx, sx_new) diff --git a/src/eko/anomalous_dimensions/as2.py b/src/eko/anomalous_dimensions/as2.py index cfea2c9f0..7c0336e7b 100644 --- a/src/eko/anomalous_dimensions/as2.py +++ b/src/eko/anomalous_dimensions/as2.py @@ -13,7 +13,7 @@ @nb.njit(cache=True) -def gamma_nsm(n, nf, sx): +def gamma_nsm(n, nf, sx, sx_new): r""" Compute the |NLO| valence-like non-singlet anomalous dimension. @@ -38,10 +38,10 @@ def gamma_nsm(n, nf, sx): S2 = sx[1] # Here, Sp refers to S' ("s-prime") (german: "s-strich" or in Pegasus language: SSTR) # of :cite:`Gluck:1989ze` and NOT to the Spence function a.k.a. dilogarithm - Sp1m = harmonics.S1((n - 1) / 2) - Sp2m = harmonics.S2((n - 1) / 2) - Sp3m = harmonics.S3((n - 1) / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + g3n = sx_new[3] + Sp1m = sx_new[4] + Sp2m = sx_new[5] + Sp3m = sx_new[6] # fmt: off gqq1m_cfca = 16*g3n - (144 + n*(1 + n)*(156 + n*(340 + n*(655 + 51*n*(2 + n)))))/(18.*np.power(n,3)*np.power(1 + n,3)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2m)/(n + np.power(n,2)) + S1*(29.77777777777778 + 16/np.power(n,2) - 16*S2 + 8*Sp2m) + 2*Sp3m + 10*zeta3 + zeta2*(16*S1 - 16*Sp1m - (16*(1 + n*log2))/n) # pylint: disable=line-too-long gqq1m_cfcf = -32*g3n + (24 - n*(-32 + 3*n*(-8 + n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(-24/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2m) + (8*Sp2m)/(n + np.power(n,2)) - 4*Sp3m - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1m + 32*(1/n + log2)) # pylint: disable=line-too-long @@ -56,7 +56,7 @@ def gamma_nsm(n, nf, sx): @nb.njit(cache=True) -def gamma_nsp(n, nf, sx): +def gamma_nsp(n, nf, sx, sx_new): r""" Compute the |NLO| singlet-like non-singlet anomalous dimension. @@ -79,10 +79,10 @@ def gamma_nsp(n, nf, sx): """ S1 = sx[0] S2 = sx[1] - Sp1p = harmonics.S1(n / 2) - Sp2p = harmonics.S2(n / 2) - Sp3p = harmonics.S3(n / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + Sp1p = sx_new[0] + Sp2p = sx_new[1] + Sp3p = sx_new[2] + g3n = sx_new[3] # fmt: off gqq1p_cfca = -16*g3n + (132 - n*(340 + n*(655 + 51*n*(2 + n))))/(18.*np.power(n,2)*np.power(1 + n,2)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2p)/(n + np.power(n,2)) + S1*(29.77777777777778 - 16/np.power(n,2) - 16*S2 + 8*Sp2p) + 2*Sp3p + 10*zeta3 + zeta2*(16*S1 - 16*Sp1p + 16*(1/n - log2)) # pylint: disable=line-too-long gqq1p_cfcf = 32*g3n - (8 + n*(32 + n*(40 + 3*n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(40/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2p) + (8*Sp2p)/(n + np.power(n,2)) - 4*Sp3p - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1p + 32*(-(1/n) + log2)) # pylint: disable=line-too-long @@ -124,7 +124,7 @@ def gamma_ps(n, nf): @nb.njit(cache=True) -def gamma_qg(n, nf, sx): +def gamma_qg(n, nf, sx, sx_new): r""" Compute the |NLO| quark-gluon singlet anomalous dimension. @@ -147,7 +147,7 @@ def gamma_qg(n, nf, sx): """ S1 = sx[0] S2 = sx[1] - Sp2p = harmonics.S2(n / 2) + Sp2p = sx_new[1] # fmt: off gqg1_nfca = (-4*(16 + n*(64 + n*(104 + n*(128 + n*(85 + n*(36 + n*(25 + n*(15 + n*(6 + n))))))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) - (16*(3 + 2*n)*S1)/np.power(2 + 3*n + np.power(n,2),2) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) - (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long gqg1_nfcf = (-2*(4 + n*(8 + n*(1 + n)*(25 + n*(26 + 5*n*(2 + n))))))/(np.power(n,3)*np.power(1 + n,3)*(2 + n)) + (8*S1)/np.power(n,2) - (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long @@ -159,7 +159,7 @@ def gamma_qg(n, nf, sx): @nb.njit(cache=True) -def gamma_gq(n, nf, sx): +def gamma_gq(n, nf, sx, sx_new): r""" Compute the |NLO| gluon-quark singlet anomalous dimension. @@ -182,7 +182,7 @@ def gamma_gq(n, nf, sx): """ S1 = sx[0] S2 = sx[1] - Sp2p = harmonics.S2(n / 2) + Sp2p = sx_new[1] # fmt: off ggq1_cfcf = (-8 + 2*n*(-12 + n*(-1 + n*(28 + n*(43 + 6*n*(5 + 2*n))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)) - (4*(10 + n*(17 + n*(8 + 5*n)))*S1)/((-1 + n)*n*np.power(1 + n,2)) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(-1 + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long ggq1_cfca = (-4*(144 + n*(432 + n*(-152 + n*(-1304 + n*(-1031 + n*(695 + n*(1678 + n*(1400 + n*(621 + 109*n))))))))))/(9.*np.power(n,3)*np.power(1 + n,3)*np.power(-2 + n + np.power(n,2),2)) + (4*(-12 + n*(-22 + 41*n + 17*np.power(n,3)))*S1)/(3.*np.power(-1 + n,2)*np.power(n,2)*(1 + n)) + ((8 + 4*n + 4*np.power(n,2))*np.power(S1,2))/(n - np.power(n,3)) + ((8 + 4*n + 4*np.power(n,2))*S2)/(n - np.power(n,3)) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long @@ -197,7 +197,7 @@ def gamma_gq(n, nf, sx): @nb.njit(cache=True) -def gamma_gg(n, nf, sx): +def gamma_gg(n, nf, sx, sx_new): r""" Compute the |NLO| gluon-gluon singlet anomalous dimension. @@ -219,10 +219,10 @@ def gamma_gg(n, nf, sx): :math:`\\gamma_{gg}^{(1)}(N)` """ S1 = sx[0] - Sp1p = harmonics.S1(n / 2) - Sp2p = harmonics.S2(n / 2) - Sp3p = harmonics.S3(n / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + Sp1p = sx_new[0] + Sp2p = sx_new[1] + Sp3p = sx_new[2] + g3n = sx_new[3] # fmt: off ggg1_caca = 16*g3n - (2*(576 + n*(1488 + n*(560 + n*(-1248 + n*(-1384 + n*(1663 + n*(4514 + n*(4744 + n*(3030 + n*(1225 + 48*n*(7 + n))))))))))))/(9.*np.power(-1 + n,2)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) + S1*(29.77777777777778 + 16/np.power(-1 + n,2) + 16/np.power(1 + n,2) - 16/np.power(2 + n,2) - 8*Sp2p) + (16*(1 + n + np.power(n,2))*Sp2p)/(n*(1 + n)*(-2 + n + np.power(n,2))) - 2*Sp3p - 10*zeta3 + zeta2*(-16*S1 + 16*Sp1p + 16*(-(1/n) + log2)) # pylint: disable=line-too-long ggg1_canf = (8*(6 + n*(1 + n)*(28 + n*(1 + n)*(13 + 3*n*(1 + n)))))/(9.*np.power(n,2)*np.power(1 + n,2)*(-2 + n + np.power(n,2))) - (40*S1)/9. # pylint: disable=line-too-long @@ -235,7 +235,7 @@ def gamma_gg(n, nf, sx): @nb.njit(cache=True) -def gamma_singlet(n, nf, sx): +def gamma_singlet(n, nf, sx, sx_new): r""" Compute the next-leading-order singlet anomalous dimension matrix. @@ -267,16 +267,19 @@ def gamma_singlet(n, nf, sx): gamma_gq : :math:`\gamma_{gq}^{(1)}` gamma_gg : :math:`\gamma_{gg}^{(1)}` """ - gamma_qq = gamma_nsp(n, nf, sx) + gamma_ps(n, nf) + gamma_qq = gamma_nsp(n, nf, sx, sx_new) + gamma_ps(n, nf) gamma_S_0 = np.array( - [[gamma_qq, gamma_qg(n, nf, sx)], [gamma_gq(n, nf, sx), gamma_gg(n, nf, sx)]], + [ + [gamma_qq, gamma_qg(n, nf, sx, sx_new)], + [gamma_gq(n, nf, sx, sx_new), gamma_gg(n, nf, sx, sx_new)], + ], np.complex_, ) return gamma_S_0 @nb.njit(cache=True) -def gamma_singlet_qed(N, nf, sx): +def gamma_singlet_qed(N, nf, sx, sx_new): r""" Compute the leading-order singlet anomalous dimension matrix. @@ -309,13 +312,18 @@ def gamma_singlet_qed(N, nf, sx): gamma_gq : :math:`\gamma_{gq}^{(0)}` gamma_gg : :math:`\gamma_{gg}^{(0)}` """ - gamma_ns_p = gamma_nsp(N, nf, sx) + gamma_ns_p = gamma_nsp(N, nf, sx, sx_new) gamma_qq = gamma_ns_p + gamma_ps(N, nf) gamma_S = np.array( [ - [gamma_gg(N, nf, sx), 0.0 + 0.0j, gamma_gq(N, nf, sx), 0.0 + 0.0j], + [ + gamma_gg(N, nf, sx, sx_new), + 0.0 + 0.0j, + gamma_gq(N, nf, sx, sx_new), + 0.0 + 0.0j, + ], [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j], - [gamma_qg(N, nf, sx), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j], + [gamma_qg(N, nf, sx, sx_new), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j], [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_ns_p], ], np.complex_, @@ -324,7 +332,7 @@ def gamma_singlet_qed(N, nf, sx): @nb.njit(cache=True) -def gamma_valence_qed(N, nf, sx): +def gamma_valence_qed(N, nf, sx, sx_new): r""" Compute the leading-order valence anomalous dimension matrix. @@ -362,4 +370,4 @@ def gamma_valence_qed(N, nf, sx): ], np.complex_, ) - return gamma_V * gamma_nsm(N, nf, sx) + return gamma_V * gamma_nsm(N, nf, sx, sx_new) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index b284c9f39..cbc84b80a 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -221,7 +221,7 @@ def compute_cache(n, max_weight, is_singlet): @nb.njit(cache=True) -def compute_qed_ns_cache(n, s1): +def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name r"""Get the harmonics sums cache needed for the qed non singlet AD. Parameters @@ -249,21 +249,29 @@ def compute_qed_ns_cache(n, s1): """ s1h = S1(n / 2.0) sx_qed_ns = [s1h] - S2h = S2(n / 2.0) - sx_qed_ns.append(S2h) - S3h = S3(n / 2.0) - sx_qed_ns.append(S3h) - S1p1h = S1((n + 1.0) / 2.0) - sx_qed_ns.append(S1p1h) - S2p1h = S2((n + 1.0) / 2.0) - sx_qed_ns.append(S2p1h) - S3p1h = S3((n + 1.0) / 2.0) - sx_qed_ns.append(S3p1h) + s2h = S2(n / 2.0) + sx_qed_ns.append(s2h) + s3h = S3(n / 2.0) + sx_qed_ns.append(s3h) g3N = g_functions.mellin_g3(n, s1) sx_qed_ns.append(g3N) - S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) - g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) - sx_qed_ns.append(g3Np2) + if not qcd_singlet: + s1_nm1_h = S1((n - 1.0) / 2.0) + sx_qed_ns.append(s1_nm1_h) + s2_nm1_h = S2((n - 1.0) / 2.0) + sx_qed_ns.append(s2_nm1_h) + s3_nm1_h = S3((n - 1.0) / 2.0) + sx_qed_ns.append(s3_nm1_h) + if qed: + S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) + g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) + sx_qed_ns.append(g3Np2) + s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1) / 2, 1, 1) + sx_qed_ns.append(s1_np1_h) + s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1) / 2, 1, 2) + sx_qed_ns.append(s2_np1_h) + s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1) / 2, 1, 3) + sx_qed_ns.append(s3_np1_h) return sx_qed_ns diff --git a/tests/eko/test_ad.py b/tests/eko/test_ad.py index 4993bfc8f..80f93d843 100644 --- a/tests/eko/test_ad.py +++ b/tests/eko/test_ad.py @@ -191,11 +191,12 @@ def test_dim_singlet(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) + sx_new = harmonics.compute_new_sx_cache(N, sx[0], False, True) gamma_singlet = ad.gamma_singlet_qed((3, 2), N, nf) assert gamma_singlet.shape == (4, 3, 4, 4) gamma_singlet_as1 = as1.gamma_singlet_qed(N, sx[0], nf) assert gamma_singlet_as1.shape == (4, 4) - gamma_singlet_as2 = as2.gamma_singlet_qed(N, nf, sx) + gamma_singlet_as2 = as2.gamma_singlet_qed(N, nf, sx, sx_new) assert gamma_singlet_as2.shape == (4, 4) gamma_singlet_as3 = as3.gamma_singlet_qed(N, nf, sx) assert gamma_singlet_as3.shape == (4, 4) @@ -205,11 +206,12 @@ def test_dim_valence(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) + sx_new = harmonics.compute_new_sx_cache(N, sx[0], False, False) gamma_valence = ad.gamma_valence_qed((3, 2), N, nf) assert gamma_valence.shape == (4, 3, 2, 2) gamma_valence_as1 = as1.gamma_valence_qed(N, sx[0]) assert gamma_valence_as1.shape == (2, 2) - gamma_valence_as2 = as2.gamma_valence_qed(N, nf, sx) + gamma_valence_as2 = as2.gamma_valence_qed(N, nf, sx, sx_new) assert gamma_valence_as2.shape == (2, 2) gamma_valence_as3 = as3.gamma_valence_qed(N, nf, sx) assert gamma_valence_as3.shape == (2, 2) diff --git a/tests/eko/test_ad_aem2.py b/tests/eko/test_ad_aem2.py index 362b46306..18b5e6cbb 100644 --- a/tests/eko/test_ad_aem2.py +++ b/tests/eko/test_ad_aem2.py @@ -10,7 +10,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_qed_ns_cache(N, sx[0]) + sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) for NF in range(2, 6 + 1): np.testing.assert_almost_equal( ad.aem2.gamma_nsmu(N, NF, sx, sx_ns_qed), 0, decimal=4 @@ -39,7 +39,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_qed_ns_cache(N, sx[0]) + sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) NF = 6 NU = constants.uplike_flavors(NF) ND = NF - NU diff --git a/tests/eko/test_ad_as1aem1.py b/tests/eko/test_ad_as1aem1.py index c0db6dd6b..dc14d8a24 100644 --- a/tests/eko/test_ad_as1aem1.py +++ b/tests/eko/test_ad_as1aem1.py @@ -11,7 +11,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_qed_ns_cache(N, sx[0]) + sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) np.testing.assert_almost_equal(ad.as1aem1.gamma_nsm(N, sx, sx_ns_qed), 0, decimal=4) @@ -53,7 +53,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_qed_ns_cache(N, sx[0]) + sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) np.testing.assert_almost_equal( ad.as1aem1.gamma_nsp(N, sx, sx_ns_qed) + ad.as1aem1.gamma_gq(N, sx) diff --git a/tests/eko/test_ad_as2.py b/tests/eko/test_ad_as2.py index 01d7f2f28..89c009e04 100644 --- a/tests/eko/test_ad_as2.py +++ b/tests/eko/test_ad_as2.py @@ -11,10 +11,12 @@ def test_gamma_1(): # number conservation sx_n1 = h.sx(1, 2) - np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1), 0.0, atol=2e-6) + sx_new1 = h.compute_new_sx_cache(1, sx_n1[0], False, False) + np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1, sx_new1), 0.0, atol=2e-6) sx_n2 = h.sx(2, 2) - gS1 = ad_as2.gamma_singlet(2, NF, sx_n2) + sx_new2 = h.compute_new_sx_cache(2, sx_n2[0], True, False) + gS1 = ad_as2.gamma_singlet(2, NF, sx_n2, sx_new2) # 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) @@ -26,7 +28,7 @@ def test_gamma_1(): # reference values are obtained from MMa # non-singlet sector np.testing.assert_allclose( - ad_as2.gamma_nsp(2, NF, sx_n2), + ad_as2.gamma_nsp(2, NF, sx_n2, sx_new2), (-112.0 * const.CF + 376.0 * const.CA - 64.0 * NF) * const.CF / 27.0, ) # singlet sector @@ -39,8 +41,9 @@ def test_gamma_1(): ) # gq # add additional point at (analytical) continuation point + sx_new2 = h.compute_new_sx_cache(2, sx_n2[0], False, False) np.testing.assert_allclose( - ad_as2.gamma_nsm(2, NF, sx_n2), + ad_as2.gamma_nsm(2, NF, sx_n2, sx_new2), ( (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.constants.zeta3) * const.CF @@ -52,8 +55,10 @@ def test_gamma_1(): ) sx_n3 = h.sx(3, 2) sx_n4 = h.sx(4, 2) + sx_new3 = h.compute_new_sx_cache(3, sx_n3[0], True, False) + sx_new4 = h.compute_new_sx_cache(4, sx_n4[0], True, False) np.testing.assert_allclose( - ad_as2.gamma_nsp(3, NF, sx_n3), + ad_as2.gamma_nsp(3, NF, sx_n3, sx_new3), ( (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.constants.zeta3) * const.CF @@ -64,7 +69,7 @@ def test_gamma_1(): * const.CF, ) np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * const.CF * NF / 5400.0) - gS1 = ad_as2.gamma_singlet(3, NF, sx_n3) + gS1 = ad_as2.gamma_singlet(3, NF, sx_n3, sx_new3) np.testing.assert_allclose( gS1[1, 0], ( @@ -84,7 +89,7 @@ def test_gamma_1(): ), rtol=6e-7, ) # gg - gS1 = ad_as2.gamma_singlet(4, NF, sx_n4) + gS1 = ad_as2.gamma_singlet(4, NF, sx_n4, sx_new4) np.testing.assert_allclose( gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF ) # qg @@ -92,7 +97,7 @@ def test_gamma_1(): const.update_colors(4) np.testing.assert_allclose(const.CA, 4.0) - gS1 = ad_as2.gamma_singlet(3, NF, sx_n3) + gS1 = ad_as2.gamma_singlet(3, NF, sx_n3, sx_new3) np.testing.assert_allclose( gS1[1, 0], ( @@ -112,7 +117,7 @@ def test_gamma_1(): ), rtol=6e-7, ) # gg - gS1 = ad_as2.gamma_singlet(4, NF, sx_n4) + gS1 = ad_as2.gamma_singlet(4, NF, sx_n4, sx_new4) np.testing.assert_allclose( gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF ) # qg From 859afd4cb81e995b986341fd84ebfa07c71350b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 2 Dec 2022 15:00:48 +0100 Subject: [PATCH 03/14] Change name of sx_cache_qed --- src/eko/anomalous_dimensions/__init__.py | 20 ++++++++++---------- src/eko/harmonics/__init__.py | 9 ++++++--- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index d0fd75d46..c4b323f21 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -337,7 +337,7 @@ def choose_ns_ad_aem1(mode, n, sx): @nb.njit(cache=True) -def choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed): +def choose_ns_ad_as1aem1(mode, n, sx, sx_new): r""" Select the non-singlet anomalous dimension at O(as1aem1) with the correct charge factor. @@ -356,17 +356,17 @@ def choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed): non-singlet anomalous dimensions """ if mode == 10102: - return constants.eu2 * as1aem1.gamma_nsp(n, sx, sx_ns_qed) + return constants.eu2 * as1aem1.gamma_nsp(n, sx, sx_new) elif mode == 10103: - return constants.ed2 * as1aem1.gamma_nsp(n, sx, sx_ns_qed) + return constants.ed2 * as1aem1.gamma_nsp(n, sx, sx_new) elif mode == 10202: - return constants.eu2 * as1aem1.gamma_nsm(n, sx, sx_ns_qed) + return constants.eu2 * as1aem1.gamma_nsm(n, sx, sx_new) elif mode == 10203: - return constants.ed2 * as1aem1.gamma_nsm(n, sx, sx_ns_qed) + return constants.ed2 * as1aem1.gamma_nsm(n, sx, sx_new) @nb.njit(cache=True) -def choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed): +def choose_ns_ad_aem2(mode, n, nf, sx, sx_new): r""" Select the non-singlet anomalous dimension at O(aem2) with the correct charge factor. @@ -385,13 +385,13 @@ def choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed): non-singlet anomalous dimensions """ if mode == 10102: - return constants.eu2 * aem2.gamma_nspu(n, nf, sx, sx_ns_qed) + return constants.eu2 * aem2.gamma_nspu(n, nf, sx, sx_new) elif mode == 10103: - return constants.ed2 * aem2.gamma_nspd(n, nf, sx, sx_ns_qed) + return constants.ed2 * aem2.gamma_nspd(n, nf, sx, sx_new) elif mode == 10202: - return constants.eu2 * aem2.gamma_nsmu(n, nf, sx, sx_ns_qed) + return constants.eu2 * aem2.gamma_nsmu(n, nf, sx, sx_new) elif mode == 10203: - return constants.ed2 * aem2.gamma_nsmd(n, nf, sx, sx_ns_qed) + return constants.ed2 * aem2.gamma_nsmd(n, nf, sx, sx_new) @nb.njit(cache=True) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index cbc84b80a..b892022ee 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -240,11 +240,14 @@ def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name [S_1(n/2), S_2(n/2), S_{3}(n/2), + g_3(n), + S_1((n-1)/2), + S_2((n-1)/2), + S_{3}((n-1)/2), + g_3(n+2), S_1((n+1)/2), S_2((n+1)/2), - S_{3}((n+1)/2), - g_3(n), - g_3(n+2)] + S_{3}((n+1)/2),] """ s1h = S1(n / 2.0) From 342de1b22caf02b813cd1a38bf659a8e293a5a81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 2 Dec 2022 15:15:48 +0100 Subject: [PATCH 04/14] Add description for compute_new_sx_cache --- src/eko/harmonics/__init__.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index b892022ee..fad651721 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -224,6 +224,15 @@ def compute_cache(n, max_weight, is_singlet): def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name r"""Get the harmonics sums cache needed for the qed non singlet AD. + We have three cases : + - 1 : for qcd singlet we need only S_1(n/2), S_2(n/2), S_{3}(n/2) and g_3(n) + - 2 : for qcd non singlet we need also S_1((n-1)/2), S_2((n-1)/2), S_{3}((n-1)/2) + - 3 : for qed we need also g_3(n+2), S_1((n-1)/2), S_2((n-1)/2), S_{3}((n-1)/2) + These three cases are called with (qcd_singlet, qed) = : + - 1 : (True, False) + - 2 : (False, False) + - 3 : (False, True) + Parameters ---------- n : complex From 9db86a03f18afc0cfddc69c767cfb5ffaa00e533 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 2 Dec 2022 15:17:11 +0100 Subject: [PATCH 05/14] Remove unnecessary functions --- src/eko/harmonics/__init__.py | 93 ----------------------------------- 1 file changed, 93 deletions(-) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index fad651721..c86d42940 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -285,96 +285,3 @@ def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1) / 2, 1, 3) sx_qed_ns.append(s3_np1_h) return sx_qed_ns - - -@nb.njit(cache=True) -def compute_cache_S_n_half(n): - r"""Get the harmonics sums cache S_i(n/2). - - Parameters - ---------- - n : complex - Mellin moment - - Returns - ------- - list - harmonic sums cache. It contains: - - .. math :: - [S_1(n/2), S_2(n/2), S_3(n/2)] - """ - return [S1(n / 2.0), S2(n / 2.0), S3(n / 2.0)] - - -@nb.njit(cache=True) -def compute_cache_S_n_half_minus_1halph(n): - r"""Get the harmonics sums cache S_i((n-1)/2). - - Parameters - ---------- - n : complex - Mellin moment - - Returns - ------- - list - harmonic sums cache. It contains: - - .. math :: - [S_1((n-1)/2), S_2((n-1)/2), S_3((n-1)/2)] - """ - return [S1((n - 1) / 2), S2((n - 1) / 2), S3((n - 1) / 2)] - - -@nb.njit(cache=True) -def compute_cache_S_n_half_plus_1halph(n, sx): - r"""Get the harmonics sums cache S_i((n+1)/2). - - Parameters - ---------- - n : complex - Mellin moment - sx : list - harmonic sums cache S_i((n-1)/2) - - Returns - ------- - list - harmonic sums cache. It contains: - - .. math :: - [S_1((n+1)/2), S_2((n+1)/2), S_3((n+1)/2)] - """ - S1_nh = polygamma.recursive_harmonic_sum(sx[0], n, 1, 1) - S2_nh = polygamma.recursive_harmonic_sum(sx[1], n, 1, 2) - S3_nh = polygamma.recursive_harmonic_sum(sx[2], n, 1, 3) - return [S1_nh, S2_nh, S3_nh] - - -@nb.njit(cache=True) -def compute_cache_g3(n, s1, qed): - r"""Get the functionis g3(n) and g3(n+2). - - Parameters - ---------- - n : complex - Mellin moment - s1 : list - harmonic sum :math:`S_1(n)` - - - Returns - ------- - list - g3 functions cache. It contains: - - .. math :: - [g3(n), g3(n+2)] - """ - g_list = [g_functions.mellin_g3(n, s1)] - if qed: - S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) - g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) - g_list.append(g3Np2) - return g_list From e9988ddd8ddf97e6ef70b27fd22c03b0a7860a52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Thu, 8 Dec 2022 12:40:31 +0100 Subject: [PATCH 06/14] Fix docstring of compute_new_sx_cache --- src/eko/harmonics/__init__.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index c86d42940..271315479 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -239,6 +239,10 @@ def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name Mellin moment s1 : float harmonic sum :math:`S_1(N)` + qcd_singlet : bool + True if singlet for pure QCD evolution + qed : bool + True if QED evolution activated Returns ------- @@ -256,7 +260,7 @@ def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name g_3(n+2), S_1((n+1)/2), S_2((n+1)/2), - S_{3}((n+1)/2),] + S_{3}((n+1)/2)] """ s1h = S1(n / 2.0) @@ -278,10 +282,10 @@ def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) sx_qed_ns.append(g3Np2) - s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1) / 2, 1, 1) + s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1.0) / 2.0, 1, 1) sx_qed_ns.append(s1_np1_h) - s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1) / 2, 1, 2) + s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1.0) / 2.0, 1, 2) sx_qed_ns.append(s2_np1_h) - s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1) / 2, 1, 3) + s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1.0) / 2.0, 1, 3) sx_qed_ns.append(s3_np1_h) return sx_qed_ns From 7a1027418328405f69542f51da2a3187725deed2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 9 Dec 2022 16:33:04 +0100 Subject: [PATCH 07/14] Rename compute_new_sx_cache -> compute_additional_sx_cache --- src/eko/anomalous_dimensions/__init__.py | 10 +++++----- src/eko/harmonics/__init__.py | 2 +- tests/eko/test_ad.py | 4 ++-- tests/eko/test_ad_aem2.py | 4 ++-- tests/eko/test_ad_as1aem1.py | 4 ++-- tests/eko/test_ad_as2.py | 10 +++++----- 6 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index c4b323f21..881d2031f 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -152,7 +152,7 @@ def gamma_ns(order, mode, n, nf): gamma_ns[0] = as1.gamma_ns(n, sx[0]) # NLO and beyond if order[0] >= 2: - sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, False) + sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, False) if mode == 10101: gamma_ns_1 = as2.gamma_nsp(n, nf, sx, sx_new) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here @@ -228,7 +228,7 @@ def gamma_singlet(order, n, nf): gamma_s = np.zeros((order[0], 2, 2), np.complex_) gamma_s[0] = as1.gamma_singlet(n, sx[0], nf) if order[0] >= 2: - sx_new = harmonics.compute_new_sx_cache(n, sx[0], True, False) + sx_new = harmonics.compute_additional_sx_cache(n, sx[0], True, False) gamma_s[1] = as2.gamma_singlet(n, nf, sx, sx_new) if order[0] >= 3: gamma_s[2] = as3.gamma_singlet(n, nf, sx) @@ -282,7 +282,7 @@ def gamma_ns_qed(order, mode, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) # now combine gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex_) if order[0] >= 1: @@ -429,7 +429,7 @@ def gamma_singlet_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex_) if order[0] >= 1: gamma_s[1, 0] = as1.gamma_singlet_qed(n, sx[0], nf) @@ -481,7 +481,7 @@ def gamma_valence_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_new_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex_) if order[0] >= 1: gamma_v[1, 0] = as1.gamma_valence_qed(n, sx[0]) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index 271315479..22e2450b8 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -221,7 +221,7 @@ def compute_cache(n, max_weight, is_singlet): @nb.njit(cache=True) -def compute_new_sx_cache(n, s1, qcd_singlet, qed): # TODO : change name +def compute_additional_sx_cache(n, s1, qcd_singlet, qed): r"""Get the harmonics sums cache needed for the qed non singlet AD. We have three cases : diff --git a/tests/eko/test_ad.py b/tests/eko/test_ad.py index 80f93d843..93d21d94f 100644 --- a/tests/eko/test_ad.py +++ b/tests/eko/test_ad.py @@ -191,7 +191,7 @@ def test_dim_singlet(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) - sx_new = harmonics.compute_new_sx_cache(N, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(N, sx[0], False, True) gamma_singlet = ad.gamma_singlet_qed((3, 2), N, nf) assert gamma_singlet.shape == (4, 3, 4, 4) gamma_singlet_as1 = as1.gamma_singlet_qed(N, sx[0], nf) @@ -206,7 +206,7 @@ def test_dim_valence(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) - sx_new = harmonics.compute_new_sx_cache(N, sx[0], False, False) + sx_new = harmonics.compute_additional_sx_cache(N, sx[0], False, False) gamma_valence = ad.gamma_valence_qed((3, 2), N, nf) assert gamma_valence.shape == (4, 3, 2, 2) gamma_valence_as1 = as1.gamma_valence_qed(N, sx[0]) diff --git a/tests/eko/test_ad_aem2.py b/tests/eko/test_ad_aem2.py index 18b5e6cbb..2661e7709 100644 --- a/tests/eko/test_ad_aem2.py +++ b/tests/eko/test_ad_aem2.py @@ -10,7 +10,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) for NF in range(2, 6 + 1): np.testing.assert_almost_equal( ad.aem2.gamma_nsmu(N, NF, sx, sx_ns_qed), 0, decimal=4 @@ -39,7 +39,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) NF = 6 NU = constants.uplike_flavors(NF) ND = NF - NU diff --git a/tests/eko/test_ad_as1aem1.py b/tests/eko/test_ad_as1aem1.py index dc14d8a24..9de4cddd2 100644 --- a/tests/eko/test_ad_as1aem1.py +++ b/tests/eko/test_ad_as1aem1.py @@ -11,7 +11,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) np.testing.assert_almost_equal(ad.as1aem1.gamma_nsm(N, sx, sx_ns_qed), 0, decimal=4) @@ -53,7 +53,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_new_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) np.testing.assert_almost_equal( ad.as1aem1.gamma_nsp(N, sx, sx_ns_qed) + ad.as1aem1.gamma_gq(N, sx) diff --git a/tests/eko/test_ad_as2.py b/tests/eko/test_ad_as2.py index 89c009e04..72dc9b32b 100644 --- a/tests/eko/test_ad_as2.py +++ b/tests/eko/test_ad_as2.py @@ -11,11 +11,11 @@ def test_gamma_1(): # number conservation sx_n1 = h.sx(1, 2) - sx_new1 = h.compute_new_sx_cache(1, sx_n1[0], False, False) + sx_new1 = h.compute_additional_sx_cache(1, sx_n1[0], False, False) np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1, sx_new1), 0.0, atol=2e-6) sx_n2 = h.sx(2, 2) - sx_new2 = h.compute_new_sx_cache(2, sx_n2[0], True, False) + sx_new2 = h.compute_additional_sx_cache(2, sx_n2[0], True, False) gS1 = ad_as2.gamma_singlet(2, NF, sx_n2, sx_new2) # gluon momentum conservation # the CA*NF term seems to be tough to compute, so raise the constraint ... @@ -41,7 +41,7 @@ def test_gamma_1(): ) # gq # add additional point at (analytical) continuation point - sx_new2 = h.compute_new_sx_cache(2, sx_n2[0], False, False) + sx_new2 = h.compute_additional_sx_cache(2, sx_n2[0], False, False) np.testing.assert_allclose( ad_as2.gamma_nsm(2, NF, sx_n2, sx_new2), ( @@ -55,8 +55,8 @@ def test_gamma_1(): ) sx_n3 = h.sx(3, 2) sx_n4 = h.sx(4, 2) - sx_new3 = h.compute_new_sx_cache(3, sx_n3[0], True, False) - sx_new4 = h.compute_new_sx_cache(4, sx_n4[0], True, False) + sx_new3 = h.compute_additional_sx_cache(3, sx_n3[0], True, False) + sx_new4 = h.compute_additional_sx_cache(4, sx_n4[0], True, False) np.testing.assert_allclose( ad_as2.gamma_nsp(3, NF, sx_n3, sx_new3), ( From 63d800a44d4616e9923f4a53be8c5d39a7648b28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Fri, 9 Dec 2022 16:50:19 +0100 Subject: [PATCH 08/14] Fix benchmarks/eko/benchmark_ad.py --- benchmarks/eko/benchmark_ad.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py index cb5c07e66..9e86323e2 100644 --- a/benchmarks/eko/benchmark_ad.py +++ b/benchmarks/eko/benchmark_ad.py @@ -138,8 +138,9 @@ def check_gamma_1_pegasus(N, NF): P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC) sx = h.sx(N, 2) - np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx), -P1NSP) - np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx), -P1NSM) + sx_new = h.compute_additional_sx_cache(N, sx[0], False, False) + np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx, sx_new), -P1NSP) + np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx, sx_new), -P1NSM) NS = N * N NT = NS * N @@ -256,7 +257,8 @@ 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_as2.gamma_singlet(N, NF, sx) + sx_new = h.compute_additional_sx_cache(N, sx[0], True, False) + gS1 = ad_as2.gamma_singlet(N, NF, sx, sx_new) np.testing.assert_allclose(gS1[0, 0], -P1Sqq) np.testing.assert_allclose(gS1[0, 1], -P1Sqg) np.testing.assert_allclose(gS1[1, 0], -P1Sgq) From fcf874451ec127e1b813dbea683db3d49d9efa06 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Sun, 18 Dec 2022 21:03:22 +0100 Subject: [PATCH 09/14] Add additional_sx_cache to compute_cache --- src/eko/anomalous_dimensions/__init__.py | 22 +++++--- src/eko/harmonics/__init__.py | 70 ++++++++++++++++-------- tests/eko/test_ad.py | 4 +- tests/eko/test_ad_aem2.py | 4 +- tests/eko/test_ad_as1aem1.py | 4 +- tests/eko/test_ad_as2.py | 10 ++-- 6 files changed, 71 insertions(+), 43 deletions(-) diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 23e8dec63..235359780 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -135,7 +135,9 @@ def gamma_ns(order, mode, n, nf): """ # cache the s-es if order[0] >= 4: - full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False) + full_sx_cache = harmonics.compute_cache( + n, 5, is_singlet=False, additional=True, qed=False + ) sx = np.array( [ full_sx_cache[0][0], @@ -144,14 +146,15 @@ def gamma_ns(order, mode, n, nf): full_sx_cache[3][0], ] ) + sx_new = full_sx_cache[-1] else: sx = harmonics.sx(n, max_weight=order[0] + 1) + sx_new = harmonics.compute_additional_sx_cache(n, False, False, sx) # now combine gamma_ns = np.zeros(order[0], np.complex_) gamma_ns[0] = as1.gamma_ns(n, sx[0]) # NLO and beyond if order[0] >= 2: - sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, False) if mode == 10101: gamma_ns_1 = as2.gamma_nsp(n, nf, sx, sx_new) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here @@ -209,7 +212,9 @@ def gamma_singlet(order, n, nf): """ # cache the s-es if order[0] >= 4: - full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False) + full_sx_cache = harmonics.compute_cache( + n, 5, is_singlet=False, additional=True, qed=False + ) sx = np.array( [ full_sx_cache[0][0], @@ -218,16 +223,17 @@ def gamma_singlet(order, n, nf): full_sx_cache[3][0], ] ) + sx_new = full_sx_cache[-1] elif order[0] >= 3: # here we need only S1,S2,S3,S4 sx = harmonics.sx(n, max_weight=order[0] + 1) + sx_new = harmonics.compute_additional_sx_cache(n, True, False, sx) else: sx = harmonics.sx(n, max_weight=order[0]) - + sx_new = harmonics.compute_additional_sx_cache(n, True, False, sx) gamma_s = np.zeros((order[0], 2, 2), np.complex_) gamma_s[0] = as1.gamma_singlet(n, sx[0], nf) if order[0] >= 2: - sx_new = harmonics.compute_additional_sx_cache(n, sx[0], True, False) gamma_s[1] = as2.gamma_singlet(n, nf, sx, sx_new) if order[0] >= 3: gamma_s[2] = as3.gamma_singlet(n, nf, sx) @@ -281,7 +287,7 @@ def gamma_ns_qed(order, mode, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx) # now combine gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex_) if order[0] >= 1: @@ -428,7 +434,7 @@ def gamma_singlet_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx) gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex_) if order[0] >= 1: gamma_s[1, 0] = as1.gamma_singlet_qed(n, sx[0], nf) @@ -480,7 +486,7 @@ def gamma_valence_qed(order, n, nf): sx = harmonics.sx(n, max_weight=max_weight + 1) else: sx = harmonics.sx(n, max_weight=3) - sx_new = harmonics.compute_additional_sx_cache(n, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx) gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex_) if order[0] >= 1: gamma_v[1, 0] = as1.gamma_valence_qed(n, sx[0]) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index 22e2450b8..cc07402d0 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -180,7 +180,7 @@ def s4x(n, sx, smx, is_singlet): @nb.njit(cache=True) -def compute_cache(n, max_weight, is_singlet): +def compute_cache(n, max_weight, is_singlet, additional=False, qed=False): r"""Get the harmonics sums cache. Parameters @@ -217,11 +217,14 @@ def compute_cache(n, max_weight, is_singlet): sx[2, 1:-2] = s3x(n, sx[:, 0], sx[:, -1], is_singlet) sx[3, 1:-1] = s4x(n, sx[:, 0], sx[:, -1], is_singlet) # return list of list keeping the non zero values - return [[el for el in sx_list if el != 0] for sx_list in sx] + tmp = [[el for el in sx_list if el != 0] for sx_list in sx] + if additional: + tmp.append(compute_additional_sx_cache(n, is_singlet, qed, sx[:, 0], sx[:, -1])) + return tmp @nb.njit(cache=True) -def compute_additional_sx_cache(n, s1, qcd_singlet, qed): +def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None): r"""Get the harmonics sums cache needed for the qed non singlet AD. We have three cases : @@ -263,29 +266,48 @@ def compute_additional_sx_cache(n, s1, qcd_singlet, qed): S_{3}((n+1)/2)] """ - s1h = S1(n / 2.0) - sx_qed_ns = [s1h] - s2h = S2(n / 2.0) - sx_qed_ns.append(s2h) - s3h = S3(n / 2.0) - sx_qed_ns.append(s3h) - g3N = g_functions.mellin_g3(n, s1) - sx_qed_ns.append(g3N) - if not qcd_singlet: - s1_nm1_h = S1((n - 1.0) / 2.0) - sx_qed_ns.append(s1_nm1_h) - s2_nm1_h = S2((n - 1.0) / 2.0) - sx_qed_ns.append(s2_nm1_h) - s3_nm1_h = S3((n - 1.0) / 2.0) - sx_qed_ns.append(s3_nm1_h) + if sminus is None: + s1h = S1(n / 2.0) + sx_new = [s1h] + s2h = S2(n / 2.0) + sx_new.append(s2h) + s3h = S3(n / 2.0) + sx_new.append(s3h) + g3N = g_functions.mellin_g3(n, sx[0]) + sx_new.append(g3N) + if not is_singlet: + s1_nm1_h = S1((n - 1.0) / 2.0) + sx_new.append(s1_nm1_h) + s2_nm1_h = S2((n - 1.0) / 2.0) + sx_new.append(s2_nm1_h) + s3_nm1_h = S3((n - 1.0) / 2.0) + sx_new.append(s3_nm1_h) + else: + if is_singlet: + s1h = sminus[0] + sx[0] + sx_new = [s1h] + s2h = 2 * (sminus[1] + sx[1]) + sx_new.append(s2h) + s3h = 4 * (sminus[2] + sx[2]) + sx_new.append(s3h) + g3N = g_functions.mellin_g3(n, sx[0]) + sx_new.append(g3N) + else: + sx_new = compute_additional_sx_cache(n, True, False, sx, None) + s1_nm1_h = sminus[0] + sx[0] + sx_new.append(s1_nm1_h) + s2_nm1_h = 2 * (sminus[1] + sx[1]) + sx_new.append(s2_nm1_h) + s3_nm1_h = 4 * (sminus[2] + sx[2]) + sx_new.append(s3_nm1_h) if qed: - S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) + S1p2 = polygamma.recursive_harmonic_sum(sx[0], n, 2, 1) g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) - sx_qed_ns.append(g3Np2) + sx_new.append(g3Np2) s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1.0) / 2.0, 1, 1) - sx_qed_ns.append(s1_np1_h) + sx_new.append(s1_np1_h) s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1.0) / 2.0, 1, 2) - sx_qed_ns.append(s2_np1_h) + sx_new.append(s2_np1_h) s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1.0) / 2.0, 1, 3) - sx_qed_ns.append(s3_np1_h) - return sx_qed_ns + sx_new.append(s3_np1_h) + return sx_new diff --git a/tests/eko/test_ad.py b/tests/eko/test_ad.py index 93d21d94f..4faa61046 100644 --- a/tests/eko/test_ad.py +++ b/tests/eko/test_ad.py @@ -191,7 +191,7 @@ def test_dim_singlet(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) - sx_new = harmonics.compute_additional_sx_cache(N, sx[0], False, True) + sx_new = harmonics.compute_additional_sx_cache(N, False, True, sx) gamma_singlet = ad.gamma_singlet_qed((3, 2), N, nf) assert gamma_singlet.shape == (4, 3, 4, 4) gamma_singlet_as1 = as1.gamma_singlet_qed(N, sx[0], nf) @@ -206,7 +206,7 @@ def test_dim_valence(): nf = 3 N = 2 sx = harmonics.sx(N, max_weight=3 + 1) - sx_new = harmonics.compute_additional_sx_cache(N, sx[0], False, False) + sx_new = harmonics.compute_additional_sx_cache(N, False, False, sx) gamma_valence = ad.gamma_valence_qed((3, 2), N, nf) assert gamma_valence.shape == (4, 3, 2, 2) gamma_valence_as1 = as1.gamma_valence_qed(N, sx[0]) diff --git a/tests/eko/test_ad_aem2.py b/tests/eko/test_ad_aem2.py index 2661e7709..92ebb0788 100644 --- a/tests/eko/test_ad_aem2.py +++ b/tests/eko/test_ad_aem2.py @@ -10,7 +10,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx) for NF in range(2, 6 + 1): np.testing.assert_almost_equal( ad.aem2.gamma_nsmu(N, NF, sx, sx_ns_qed), 0, decimal=4 @@ -39,7 +39,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx) NF = 6 NU = constants.uplike_flavors(NF) ND = NF - NU diff --git a/tests/eko/test_ad_as1aem1.py b/tests/eko/test_ad_as1aem1.py index 9de4cddd2..58cb11e05 100644 --- a/tests/eko/test_ad_as1aem1.py +++ b/tests/eko/test_ad_as1aem1.py @@ -11,7 +11,7 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx) np.testing.assert_almost_equal(ad.as1aem1.gamma_nsm(N, sx, sx_ns_qed), 0, decimal=4) @@ -53,7 +53,7 @@ def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) sx = h.sx(N, 3) - sx_ns_qed = h.compute_additional_sx_cache(N, sx[0], False, True) + sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx) np.testing.assert_almost_equal( ad.as1aem1.gamma_nsp(N, sx, sx_ns_qed) + ad.as1aem1.gamma_gq(N, sx) diff --git a/tests/eko/test_ad_as2.py b/tests/eko/test_ad_as2.py index 72dc9b32b..6562833ea 100644 --- a/tests/eko/test_ad_as2.py +++ b/tests/eko/test_ad_as2.py @@ -11,11 +11,11 @@ def test_gamma_1(): # number conservation sx_n1 = h.sx(1, 2) - sx_new1 = h.compute_additional_sx_cache(1, sx_n1[0], False, False) + sx_new1 = h.compute_additional_sx_cache(1, False, False, sx_n1) np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1, sx_new1), 0.0, atol=2e-6) sx_n2 = h.sx(2, 2) - sx_new2 = h.compute_additional_sx_cache(2, sx_n2[0], True, False) + sx_new2 = h.compute_additional_sx_cache(2, True, False, sx_n2) gS1 = ad_as2.gamma_singlet(2, NF, sx_n2, sx_new2) # gluon momentum conservation # the CA*NF term seems to be tough to compute, so raise the constraint ... @@ -41,7 +41,7 @@ def test_gamma_1(): ) # gq # add additional point at (analytical) continuation point - sx_new2 = h.compute_additional_sx_cache(2, sx_n2[0], False, False) + sx_new2 = h.compute_additional_sx_cache(2, False, False, sx_n2) np.testing.assert_allclose( ad_as2.gamma_nsm(2, NF, sx_n2, sx_new2), ( @@ -55,8 +55,8 @@ def test_gamma_1(): ) sx_n3 = h.sx(3, 2) sx_n4 = h.sx(4, 2) - sx_new3 = h.compute_additional_sx_cache(3, sx_n3[0], True, False) - sx_new4 = h.compute_additional_sx_cache(4, sx_n4[0], True, False) + sx_new3 = h.compute_additional_sx_cache(3, True, False, sx_n3) + sx_new4 = h.compute_additional_sx_cache(4, True, False, sx_n4) np.testing.assert_allclose( ad_as2.gamma_nsp(3, NF, sx_n3, sx_new3), ( From f0283bc952de642c7382f406180ca3fb71151dcf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Mon, 19 Dec 2022 15:30:47 +0100 Subject: [PATCH 10/14] Fix benchmark_ad.py --- benchmarks/eko/benchmark_ad.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py index 9e86323e2..f8ba47bd2 100644 --- a/benchmarks/eko/benchmark_ad.py +++ b/benchmarks/eko/benchmark_ad.py @@ -138,7 +138,7 @@ def check_gamma_1_pegasus(N, NF): P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC) sx = h.sx(N, 2) - sx_new = h.compute_additional_sx_cache(N, sx[0], False, False) + sx_new = h.compute_additional_sx_cache(N, False, False, sx) np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx, sx_new), -P1NSP) np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx, sx_new), -P1NSM) @@ -257,7 +257,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 - sx_new = h.compute_additional_sx_cache(N, sx[0], True, False) + sx_new = h.compute_additional_sx_cache(N, True, False, sx) gS1 = ad_as2.gamma_singlet(N, NF, sx, sx_new) np.testing.assert_allclose(gS1[0, 0], -P1Sqq) np.testing.assert_allclose(gS1[0, 1], -P1Sqg) From 9338c4cd1e074ffbf1a2df6e2daac41ece229113 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Mon, 19 Dec 2022 16:07:20 +0100 Subject: [PATCH 11/14] Implement g3 as a function of Sm21 --- src/eko/harmonics/__init__.py | 98 +++++++++++++++++++++++++---------- 1 file changed, 72 insertions(+), 26 deletions(-) diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py index cc07402d0..5d07de63d 100644 --- a/src/eko/harmonics/__init__.py +++ b/src/eko/harmonics/__init__.py @@ -6,6 +6,7 @@ import numpy as np from . import g_functions, polygamma +from .constants import log2, zeta2, zeta3 from .w1 import S1, Sm1 from .w2 import S2, Sm2 from .w3 import S3, S21, S2m1, Sm2m1, Sm3, Sm21 @@ -219,12 +220,21 @@ def compute_cache(n, max_weight, is_singlet, additional=False, qed=False): # return list of list keeping the non zero values tmp = [[el for el in sx_list if el != 0] for sx_list in sx] if additional: - tmp.append(compute_additional_sx_cache(n, is_singlet, qed, sx[:, 0], sx[:, -1])) + tmp.append( + compute_additional_sx_cache( + n, + is_singlet, + qed, + sx[:, 0], + sx[:, -1], + sx[2, 1 if max_weight == 3 else 3], + ) + ) return tmp @nb.njit(cache=True) -def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None): +def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None, sm21=None): r"""Get the harmonics sums cache needed for the qed non singlet AD. We have three cases : @@ -240,12 +250,16 @@ def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None): ---------- n : complex Mellin moment - s1 : float - harmonic sum :math:`S_1(N)` qcd_singlet : bool True if singlet for pure QCD evolution qed : bool True if QED evolution activated + s1 : list + harmonic sum :math:`S_i(N)` + s1 : list + harmonic sum :math:`S_{-i}(N)` + sm21 : complex + harmonic sum :math:`S_{-2,1}(N)` Returns ------- @@ -267,47 +281,79 @@ def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None): """ if sminus is None: + sm1 = None # needed later s1h = S1(n / 2.0) - sx_new = [s1h] s2h = S2(n / 2.0) - sx_new.append(s2h) s3h = S3(n / 2.0) - sx_new.append(s3h) g3N = g_functions.mellin_g3(n, sx[0]) - sx_new.append(g3N) + sx_new = [s1h, s2h, s3h, g3N] if not is_singlet: s1_nm1_h = S1((n - 1.0) / 2.0) - sx_new.append(s1_nm1_h) s2_nm1_h = S2((n - 1.0) / 2.0) - sx_new.append(s2_nm1_h) s3_nm1_h = S3((n - 1.0) / 2.0) - sx_new.append(s3_nm1_h) + sx_new.extend([s1_nm1_h, s2_nm1_h, s3_nm1_h]) else: + sm1 = sminus[0] # needed later if is_singlet: s1h = sminus[0] + sx[0] - sx_new = [s1h] s2h = 2 * (sminus[1] + sx[1]) - sx_new.append(s2h) s3h = 4 * (sminus[2] + sx[2]) - sx_new.append(s3h) - g3N = g_functions.mellin_g3(n, sx[0]) - sx_new.append(g3N) + g3N = -( + zeta2 * (sminus[0] - 1 / n) + - 5.0 / 8 * zeta3 + + zeta2 * log2 + - (sm21 - 1 / n**2 * sx[0]) + ) + sx_new = [s1h, s2h, s3h, g3N] else: - sx_new = compute_additional_sx_cache(n, True, False, sx, None) + s1h = S1(n / 2.0) + s2h = S2(n / 2.0) + s3h = S3(n / 2.0) + g3N = ( + zeta2 * (sminus[0] + 1 / n) + - 5.0 / 8 * zeta3 + + zeta2 * log2 + - (sm21 + 1 / n**2 * sx[0]) + ) s1_nm1_h = sminus[0] + sx[0] - sx_new.append(s1_nm1_h) s2_nm1_h = 2 * (sminus[1] + sx[1]) - sx_new.append(s2_nm1_h) s3_nm1_h = 4 * (sminus[2] + sx[2]) - sx_new.append(s3_nm1_h) + sx_new = [s1h, s2h, s3h, g3N, s1_nm1_h, s2_nm1_h, s3_nm1_h] if qed: - S1p2 = polygamma.recursive_harmonic_sum(sx[0], n, 2, 1) - g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2) - sx_new.append(g3Np2) + g3Np2 = g3np2(n, is_singlet, sx[0], sm21, sm1) s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1.0) / 2.0, 1, 1) - sx_new.append(s1_np1_h) s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1.0) / 2.0, 1, 2) - sx_new.append(s2_np1_h) s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1.0) / 2.0, 1, 3) - sx_new.append(s3_np1_h) + sx_new.extend([g3Np2, s1_np1_h, s2_np1_h, s3_np1_h]) return sx_new + + +@nb.njit(cache=True) +def g3np2(n, is_singlet, s1, sm21=None, sm1=None): + """ + Compute g3(n+2). + + Parameters + ---------- + n : complex + Mellin moment + is_singlet : bool + True if singlet for pure QCD evolution + qed : bool + True if QED evolution activated + s1 : complex + harmonic sum :math:`S_1(N)` + + Returns + ------- + g3(n+2) : complex + """ + eta = 1 if is_singlet else -1 + if sm21 is None: + s1_np2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1) + return g_functions.mellin_g3(n + 2.0, s1_np2) + else: + s1_np1 = polygamma.recursive_harmonic_sum(s1, n, 1, 1) + sm1_np1 = sm1 - eta / (n + 1) + sm21_np1 = sm21 - eta / (n + 1) ** 2 * s1_np1 + return -1 / eta * (zeta2 * sm1_np1 - 5.0 / 8 * zeta3 + zeta2 * log2 - sm21_np1) From 4a583abbc670406eb43d31c265b15cc6fad0357a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 10 Jan 2023 12:04:31 +0100 Subject: [PATCH 12/14] Simplify test_init.py --- tests/eko/evolution_operator/test_init.py | 43 +++++++++++------------ 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index db9174c91..92ac5ff51 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -287,14 +287,23 @@ def test_quad_ker(monkeypatch): np.testing.assert_allclose(res_ns, 0.0) +class FakeCoupling: + def __init__(self): + self.alphaem_running = None + self.q2_ref = 0.0 + + def a(self, scale_to=None, fact_scale=None, nf_to=None): + return (0.1, 0.01) + + def compute_aem_as(self, aem_ref, as_from, as_to, nf): + return aem_ref + + +fake_managers = {"couplings": FakeCoupling()} + + class TestOperator: def test_labels(self, theory_ffns, operator_card, tmp_path): - tcard: TheoryCard = theory_ffns(3) - tcard.xif = 2.0 - ocard: OperatorCard = operator_card - ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED - r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar") - g = r.op_grid o = Operator( dict( order=(3, 0), @@ -303,7 +312,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path): n_integration_cores=1, ModSV=None, ), - g.managers, + fake_managers, 3, 1, 2, @@ -317,7 +326,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path): n_integration_cores=1, ModSV=None, ), - g.managers, + fake_managers, 3, 1, 2, @@ -325,12 +334,6 @@ def test_labels(self, theory_ffns, operator_card, tmp_path): assert sorted(o.labels) == [] def test_labels_qed(self, theory_ffns, operator_card, tmp_path): - tcard: TheoryCard = theory_ffns(3) - tcard.xif = 2.0 - ocard: OperatorCard = operator_card - ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED - r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar") - g = r.op_grid o = Operator( dict( order=(3, 1), @@ -340,7 +343,7 @@ def test_labels_qed(self, theory_ffns, operator_card, tmp_path): ModSV=None, ev_op_iterations=1, ), - g.managers, + fake_managers, 3, 1, 2, @@ -355,7 +358,7 @@ def test_labels_qed(self, theory_ffns, operator_card, tmp_path): ModSV=None, ev_op_iterations=1, ), - g.managers, + fake_managers, 3, 1, 2, @@ -367,12 +370,6 @@ def test_n_pools(self, theory_ffns, operator_card, tmp_path): # make sure we actually have more the those cores (e.g. on github we don't) if os.cpu_count() <= excluded_cores: return - tcard: TheoryCard = theory_ffns(3) - tcard.xif = 2.0 - ocard: OperatorCard = operator_card - ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED - r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar") - g = r.op_grid o = Operator( dict( order=(2, 0), @@ -381,7 +378,7 @@ def test_n_pools(self, theory_ffns, operator_card, tmp_path): n_integration_cores=-excluded_cores, ModSV=None, ), - g.managers, + fake_managers, 3, 1, 10, From 8449d4571371871dea01e6afcdd86308967d5fc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 10 Jan 2023 15:16:48 +0100 Subject: [PATCH 13/14] Refactor tests in test_init.py --- tests/eko/evolution_operator/test_init.py | 211 ++++------------------ 1 file changed, 34 insertions(+), 177 deletions(-) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index 92ac5ff51..ba2cb90ba 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -28,183 +28,40 @@ def test_quad_ker(monkeypatch): monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0) monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0) monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2)) - for is_log in [True, False]: - res_ns = quad_ker( - u=0, - order=(1, 0), - mode0=br.non_singlet_pids_map["ns+"], - mode1=0, - method="", - is_log=is_log, - logx=0.0, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_ns, 0.0) - res_ns = quad_ker( - u=0, - order=(3, 1), - mode0=br.non_singlet_pids_map["ns+u"], - mode1=0, - method="", - is_log=is_log, - logx=0.0, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_ns, 0.0) - res_s = quad_ker( - u=0, - order=(1, 0), - mode0=100, - mode1=100, - method="", - is_log=is_log, - logx=0.123, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_s, 1.0) - res_s = quad_ker( - u=0, - order=(1, 1), - mode0=100, - mode1=100, - method="iterate-exact", - is_log=is_log, - logx=0.123, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_s, 1.0) - res_s = quad_ker( - u=0, - order=(1, 0), - mode0=100, - mode1=21, - method="", - is_log=is_log, - logx=0.0, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_s, 0.0) - res_s = quad_ker( - u=0, - order=(1, 1), - mode0=100, - mode1=21, - method="iterate-exact", - is_log=is_log, - logx=0.0, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_s, 0.0) - res_v = quad_ker( - u=0, - order=(1, 1), - mode0=10200, - mode1=10200, - method="iterate-exact", - is_log=is_log, - logx=0.123, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_v, 1.0) - res_v = quad_ker( - u=0, - order=(1, 1), - mode0=10200, - mode1=10204, - method="iterate-exact", - is_log=is_log, - logx=0.123, - areas=np.zeros(3), - as1=1, - as0=2, - as_raw=1, - aem_list=[0.00058], - alphaem_running=False, - nf=3, - L=0, - ev_op_iterations=0, - ev_op_max_order=(0, 0), - sv_mode=1, - is_threshold=False, - ) - np.testing.assert_allclose(res_v, 0.0) + params = [ + ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0), + ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0), + ((1, 0), 100, 100, "", 0.123, 1.0), + ((1, 0), 100, 21, "", 0.0, 0.0), + ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0), + ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0), + ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0), + ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0), + ] + for p in params: + for is_log in [True, False]: + res_ns = quad_ker( + u=0, + order=p[0], + mode0=p[1], + mode1=p[2], + method=p[3], + is_log=is_log, + logx=p[4], + areas=np.zeros(3), + as1=1, + as0=2, + as_raw=1, + aem_list=[0.00058], + alphaem_running=False, + nf=3, + L=0, + ev_op_iterations=0, + ev_op_max_order=(0, 0), + sv_mode=1, + is_threshold=False, + ) + np.testing.assert_allclose(res_ns, p[5]) for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: for sv in [2, 3]: res_sv = quad_ker( From 018b1216b4c743f2480be61fbe9ef28283b30f1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 10 Jan 2023 15:28:24 +0100 Subject: [PATCH 14/14] REmove unnecessary parameters from test --- tests/eko/evolution_operator/test_init.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index ba2cb90ba..59f4fbacd 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -160,7 +160,7 @@ def compute_aem_as(self, aem_ref, as_from, as_to, nf): class TestOperator: - def test_labels(self, theory_ffns, operator_card, tmp_path): + def test_labels(self): o = Operator( dict( order=(3, 0), @@ -190,7 +190,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path): ) assert sorted(o.labels) == [] - def test_labels_qed(self, theory_ffns, operator_card, tmp_path): + def test_labels_qed(self): o = Operator( dict( order=(3, 1), @@ -222,7 +222,7 @@ def test_labels_qed(self, theory_ffns, operator_card, tmp_path): ) assert sorted(o.labels) == [] - def test_n_pools(self, theory_ffns, operator_card, tmp_path): + def test_n_pools(self): excluded_cores = 3 # make sure we actually have more the those cores (e.g. on github we don't) if os.cpu_count() <= excluded_cores: