From 4f633d4e507ad3d31dd57080ee6b8750f47e22ee Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 13 Feb 2024 13:25:39 +0100 Subject: [PATCH 01/11] init LO test solution --- tests/eko/scale_variations/test_expanded.py | 165 +++++++++++--------- 1 file changed, 91 insertions(+), 74 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index d85b8926a..6fed95f7d 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,10 @@ from eko import basis_rotation as br from eko.beta import beta_qcd_as2, beta_qcd_as3 +from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo +from eko.io.types import ScaleVariationsMethod as svm from eko.kernels import non_singlet, singlet +from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import Modes, expanded, exponentiated from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet @@ -92,72 +95,64 @@ def test_valence_sv_dispacher_qed(): def test_scale_variation_a_vs_b(): - r""" - Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. - We test that the quantity :math:`(ker_A - ker_B)/ker_{unv}` - - Note this is NOT the real difference between scheme expanded - and exponentiated since here we don't take into account the - shifts in path length and :math:`\alpha_s` values. - The real difference is always higher order. + r"""Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. + + We test that the quantity :math:`(ker_A - ker_B)` for truncated solution, + is always higher order difference. + + For simplicity we do FFNS nf=4. """ nf = 5 n = 10 - a1 = 0.118 / (4 * np.pi) - a0 = 0.2 / (4 * np.pi) - method = "truncated" - - def scheme_diff(g, k, pto, is_singlet): - """ - :math:`(ker_A - ker_B)/ker_{unv}` for truncated expansion - Effects due to non commutativity are neglected thus, - the accuracy of singlet quantities is slightly worse. - """ - if pto[0] >= 2: - diff = -g[0] * k * a0 # - 2 * a1 * k * g[0] - # if pto[0] >= 3: - # b0 = beta_qcd_as2(nf) - # g02 = g[0] @ g[0] if is_singlet else g[0] ** 2 - # diff += ( - # -2 * a1**2 * g[1] * k - # + a0**2 * g[1] * k - # + a1**2 * b0 * g[0] * k**2 - # - 0.5 * a0**2 * b0 * g[0] * k**2 - # - a1 * a0 * g02 * k**2 - # + 0.5 * a0**2 * g02 * k**2 - # + a1**2 * g02 * k**2 - # ) - # if pto[0] >= 4: - # b1 = beta_qcd_as3(nf) - # g0g1 = g[0] @ g[1] if is_singlet else g[0] * g[1] - # g03 = g02 @ g[0] if is_singlet else g02 * g[0] - # diff += ( - # a0**3 * g[2] * k - # - 2 * a1**3 * g[2] * k - # - 1 / 2 * a0**3 * b1 * g[0] * k**2 - # + a1**3 * b1 * g[0] * k**2 - # - a0**3 * b0 * g[1] * k**2 - # + 2 * a1**3 * b0 * g[1] * k**2 - # + a0**3 * g0g1 * k**2 - # - a0**2 * a1 * g0g1 * k**2 - # - a0 * a1**2 * g0g1 * k**2 - # + 2 * a1**3 * g0g1 * k**2 - # + 1 / 3 * a0**3 * b0**2 * g[0] * k**3 - # - 2 / 3 * a1**3 * b0**2 * g[0] * k**3 - # - 1 / 2 * a0**3 * b0 * g02 * k**3 - # + 1 / 2 * a0**2 * a1 * b0 * g02 * k**3 - # + 1 / 2 * a0 * a1**2 * b0 * g02 * k**3 - # - a1**3 * b0 * g02 * k**3 - # + 1 / 6 * a0**3 * g03 * k**3 - # - 1 / 2 * a0**2 * a1 * g03 * k**3 - # + 1 / 2 * a0 * a1**2 * g03 * k**3 - # - 1 / 3 * a1**3 * g03 * k**3 - # ) + q02 = 1.65**2 + q12 = 100**2 + ev_method = "truncated" + + ref = CouplingsInfo( + alphas=0.1181, + alphaem=0.007496, + scale=91.00, + max_num_flavs=6, + num_flavs_ref=4, + ) + + def compute_a_s(q2, nf, order): + sc = Couplings( + couplings=ref, + order=order, + method=CouplingEvolutionMethod.EXPANDED, + masses=np.array([1.51, 4.92, 172.5]) ** 2, + hqm_scheme=QuarkMassScheme.POLE, + # thresholds_ratios=np.array([1.0, 1.0, 1.0]) ** 2 * ( + # xif2 if scvar_method == svm.EXPONENTIATED + # else 1.0 + # ) + # Let's do FFNS nf=4 for simplicity + thresholds_ratios=np.array([0, np.inf, np.inf]), + ) + # the multiplication for xif2 here it's done explicitly above + return sc.a_s(scale_to=q2, nf_to=nf) + + def scheme_diff(g, L, order): + """:math:`(ker_A - ker_B)/ker_{unv}` for truncated expansion.""" + if order == (1, 0): + a0 = compute_a_s(q02, nf, order) + diff = g[0] * L * a0 + # TODO: add higher order expressions return diff - for L in [np.log(0.5), np.log(2)]: + for xif2 in [0.5**2, 2**2]: + L = np.log(xif2) # for order in [(2, 0), (3, 0), (4, 0)]: - for order in [(2, 0)]: + for order in [(1, 0)]: + # compute values of alphas + a0_b = compute_a_s(q02, nf, order) + a1_b = compute_a_s(q12 * xif2, nf, order) + + a0_a = compute_a_s(q02 * xif2, nf, order) + # for FFNS these 2 will coincide + a1_a = a1_b + # Non singlet kernels gns = gamma_ns( order, @@ -166,36 +161,58 @@ def scheme_diff(g, k, pto, is_singlet): nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), ) - ker = non_singlet.dispatcher( - order, method, gns, a1, a0, nf, ev_op_iterations=1 + + # build scheme B solution + ker_b = non_singlet.dispatcher( + order, ev_method, gns, a1_b, a0_b, nf, ev_op_iterations=1 ) + ker_b = ker_b * expanded.non_singlet_variation(gns, a1_b, order, nf, L) + + # build scheme A solution gns_a = exponentiated.gamma_variation(gns.copy(), order, nf, L) ker_a = non_singlet.dispatcher( - order, method, gns_a, a1, a0, nf, ev_op_iterations=1 + order, ev_method, gns_a, a1_a, a0_a, nf, ev_op_iterations=1 ) - ker_b = ker * expanded.non_singlet_variation(gns, a1, order, nf, L) - ns_diff = scheme_diff(gns, L, order, False) + + ns_diff = scheme_diff(gns, L, order) np.testing.assert_allclose( - (ker_a - ker_b) / ker, + (ker_a - ker_b), ns_diff, - atol=1e-3, err_msg=f"L={L},order={order},non-singlet", ) # Singlet kernels gs = gamma_singlet(order, n, nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0)) - ker = singlet.dispatcher( - order, method, gs, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1 + + # build scheme B solution + ker_b = singlet.dispatcher( + order, + ev_method, + gs, + a1_b, + a0_b, + nf, + ev_op_iterations=1, + ev_op_max_order=1, ) + ker_b = ker_b @ expanded.singlet_variation(gs, a1_b, order, nf, L, 2) + + # build scheme A solution gs_a = exponentiated.gamma_variation(gs.copy(), order, nf, L) ker_a = singlet.dispatcher( - order, method, gs_a, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1 + order, + ev_method, + gs_a, + a1_a, + a0_a, + nf, + ev_op_iterations=1, + ev_op_max_order=1, ) - ker_b = ker @ expanded.singlet_variation(gs, a1, order, nf, L, 2) - s_diff = scheme_diff(gs, L, order, True) + + s_diff = scheme_diff(gs, L, order) np.testing.assert_allclose( - (ker_a - ker_b) @ np.linalg.inv(ker), + (ker_a - ker_b), s_diff, - atol=5e-3, err_msg=f"L={L},order={order},singlet", ) From 5725ce07ce9d80b391e603e2635822b6cb023bd0 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 13 Feb 2024 13:55:19 +0100 Subject: [PATCH 02/11] attempt 2 to fix LO test --- tests/eko/scale_variations/test_expanded.py | 28 ++++++++++++--------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 6fed95f7d..83b2d0668 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -97,12 +97,12 @@ def test_valence_sv_dispacher_qed(): def test_scale_variation_a_vs_b(): r"""Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. - We test that the quantity :math:`(ker_A - ker_B)` for truncated solution, + We test that the quantity :math:`ker_A / ker_B` for truncated solution, is always higher order difference. For simplicity we do FFNS nf=4. """ - nf = 5 + nf = 4 n = 10 q02 = 1.65**2 q12 = 100**2 @@ -121,23 +121,26 @@ def compute_a_s(q2, nf, order): couplings=ref, order=order, method=CouplingEvolutionMethod.EXPANDED, - masses=np.array([1.51, 4.92, 172.5]) ** 2, + masses=np.array([1.51, 1e4, 1e5]) ** 2, hqm_scheme=QuarkMassScheme.POLE, # thresholds_ratios=np.array([1.0, 1.0, 1.0]) ** 2 * ( # xif2 if scvar_method == svm.EXPONENTIATED # else 1.0 # ) # Let's do FFNS nf=4 for simplicity - thresholds_ratios=np.array([0, np.inf, np.inf]), + thresholds_ratios=np.array([1.0, 1.0, 1.0]), ) # the multiplication for xif2 here it's done explicitly above return sc.a_s(scale_to=q2, nf_to=nf) - def scheme_diff(g, L, order): - """:math:`(ker_A - ker_B)/ker_{unv}` for truncated expansion.""" + def scheme_diff(g, L, order, is_singlet): + """:math:`ker_A / ker_B` for truncated expansion.""" + + idx = np.eye(2) if is_singlet else 1 if order == (1, 0): a0 = compute_a_s(q02, nf, order) - diff = g[0] * L * a0 + b0 = beta_qcd_as2(nf) + diff = (idx + b0 * L * a0) ** (g[0] / b0) # TODO: add higher order expressions return diff @@ -174,9 +177,9 @@ def scheme_diff(g, L, order): order, ev_method, gns_a, a1_a, a0_a, nf, ev_op_iterations=1 ) - ns_diff = scheme_diff(gns, L, order) + ns_diff = scheme_diff(gns, L, order, False) np.testing.assert_allclose( - (ker_a - ker_b), + ker_a / ker_b, ns_diff, err_msg=f"L={L},order={order},non-singlet", ) @@ -210,9 +213,10 @@ def scheme_diff(g, L, order): ev_op_max_order=1, ) - s_diff = scheme_diff(gs, L, order) + s_diff = scheme_diff(gs, L, order, True) np.testing.assert_allclose( - (ker_a - ker_b), - s_diff, + np.diag(ker_a / ker_b), + np.diag(s_diff), err_msg=f"L={L},order={order},singlet", + rtol=4e-3, ) From 7815a09e1b5b57280b7941b26d69aaa09d6281a0 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 13 Feb 2024 14:51:07 +0100 Subject: [PATCH 03/11] fix LO singlet test --- tests/eko/scale_variations/test_expanded.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 83b2d0668..873ec4b76 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -7,6 +7,7 @@ from eko.kernels import non_singlet, singlet from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import Modes, expanded, exponentiated +from ekore.anomalous_dimensions import exp_matrix_2D from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet @@ -136,11 +137,12 @@ def compute_a_s(q2, nf, order): def scheme_diff(g, L, order, is_singlet): """:math:`ker_A / ker_B` for truncated expansion.""" - idx = np.eye(2) if is_singlet else 1 if order == (1, 0): a0 = compute_a_s(q02, nf, order) b0 = beta_qcd_as2(nf) - diff = (idx + b0 * L * a0) ** (g[0] / b0) + diff = (1.0 + b0 * L * a0) ** (g[0] / b0) + if is_singlet: + diff = exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0] # TODO: add higher order expressions return diff @@ -215,8 +217,7 @@ def scheme_diff(g, L, order, is_singlet): s_diff = scheme_diff(gs, L, order, True) np.testing.assert_allclose( - np.diag(ker_a / ker_b), - np.diag(s_diff), + ker_a @ np.linalg.inv(ker_b), + s_diff, err_msg=f"L={L},order={order},singlet", - rtol=4e-3, ) From acdff1038cda5fc80c214f35ab0cd59d4d0c3ce3 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 14 Feb 2024 08:54:42 +0100 Subject: [PATCH 04/11] expand analytic LO resut --- tests/eko/scale_variations/test_expanded.py | 39 ++++++++++++++++----- 1 file changed, 30 insertions(+), 9 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 873ec4b76..696f5e131 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -137,21 +137,40 @@ def compute_a_s(q2, nf, order): def scheme_diff(g, L, order, is_singlet): """:math:`ker_A / ker_B` for truncated expansion.""" + b0 = beta_qcd_as2(nf) if order == (1, 0): - a0 = compute_a_s(q02, nf, order) - b0 = beta_qcd_as2(nf) - diff = (1.0 + b0 * L * a0) ** (g[0] / b0) + # series of (1.0 + b0 * L * a0) ** g[0] / b0), L->0 + diff = ( + 1 + + a0 * g[0] * L + + 1 / 2 * a0**2 * g[0] * (-b0 + g[0]) * L**2 + + 1 / 6 * a0**3 * g[0] * (2 * b0**2 - 3 * b0 * g[0] + g[0] ** 2) * L**3 + ) if is_singlet: - diff = exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0] + # series of exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0] + diff = ( + np.eye(2) + + a0 * g[0] * L + + 1 / 2 * a0**2 * g[0] @ (-b0 + g[0]) * L**2 + + 1 + / 6 + * a0**3 + * g[0] + @ (2 * b0**2 - 3 * b0 * g[0] + g[0] @ g[0]) + * L**3 + ) # TODO: add higher order expressions return diff - for xif2 in [0.5**2, 2**2]: + # let's use smaller scale variation to + # keep the expansions under control + for xif2 in [0.7**2, 1.4**2]: L = np.log(xif2) # for order in [(2, 0), (3, 0), (4, 0)]: for order in [(1, 0)]: # compute values of alphas - a0_b = compute_a_s(q02, nf, order) + a0 = compute_a_s(q02, nf, order) + a0_b = a0 a1_b = compute_a_s(q12 * xif2, nf, order) a0_a = compute_a_s(q02 * xif2, nf, order) @@ -184,6 +203,7 @@ def scheme_diff(g, L, order, is_singlet): ker_a / ker_b, ns_diff, err_msg=f"L={L},order={order},non-singlet", + rtol=2e-5, ) # Singlet kernels @@ -200,7 +220,7 @@ def scheme_diff(g, L, order, is_singlet): ev_op_iterations=1, ev_op_max_order=1, ) - ker_b = ker_b @ expanded.singlet_variation(gs, a1_b, order, nf, L, 2) + ker_b = expanded.singlet_variation(gs, a1_b, order, nf, L, 2) @ ker_b # build scheme A solution gs_a = exponentiated.gamma_variation(gs.copy(), order, nf, L) @@ -217,7 +237,8 @@ def scheme_diff(g, L, order, is_singlet): s_diff = scheme_diff(gs, L, order, True) np.testing.assert_allclose( - ker_a @ np.linalg.inv(ker_b), - s_diff, + np.diag(ker_a @ np.linalg.inv(ker_b)), + np.diag(s_diff), err_msg=f"L={L},order={order},singlet", + rtol=2e-3, ) From cb9e6626e59f16190628bf230a42e81c6c1e2fc6 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 14 Feb 2024 08:58:00 +0100 Subject: [PATCH 05/11] black --- tests/eko/scale_variations/test_expanded.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 696f5e131..8ee34cecd 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -147,7 +147,7 @@ def scheme_diff(g, L, order, is_singlet): + 1 / 6 * a0**3 * g[0] * (2 * b0**2 - 3 * b0 * g[0] + g[0] ** 2) * L**3 ) if is_singlet: - # series of exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0] + # series of exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0], L->0 diff = ( np.eye(2) + a0 * g[0] * L From 265dc94b39f67f4c3369aeb4672b1ad9ec271130 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 14 Feb 2024 11:28:55 +0200 Subject: [PATCH 06/11] Reshuffle SV test --- tests/eko/scale_variations/test_diff.py | 150 +++++++++++++++++++ tests/eko/scale_variations/test_expanded.py | 153 +------------------- 2 files changed, 151 insertions(+), 152 deletions(-) create mode 100644 tests/eko/scale_variations/test_diff.py diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py new file mode 100644 index 000000000..802000e12 --- /dev/null +++ b/tests/eko/scale_variations/test_diff.py @@ -0,0 +1,150 @@ +"""Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. + +We test that the quantity :math:`ker_A / ker_B` for truncated solution, +is always higher order difference. + +For simplicity we do FFNS nf=4. +""" + +import numpy as np + +from eko import basis_rotation as br +from eko.beta import beta_qcd_as2 +from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo +from eko.kernels import non_singlet, singlet +from eko.quantities.heavy_quarks import QuarkMassScheme +from eko.scale_variations import expanded, exponentiated +from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet + +NF = 4 +Q02 = 1.65**2 +Q12 = 100**2 +EV_METHOD = "truncated" + + +def compute_a_s(q2, order): + sc = Couplings( + couplings=CouplingsInfo( + alphas=0.1181, + alphaem=0.007496, + scale=91.00, + max_num_flavs=4, + num_flavs_ref=4, + ), + order=order, + method=CouplingEvolutionMethod.EXPANDED, + masses=np.array([0.0, np.inf, np.inf]), + hqm_scheme=QuarkMassScheme.POLE, + thresholds_ratios=np.array([1.0, 1.0, 1.0]), + ) + # the multiplication for xif2 here it's done explicitly outside + return sc.a_s(scale_to=q2) + + +def scheme_diff_ns(g, a0, L, order): + """:math:`ker_A / ker_B` for truncated non-singlet expansion.""" + + b0 = beta_qcd_as2(NF) + if order == (1, 0): + # series of (1.0 + b0 * L * a0) ** (g[0] / b0), L->0 + diff = 1 + a0 * g[0] * L + 1 / 2 * a0**2 * g[0] * (-b0 + g[0]) * L**2 + # TODO: add higher order expressions + return diff + + +def scheme_diff_s(g, a0, L, order): + """:math:`ker_A / ker_B` for truncated singlet expansion.""" + + b0 = beta_qcd_as2(NF) + if order == (1, 0): + # series of exp(log(1.0 + b0 * L * a0) * g[0] / b0)[0], L->0 + diff = np.eye(2) + a0 * g[0] * L + 1 / 2 * a0**2 * g[0] @ (-b0 + g[0]) * L**2 + # TODO: add higher order expressions + return diff + + +def test_scale_variation_a_vs_b(): + r"""Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``.""" + + # let's use smaller scale variation to + # keep the expansions under control + for xif2 in [0.9, 1.1]: + L = np.log(xif2) + # for order in [(2, 0), (3, 0), (4, 0)]: + for order in [(1, 0)]: + # compute values of alphas + a0 = compute_a_s(Q02, order) + a0_b = a0 + a1_b = compute_a_s(Q12 * xif2, order) + + a0_a = compute_a_s(Q02 * xif2, order) + # for FFNS these 2 will coincide + a1_a = a1_b + for n in [2.0, 3.0, 10.0]: + # Non singlet kernels + gns = gamma_ns( + order, + br.non_singlet_pids_map["ns+"], + n, + NF, + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + ) + + # build scheme B solution + ker_b = non_singlet.dispatcher( + order, EV_METHOD, gns, a1_b, a0_b, NF, ev_op_iterations=1 + ) + ker_b = ker_b * expanded.non_singlet_variation(gns, a1_b, order, NF, L) + + # build scheme A solution + gns_a = exponentiated.gamma_variation(gns.copy(), order, NF, L) + ker_a = non_singlet.dispatcher( + order, EV_METHOD, gns_a, a1_a, a0_a, NF, ev_op_iterations=1 + ) + + ns_diff = scheme_diff_ns(gns, a0, L, order) + np.testing.assert_allclose( + ker_a / ker_b, + ns_diff, + err_msg=f"{L=},{order=},{n=},non-singlet", + rtol=2e-5, + ) + + # Singlet kernels + gs = gamma_singlet( + order, n, NF, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0) + ) + + # build scheme B solution + ker_b = singlet.dispatcher( + order, + EV_METHOD, + gs, + a1_b, + a0_b, + NF, + ev_op_iterations=1, + ev_op_max_order=1, + ) + ker_b = expanded.singlet_variation(gs, a1_b, order, NF, L, 2) @ ker_b + + # build scheme A solution + gs_a = exponentiated.gamma_variation(gs.copy(), order, NF, L) + ker_a = singlet.dispatcher( + order, + EV_METHOD, + gs_a, + a1_a, + a0_a, + NF, + ev_op_iterations=1, + ev_op_max_order=1, + ) + + s_diff = scheme_diff_s(gs, a0, L, order) + np.testing.assert_allclose( + np.diag(ker_a @ np.linalg.inv(ker_b)), + np.diag(s_diff), + err_msg=f"{L=},{order=},{n=},singlet", + rtol=2e-4, + ) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 8ee34cecd..4f21cd314 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -1,13 +1,11 @@ import numpy as np from eko import basis_rotation as br -from eko.beta import beta_qcd_as2, beta_qcd_as3 +from eko.beta import beta_qcd_as2 from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo -from eko.io.types import ScaleVariationsMethod as svm from eko.kernels import non_singlet, singlet from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import Modes, expanded, exponentiated -from ekore.anomalous_dimensions import exp_matrix_2D from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet @@ -93,152 +91,3 @@ def test_valence_sv_dispacher_qed(): ), np.eye(2), ) - - -def test_scale_variation_a_vs_b(): - r"""Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``. - - We test that the quantity :math:`ker_A / ker_B` for truncated solution, - is always higher order difference. - - For simplicity we do FFNS nf=4. - """ - nf = 4 - n = 10 - q02 = 1.65**2 - q12 = 100**2 - ev_method = "truncated" - - ref = CouplingsInfo( - alphas=0.1181, - alphaem=0.007496, - scale=91.00, - max_num_flavs=6, - num_flavs_ref=4, - ) - - def compute_a_s(q2, nf, order): - sc = Couplings( - couplings=ref, - order=order, - method=CouplingEvolutionMethod.EXPANDED, - masses=np.array([1.51, 1e4, 1e5]) ** 2, - hqm_scheme=QuarkMassScheme.POLE, - # thresholds_ratios=np.array([1.0, 1.0, 1.0]) ** 2 * ( - # xif2 if scvar_method == svm.EXPONENTIATED - # else 1.0 - # ) - # Let's do FFNS nf=4 for simplicity - thresholds_ratios=np.array([1.0, 1.0, 1.0]), - ) - # the multiplication for xif2 here it's done explicitly above - return sc.a_s(scale_to=q2, nf_to=nf) - - def scheme_diff(g, L, order, is_singlet): - """:math:`ker_A / ker_B` for truncated expansion.""" - - b0 = beta_qcd_as2(nf) - if order == (1, 0): - # series of (1.0 + b0 * L * a0) ** g[0] / b0), L->0 - diff = ( - 1 - + a0 * g[0] * L - + 1 / 2 * a0**2 * g[0] * (-b0 + g[0]) * L**2 - + 1 / 6 * a0**3 * g[0] * (2 * b0**2 - 3 * b0 * g[0] + g[0] ** 2) * L**3 - ) - if is_singlet: - # series of exp_matrix_2D(np.log(1.0 + b0 * L * a0) * g[0] / b0)[0], L->0 - diff = ( - np.eye(2) - + a0 * g[0] * L - + 1 / 2 * a0**2 * g[0] @ (-b0 + g[0]) * L**2 - + 1 - / 6 - * a0**3 - * g[0] - @ (2 * b0**2 - 3 * b0 * g[0] + g[0] @ g[0]) - * L**3 - ) - # TODO: add higher order expressions - return diff - - # let's use smaller scale variation to - # keep the expansions under control - for xif2 in [0.7**2, 1.4**2]: - L = np.log(xif2) - # for order in [(2, 0), (3, 0), (4, 0)]: - for order in [(1, 0)]: - # compute values of alphas - a0 = compute_a_s(q02, nf, order) - a0_b = a0 - a1_b = compute_a_s(q12 * xif2, nf, order) - - a0_a = compute_a_s(q02 * xif2, nf, order) - # for FFNS these 2 will coincide - a1_a = a1_b - - # Non singlet kernels - gns = gamma_ns( - order, - br.non_singlet_pids_map["ns+"], - n, - nf, - n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), - ) - - # build scheme B solution - ker_b = non_singlet.dispatcher( - order, ev_method, gns, a1_b, a0_b, nf, ev_op_iterations=1 - ) - ker_b = ker_b * expanded.non_singlet_variation(gns, a1_b, order, nf, L) - - # build scheme A solution - gns_a = exponentiated.gamma_variation(gns.copy(), order, nf, L) - ker_a = non_singlet.dispatcher( - order, ev_method, gns_a, a1_a, a0_a, nf, ev_op_iterations=1 - ) - - ns_diff = scheme_diff(gns, L, order, False) - np.testing.assert_allclose( - ker_a / ker_b, - ns_diff, - err_msg=f"L={L},order={order},non-singlet", - rtol=2e-5, - ) - - # Singlet kernels - gs = gamma_singlet(order, n, nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0)) - - # build scheme B solution - ker_b = singlet.dispatcher( - order, - ev_method, - gs, - a1_b, - a0_b, - nf, - ev_op_iterations=1, - ev_op_max_order=1, - ) - ker_b = expanded.singlet_variation(gs, a1_b, order, nf, L, 2) @ ker_b - - # build scheme A solution - gs_a = exponentiated.gamma_variation(gs.copy(), order, nf, L) - ker_a = singlet.dispatcher( - order, - ev_method, - gs_a, - a1_a, - a0_a, - nf, - ev_op_iterations=1, - ev_op_max_order=1, - ) - - s_diff = scheme_diff(gs, L, order, True) - np.testing.assert_allclose( - np.diag(ker_a @ np.linalg.inv(ker_b)), - np.diag(s_diff), - err_msg=f"L={L},order={order},singlet", - rtol=2e-3, - ) From 8d5915cca6c281b40a4ee4b23d54ef9df1cbe461 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 14 Feb 2024 12:59:32 +0200 Subject: [PATCH 07/11] Clean up test/eko/sv/expo --- tests/eko/scale_variations/test_expanded.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 4f21cd314..7be181226 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -1,12 +1,6 @@ import numpy as np -from eko import basis_rotation as br -from eko.beta import beta_qcd_as2 -from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo -from eko.kernels import non_singlet, singlet -from eko.quantities.heavy_quarks import QuarkMassScheme -from eko.scale_variations import Modes, expanded, exponentiated -from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet +from eko.scale_variations import Modes, expanded def test_modes(): From d297c7d1150f198901ab41c68c128cf0604568f4 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 14 Feb 2024 13:37:01 +0200 Subject: [PATCH 08/11] Add test_expanded_is_linear --- tests/eko/scale_variations/test_expanded.py | 108 ++++++++++++++++++++ 1 file changed, 108 insertions(+) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 7be181226..6af22dd51 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -1,6 +1,35 @@ import numpy as np +from eko import basis_rotation as br +from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo +from eko.kernels import non_singlet, singlet +from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import Modes, expanded +from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet + +NF = 4 +Q02 = 1.65**2 +Q12 = 100**2 +EV_METHOD = "truncated" + + +def compute_a_s(q2, order): + sc = Couplings( + couplings=CouplingsInfo( + alphas=0.1181, + alphaem=0.007496, + scale=91.00, + max_num_flavs=4, + num_flavs_ref=4, + ), + order=order, + method=CouplingEvolutionMethod.EXPANDED, + masses=np.array([0.0, np.inf, np.inf]), + hqm_scheme=QuarkMassScheme.POLE, + thresholds_ratios=np.array([1.0, 1.0, 1.0]), + ) + # the multiplication for xif2 here it's done explicitly outside + return sc.a_s(scale_to=q2) def test_modes(): @@ -85,3 +114,82 @@ def test_valence_sv_dispacher_qed(): ), np.eye(2), ) + + +def test_expanded_is_linear(): + r"""Test is linear.""" + for order in [(1, 0), (2, 0), (3, 0), (4, 0)]: + for n in [2.0, 3.0, 10.0]: + rel_err_ns = [] + rel_err_s = [] + for L in [0.3, 0.5, 0.7]: + xif2 = np.exp(L) + # compute values of alphas + a0 = compute_a_s(Q02, order) + a1 = compute_a_s(Q12, order) + a1_b = compute_a_s(Q12 * xif2, order) + # Non singlet kernels + gns = gamma_ns( + order, + br.non_singlet_pids_map["ns+"], + n, + NF, + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + ) + + # build scheme B solution + ker_b = non_singlet.dispatcher( + order, EV_METHOD, gns, a1, a0, NF, ev_op_iterations=1 + ) + sv_b = non_singlet.dispatcher( + order, EV_METHOD, gns, a1_b, a0, NF, ev_op_iterations=1 + ) + sv_b = sv_b * expanded.non_singlet_variation(gns, a1_b, order, NF, L) + + rel_err_ns.append(sv_b / ker_b) + + # Singlet kernels + gs = gamma_singlet( + order, n, NF, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0) + ) + + # build scheme B solution + ker_b = singlet.dispatcher( + order, + EV_METHOD, + gs, + a1, + a0, + NF, + ev_op_iterations=1, + ev_op_max_order=1, + ) + sv_b = singlet.dispatcher( + order, + EV_METHOD, + gs, + a1_b, + a0, + NF, + ev_op_iterations=1, + ev_op_max_order=1, + ) + sv_b = expanded.singlet_variation(gs, a1_b, order, NF, L, 2) @ sv_b + rel_err_s.append(sv_b @ np.linalg.inv(ker_b)) + + # there must be something + for err in rel_err_ns: + assert np.abs(err) != 1.0 + for err in np.array(rel_err_s).flatten(): + assert np.abs(err) != 1.0 + # error has to increase + np.testing.assert_allclose( + rel_err_ns, + sorted(rel_err_ns, reverse=True), + err_msg=f"{order=},{n=},non-singlet", + ) + np.testing.assert_allclose( + rel_err_s, + sorted(rel_err_s, key=np.max, reverse=True), + err_msg=f"{order=},{n=},singlet", + ) From 6118ca0cb5c8c1ea8211256575608db738fd94f1 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 14 Feb 2024 18:13:50 +0100 Subject: [PATCH 09/11] add NLO diff --- tests/eko/scale_variations/test_diff.py | 101 +++++++++++++++++++++--- 1 file changed, 91 insertions(+), 10 deletions(-) diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index 802000e12..b2b9c9aa5 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -9,7 +9,7 @@ import numpy as np from eko import basis_rotation as br -from eko.beta import beta_qcd_as2 +from eko.beta import beta_qcd_as2, beta_qcd_as3 from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo from eko.kernels import non_singlet, singlet from eko.quantities.heavy_quarks import QuarkMassScheme @@ -41,25 +41,105 @@ def compute_a_s(q2, order): return sc.a_s(scale_to=q2) -def scheme_diff_ns(g, a0, L, order): +def scheme_diff_ns(g, a0, a1, L, order): """:math:`ker_A / ker_B` for truncated non-singlet expansion.""" b0 = beta_qcd_as2(NF) + b1 = beta_qcd_as3(NF) if order == (1, 0): # series of (1.0 + b0 * L * a0) ** (g[0] / b0), L->0 diff = 1 + a0 * g[0] * L + 1 / 2 * a0**2 * g[0] * (-b0 + g[0]) * L**2 - # TODO: add higher order expressions + if order == (2, 0): + # this term is formally 1 + as^2 + diff = ( + 1 + - (a1**2 * g[0] * L * (-b1 * g[0] + b0 * (g[1] + b0 * g[0] * L))) / b0**2 + + ( + (a0 * a1 * g[0] * L) + * ( + -2 * b1 * (g[0] + 2 * a1 * b0 * g[0] * L) + + b0 * (3 * b0 * g[0] * L + g[1] * (2 + 4 * a1 * b0 * L)) + ) + ) + / b0**2 + + (a0**2 * L / (2 * b0**5)) + * ( + -2 * a1**2 * b1**3 * g[0] ** 3 + + 6 * a1**2 * b0 * b1**2 * g[0] ** 2 * g[1] + + 3 * b0**6 * g[0] * L + + b0**5 + * ( + -3 * g[0] ** 2 * L + + 2 * a1**2 * g[1] ** 2 * L + + g[1] * (2 - 20 * a1 * g[0] * L) + ) + + (2 * a1 * b0**2 * b1 * g[0]) + * (-3 * a1 * g[1] ** 2 + b1 * g[0] * (-1 + 3 * a1 * g[0] * L)) + + (2 * b0**3) + * ( + a1**2 * g[1] ** 3 + - 5 * a1**2 * b1**2 * g[0] ** 2 * L + + b1 * g[0] * (g[0] + 2 * a1 * g[1] - 5 * a1**2 * g[0] * g[1] * L) + ) + + b0**4 + * ( + -2 * a1 * g[1] ** 2 + + 20 * a1 * b1 * g[0] ** 2 * L + + 2 * g[0] * g[1] * (-1 + 2 * a1**2 * (2 * b1 + g[1]) * L) + ) + ) + ) return diff -def scheme_diff_s(g, a0, L, order): +def scheme_diff_s(g, a0, a1, L, order): """:math:`ker_A / ker_B` for truncated singlet expansion.""" b0 = beta_qcd_as2(NF) + b1 = beta_qcd_as3(NF) if order == (1, 0): # series of exp(log(1.0 + b0 * L * a0) * g[0] / b0)[0], L->0 diff = np.eye(2) + a0 * g[0] * L + 1 / 2 * a0**2 * g[0] @ (-b0 + g[0]) * L**2 - # TODO: add higher order expressions + elif order == (2, 0): + # this term is formally 1 + as^2 + diff = ( + np.eye(2) + - (a1**2 * g[0] * L @ (-b1 * g[0] + b0 * (g[1] + b0 * g[0] * L))) / b0**2 + + ( + (a0 * a1 * g[0] * L) + @ ( + -2 * b1 * (g[0] + 2 * a1 * b0 * g[0] * L) + + b0 * (3 * b0 * g[0] * L + g[1] * (2 + 4 * a1 * b0 * L)) + ) + ) + / b0**2 + + (a0**2 * L / (2 * b0**5)) + * ( + -2 * a1**2 * b1**3 * g[0] @ g[0] @ g[0] + + 6 * a1**2 * b0 * b1**2 * g[0] @ g[0] @ g[1] + + 3 * b0**6 * g[0] * L + + b0**5 + * ( + -3 * g[0] @ g[0] * L + + 2 * a1**2 * g[1] @ g[1] * L + + g[1] @ (2 - 20 * a1 * g[0] * L) + ) + + (2 * a1 * b0**2 * b1 * g[0]) + @ (-3 * a1 * g[1] @ g[1] + b1 * g[0] @ (-1 + 3 * a1 * g[0] * L)) + + (2 * b0**3) + * ( + a1**2 * g[1] @ g[1] @ g[1] + - 5 * a1**2 * b1**2 * g[0] @ g[0] * L + + b1 * g[0] @ (g[0] + 2 * a1 * g[1] - 5 * a1**2 * g[0] @ g[1] * L) + ) + + b0**4 + * ( + -2 * a1 * g[1] @ g[1] + + 20 * a1 * b1 * g[0] @ g[0] * L + + 2 * g[0] @ g[1] @ (-1 + 2 * a1**2 * (2 * b1 + g[1]) * L) + ) + ) + ) return diff @@ -71,9 +151,10 @@ def test_scale_variation_a_vs_b(): for xif2 in [0.9, 1.1]: L = np.log(xif2) # for order in [(2, 0), (3, 0), (4, 0)]: - for order in [(1, 0)]: + for order in [(1, 0), (2, 0)]: # compute values of alphas a0 = compute_a_s(Q02, order) + a1 = compute_a_s(Q12, order) a0_b = a0 a1_b = compute_a_s(Q12 * xif2, order) @@ -102,12 +183,12 @@ def test_scale_variation_a_vs_b(): order, EV_METHOD, gns_a, a1_a, a0_a, NF, ev_op_iterations=1 ) - ns_diff = scheme_diff_ns(gns, a0, L, order) + ns_diff = scheme_diff_ns(gns, a0, a1, L, order) np.testing.assert_allclose( ker_a / ker_b, ns_diff, err_msg=f"{L=},{order=},{n=},non-singlet", - rtol=2e-5, + rtol=2e-5 if order == (1, 0) else 3e-3, ) # Singlet kernels @@ -141,10 +222,10 @@ def test_scale_variation_a_vs_b(): ev_op_max_order=1, ) - s_diff = scheme_diff_s(gs, a0, L, order) + s_diff = scheme_diff_s(gs, a0, a1, L, order) np.testing.assert_allclose( np.diag(ker_a @ np.linalg.inv(ker_b)), np.diag(s_diff), err_msg=f"{L=},{order=},{n=},singlet", - rtol=2e-4, + rtol=2e-4 if order == (1, 0) else 9e-3, ) From b958638f3abed0929c904af962542ab2644bc4c4 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 14 Feb 2024 18:32:21 +0100 Subject: [PATCH 10/11] cleaning --- tests/eko/scale_variations/test_diff.py | 64 ++++--------------------- 1 file changed, 10 insertions(+), 54 deletions(-) diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index b2b9c9aa5..eb66ccd3a 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -56,37 +56,15 @@ def scheme_diff_ns(g, a0, a1, L, order): - (a1**2 * g[0] * L * (-b1 * g[0] + b0 * (g[1] + b0 * g[0] * L))) / b0**2 + ( (a0 * a1 * g[0] * L) - * ( - -2 * b1 * (g[0] + 2 * a1 * b0 * g[0] * L) - + b0 * (3 * b0 * g[0] * L + g[1] * (2 + 4 * a1 * b0 * L)) - ) + * (-2 * b1 * g[0] + b0 * (3 * b0 * g[0] * L + g[1] * 2)) ) / b0**2 + (a0**2 * L / (2 * b0**5)) * ( - -2 * a1**2 * b1**3 * g[0] ** 3 - + 6 * a1**2 * b0 * b1**2 * g[0] ** 2 * g[1] - + 3 * b0**6 * g[0] * L - + b0**5 - * ( - -3 * g[0] ** 2 * L - + 2 * a1**2 * g[1] ** 2 * L - + g[1] * (2 - 20 * a1 * g[0] * L) - ) - + (2 * a1 * b0**2 * b1 * g[0]) - * (-3 * a1 * g[1] ** 2 + b1 * g[0] * (-1 + 3 * a1 * g[0] * L)) - + (2 * b0**3) - * ( - a1**2 * g[1] ** 3 - - 5 * a1**2 * b1**2 * g[0] ** 2 * L - + b1 * g[0] * (g[0] + 2 * a1 * g[1] - 5 * a1**2 * g[0] * g[1] * L) - ) - + b0**4 - * ( - -2 * a1 * g[1] ** 2 - + 20 * a1 * b1 * g[0] ** 2 * L - + 2 * g[0] * g[1] * (-1 + 2 * a1**2 * (2 * b1 + g[1]) * L) - ) + +3 * b0**6 * g[0] * L + + b0**5 * (-3 * g[0] * g[0] * L + g[1] * 2) + + (2 * b0**3) * (+b1 * g[0] * (g[0])) + + b0**4 * (-2 * g[0] * g[1]) ) ) return diff @@ -107,37 +85,15 @@ def scheme_diff_s(g, a0, a1, L, order): - (a1**2 * g[0] * L @ (-b1 * g[0] + b0 * (g[1] + b0 * g[0] * L))) / b0**2 + ( (a0 * a1 * g[0] * L) - @ ( - -2 * b1 * (g[0] + 2 * a1 * b0 * g[0] * L) - + b0 * (3 * b0 * g[0] * L + g[1] * (2 + 4 * a1 * b0 * L)) - ) + @ (-2 * b1 * g[0] + b0 * (3 * b0 * g[0] * L + g[1] * 2)) ) / b0**2 + (a0**2 * L / (2 * b0**5)) * ( - -2 * a1**2 * b1**3 * g[0] @ g[0] @ g[0] - + 6 * a1**2 * b0 * b1**2 * g[0] @ g[0] @ g[1] - + 3 * b0**6 * g[0] * L - + b0**5 - * ( - -3 * g[0] @ g[0] * L - + 2 * a1**2 * g[1] @ g[1] * L - + g[1] @ (2 - 20 * a1 * g[0] * L) - ) - + (2 * a1 * b0**2 * b1 * g[0]) - @ (-3 * a1 * g[1] @ g[1] + b1 * g[0] @ (-1 + 3 * a1 * g[0] * L)) - + (2 * b0**3) - * ( - a1**2 * g[1] @ g[1] @ g[1] - - 5 * a1**2 * b1**2 * g[0] @ g[0] * L - + b1 * g[0] @ (g[0] + 2 * a1 * g[1] - 5 * a1**2 * g[0] @ g[1] * L) - ) - + b0**4 - * ( - -2 * a1 * g[1] @ g[1] - + 20 * a1 * b1 * g[0] @ g[0] * L - + 2 * g[0] @ g[1] @ (-1 + 2 * a1**2 * (2 * b1 + g[1]) * L) - ) + +3 * b0**6 * g[0] * L + + b0**5 * (-3 * g[0] @ g[0] * L + g[1] * 2) + + (2 * b0**3) * (+b1 * g[0] @ (g[0])) + + b0**4 * (-2 * g[0] @ g[1]) ) ) return diff From 02a6855227545fed3686288d746437a692a409e5 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 15 Feb 2024 11:28:57 +0200 Subject: [PATCH 11/11] Apply cosmetic changes to test_diff --- tests/eko/scale_variations/test_diff.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index eb66ccd3a..c0dd67236 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -49,7 +49,7 @@ def scheme_diff_ns(g, a0, a1, L, order): if order == (1, 0): # series of (1.0 + b0 * L * a0) ** (g[0] / b0), L->0 diff = 1 + a0 * g[0] * L + 1 / 2 * a0**2 * g[0] * (-b0 + g[0]) * L**2 - if order == (2, 0): + elif order == (2, 0): # this term is formally 1 + as^2 diff = ( 1 @@ -106,17 +106,14 @@ def test_scale_variation_a_vs_b(): # keep the expansions under control for xif2 in [0.9, 1.1]: L = np.log(xif2) - # for order in [(2, 0), (3, 0), (4, 0)]: for order in [(1, 0), (2, 0)]: # compute values of alphas a0 = compute_a_s(Q02, order) a1 = compute_a_s(Q12, order) a0_b = a0 a1_b = compute_a_s(Q12 * xif2, order) - a0_a = compute_a_s(Q02 * xif2, order) - # for FFNS these 2 will coincide - a1_a = a1_b + a1_a = a1_b # for FFNS these 2 will coincide for n in [2.0, 3.0, 10.0]: # Non singlet kernels gns = gamma_ns(