From 975c0971b7f4bc21f55e4797cabd1fe8e960d292 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 9 Jan 2023 15:10:35 +0100 Subject: [PATCH 01/17] Add minus signs --- src/eko/scale_variations/expanded.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 7eb35750e..bef867ea7 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -1,4 +1,4 @@ -r"""This module contains the scale variation operator for the expanded scheme (``ModSV=expanded``). +r"""Contains the scale variation operator for the expanded scheme (``ModSV=expanded``). The expressions can be obtained using Eqs. (3.33) and (3.38) of :cite:`AbdulKhalek:2019ihb`. Be aware that corresponding the signs of the ingredients there are a number of differences. @@ -14,7 +14,7 @@ @nb.njit(cache=True) def variation_as1(gamma, L): - r"""Computes the |NLO| anomalous dimension variation. + r"""Compute the |NLO| anomalous dimension variation. Parameters ---------- @@ -33,7 +33,7 @@ def variation_as1(gamma, L): @nb.njit(cache=True) def variation_as2(gamma, L, beta0, g0e2): - r"""Computes the |NNLO| anomalous dimension variation. + r"""Compute the |NNLO| anomalous dimension variation. Parameters ---------- @@ -56,7 +56,7 @@ def variation_as2(gamma, L, beta0, g0e2): @nb.njit(cache=True) def variation_as3(gamma, L, beta0, beta1, g0e2, g0e3, g1g0, g0g1): - r"""Computes the |N3LO| anomalous dimension variation. + r"""Compute the |N3LO| anomalous dimension variation. Parameters ---------- @@ -117,14 +117,14 @@ def non_singlet_variation(gamma, a_s, order, nf, L): """ sv_ker = 1.0 if order[0] >= 2: - sv_ker += a_s * variation_as1(gamma, L) + sv_ker -= a_s * variation_as1(gamma, L) if order[0] >= 3: beta0 = beta.beta_qcd_as2(nf) - sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma[0] ** 2) + sv_ker -= a_s**2 * variation_as2(gamma, L, beta0, gamma[0] ** 2) if order[0] >= 4: beta1 = beta.beta_qcd((3, 0), nf) g0g1 = gamma[0] * gamma[1] - sv_ker += a_s**3 * variation_as3( + sv_ker -= a_s**3 * variation_as3( gamma, L, beta0, beta1, gamma[0] ** 2, gamma[0] ** 3, g0g1, g0g1 ) return sv_ker @@ -155,18 +155,18 @@ def singlet_variation(gamma, a_s, order, nf, L): sv_ker = np.eye(2, dtype=np.complex_) gamma = np.ascontiguousarray(gamma) if order[0] >= 2: - sv_ker += a_s * variation_as1(gamma, L) + sv_ker -= a_s * variation_as1(gamma, L) if order[0] >= 3: beta0 = beta.beta_qcd_as2(nf) gamma0e2 = gamma[0] @ gamma[0] - sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma0e2) + sv_ker -= a_s**2 * variation_as2(gamma, L, beta0, gamma0e2) if order[0] >= 4: beta1 = beta.beta_qcd((3, 0), nf) gamma0e3 = gamma0e2 @ gamma[0] # here the product is not commutative g1g0 = gamma[1] @ gamma[0] g0g1 = gamma[0] @ gamma[1] - sv_ker += a_s**3 * variation_as3( + sv_ker -= a_s**3 * variation_as3( gamma, L, beta0, beta1, gamma0e2, gamma0e3, g1g0, g0g1 ) return sv_ker From 20008591d5a413b270f4d9a5cf890a4d5608e09f Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 9 Jan 2023 15:25:50 +0100 Subject: [PATCH 02/17] Remove a_raw --- src/eko/evolution_operator/__init__.py | 13 +++---------- tests/eko/evolution_operator/test_init.py | 8 +------- 2 files changed, 4 insertions(+), 17 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 323b117b4..74111859e 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -123,7 +123,6 @@ def quad_ker( areas, as1, as0, - as_raw, nf, L, ev_op_iterations, @@ -155,8 +154,6 @@ def quad_ker( target coupling value as0 : float initial coupling value - as_raw : float - coupling value at the process scale nf : int number of active flavors L : float @@ -201,7 +198,7 @@ def quad_ker( # scale var expanded is applied on the kernel if sv_mode == sv.Modes.expanded and not is_threshold: ker = np.ascontiguousarray( - sv.expanded.singlet_variation(gamma_singlet, as_raw, order, nf, L) + sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L) ) @ np.ascontiguousarray(ker) ker = select_singlet_element(ker, mode0, mode1) else: @@ -218,9 +215,7 @@ def quad_ker( ev_op_iterations, ) if sv_mode == sv.Modes.expanded and not is_threshold: - ker = ( - sv.expanded.non_singlet_variation(gamma_ns, as_raw, order, nf, L) * ker - ) + ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker # recombine everything return np.real(ker * integrand) @@ -316,8 +311,7 @@ def a_s(self): self.mur2_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf ) a1 = sc.a_s(self.mur2_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf) - a_raw = sc.a_s(self.q2_to, fact_scale=self.q2_to, nf_to=self.nf) - return (a0, a1, a_raw) + return (a0, a1) @property def labels(self): @@ -375,7 +369,6 @@ def quad_ker(self, label, logx, areas): areas=areas, as1=self.a_s[1], as0=self.a_s[0], - as_raw=self.a_s[2], nf=self.nf, L=np.log(self.xif2), ev_op_iterations=self.config["ev_op_iterations"], diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index dc8ddbd8e..1fc5a0fcd 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -37,7 +37,6 @@ def test_quad_ker(monkeypatch): areas=np.zeros(3), as1=1, as0=2, - as_raw=1, nf=3, L=0, ev_op_iterations=0, @@ -57,7 +56,6 @@ def test_quad_ker(monkeypatch): areas=np.zeros(3), as1=1, as0=2, - as_raw=1, nf=3, L=0, ev_op_iterations=0, @@ -77,7 +75,6 @@ def test_quad_ker(monkeypatch): areas=np.zeros(3), as1=1, as0=2, - as_raw=1, nf=3, L=0, ev_op_iterations=0, @@ -99,7 +96,6 @@ def test_quad_ker(monkeypatch): areas=np.zeros(3), as1=1, as0=2, - as_raw=1, nf=3, L=0, ev_op_iterations=0, @@ -121,7 +117,6 @@ def test_quad_ker(monkeypatch): areas=np.zeros(3), as1=1, as0=2, - as_raw=1, nf=3, L=0, ev_op_iterations=0, @@ -314,7 +309,7 @@ def quad_ker_pegasus( mode1 = 0 method = "" logxs = np.log(int_disp.xgrid.raw) - as_raw = a1 = 1 + a1 = 1 a0 = 2 nf = 3 L = 0 @@ -335,7 +330,6 @@ def quad_ker_pegasus( bf.areas_representation, a1, a0, - as_raw, nf, L, ev_op_iterations, From 74938d4099b9c7ab0ee979618d8adc89552aa084 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 9 Jan 2023 17:50:20 +0100 Subject: [PATCH 03/17] Fix signs and tests --- src/eko/scale_variations/expanded.py | 3 ++- tests/eko/scale_variations/test_expanded.py | 11 ++++++++--- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index bef867ea7..4faedc7bf 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -51,7 +51,8 @@ def variation_as2(gamma, L, beta0, g0e2): complex variation at |NNLO| """ - return -gamma[1] * L + 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 + return -gamma[1] * L - 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 + # changing alpha will change the sign of beta*gamma0 @nb.njit(cache=True) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index bce253f9c..d915f3db0 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -64,12 +64,17 @@ def scheme_diff(g, k, pto, is_singlet): the accuracy of singlet quantities is slightly worse. """ if pto[0] >= 2: - diff = g[0] * k * a0 + 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 += a0**2 * g[1] * k - k**2 * ( - 1 / 2 * a0**2 * b0 * g[0] + a1 * a0 * g02 - 1 / 2 * a0**2 * g02 + 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 ) return diff From b4bfd1f6626c0e44ec8c5cad4cddb941ead8f340 Mon Sep 17 00:00:00 2001 From: Andrea Barontini Date: Mon, 9 Jan 2023 17:57:44 +0100 Subject: [PATCH 04/17] Update src/eko/scale_variations/expanded.py Co-authored-by: Alessandro Candido --- src/eko/scale_variations/expanded.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 4faedc7bf..00185e3db 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -1,4 +1,4 @@ -r"""Contains the scale variation operator for the expanded scheme (``ModSV=expanded``). +r"""Scale variation operator for the expanded scheme (``ModSV=expanded``). The expressions can be obtained using Eqs. (3.33) and (3.38) of :cite:`AbdulKhalek:2019ihb`. Be aware that corresponding the signs of the ingredients there are a number of differences. From e790fd443d8ce3ec701e8d2c553b7de45df5528e Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 11 Jan 2023 14:31:30 +0100 Subject: [PATCH 05/17] Fix test --- tests/eko/scale_variations/test_expanded.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index d915f3db0..170004d04 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -68,10 +68,20 @@ def scheme_diff(g, k, pto, is_singlet): if pto[0] >= 3: b0 = beta_qcd_as2(nf) g02 = g[0] @ g[0] if is_singlet else g[0] ** 2 + # This would be the difference if scheme B was expanded in terms + # of alpha(Q). However, it is expanded in terms of alpha(xi_F*Q) + # so one of the terms gets canceled by the change of alpha. + # 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 + # ) 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 From 462a80a67bc498c2d4f3d7576db9812d290a7e68 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Wed, 11 Jan 2023 14:53:39 +0100 Subject: [PATCH 06/17] Make tests more verbose --- tests/eko/scale_variations/test_expanded.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 170004d04..c1fb72bc6 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -104,7 +104,12 @@ def scheme_diff(g, k, pto, is_singlet): ) ker_b = ker * expanded.non_singlet_variation(gns, a1, order, nf, L) ns_diff = scheme_diff(gns, L, order, False) - np.testing.assert_allclose((ker_a - ker_b) / ker, ns_diff, atol=1e-3) + np.testing.assert_allclose( + (ker_a - ker_b) / ker, + ns_diff, + atol=1e-3, + err_msg=f"L={L},order={order},non-singlet", + ) # Singlet kernels gs = gamma_singlet(order, n, nf) @@ -118,5 +123,8 @@ def scheme_diff(g, k, pto, is_singlet): ker_b = ker @ expanded.singlet_variation(gs, a1, order, nf, L) s_diff = scheme_diff(gs, L, order, True) np.testing.assert_allclose( - (ker_a - ker_b) @ np.linalg.inv(ker), s_diff, atol=5e-3 + (ker_a - ker_b) @ np.linalg.inv(ker), + s_diff, + atol=5e-3, + err_msg=f"L={L},order={order},singlet", ) From 1d233239106b83dfc71095254ac7f8344b994f3b Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 11 Jan 2023 17:08:40 +0100 Subject: [PATCH 07/17] first attempt on N3LO expanded test --- tests/eko/scale_variations/test_expanded.py | 47 ++++++++++++++------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index c1fb72bc6..1b850e489 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,7 @@ from eko import basis_rotation as br from eko.anomalous_dimensions import gamma_ns, gamma_singlet -from eko.beta import beta_qcd_as2 +from eko.beta import beta_qcd_as2, beta_qcd_as3 from eko.kernels import non_singlet, singlet from eko.scale_variations import Modes, expanded, exponentiated @@ -63,22 +63,15 @@ def scheme_diff(g, k, pto, is_singlet): Effects due to non commutativity are neglected thus, the accuracy of singlet quantities is slightly worse. """ + # Note that scheme B is expended in terms of alpha(xi_F*Q) + # so once you have computed the difference A -B + # you need to subtract some terms coming from the variation of + # alpha 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 - # This would be the difference if scheme B was expanded in terms - # of alpha(Q). However, it is expanded in terms of alpha(xi_F*Q) - # so one of the terms gets canceled by the change of alpha. - # 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 - # ) diff += ( -2 * a1**2 * g[1] * k + a0**2 * g[1] * k @@ -86,13 +79,35 @@ def scheme_diff(g, k, pto, is_singlet): - a1 * a0 * g02 * k**2 + 0.5 * a0**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] + # TODO: double check this + 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 + + 1 / 3 * a0**3 * b0**2 * g[0] * k**3 + - 5 / 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 + + 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 + ) return diff - # TODO: perform this test also at N3LO, once evolution kernels - # will be implemented - for L in [np.log(0.5), np.log(2)]: - for order in [(2, 0), (3, 0)]: + for order in [(2, 0), (3, 0), (4, 0)]: # Non singlet kernels gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf) ker = non_singlet.dispatcher( From 45c266a6c47098cb3a3388f181b3618b003b5b6e Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 11 Jan 2023 17:17:09 +0100 Subject: [PATCH 08/17] Revert "first attempt on N3LO expanded test" This reverts commit 1d233239106b83dfc71095254ac7f8344b994f3b. --- tests/eko/scale_variations/test_expanded.py | 47 +++++++-------------- 1 file changed, 16 insertions(+), 31 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 1b850e489..c1fb72bc6 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,7 @@ from eko import basis_rotation as br from eko.anomalous_dimensions import gamma_ns, gamma_singlet -from eko.beta import beta_qcd_as2, beta_qcd_as3 +from eko.beta import beta_qcd_as2 from eko.kernels import non_singlet, singlet from eko.scale_variations import Modes, expanded, exponentiated @@ -63,15 +63,22 @@ def scheme_diff(g, k, pto, is_singlet): Effects due to non commutativity are neglected thus, the accuracy of singlet quantities is slightly worse. """ - # Note that scheme B is expended in terms of alpha(xi_F*Q) - # so once you have computed the difference A -B - # you need to subtract some terms coming from the variation of - # alpha 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 + # This would be the difference if scheme B was expanded in terms + # of alpha(Q). However, it is expanded in terms of alpha(xi_F*Q) + # so one of the terms gets canceled by the change of alpha. + # 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 + # ) diff += ( -2 * a1**2 * g[1] * k + a0**2 * g[1] * k @@ -79,35 +86,13 @@ def scheme_diff(g, k, pto, is_singlet): - a1 * a0 * g02 * k**2 + 0.5 * a0**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] - # TODO: double check this - 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 - + 1 / 3 * a0**3 * b0**2 * g[0] * k**3 - - 5 / 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 - + 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 - ) return diff + # TODO: perform this test also at N3LO, once evolution kernels + # will be implemented + for L in [np.log(0.5), np.log(2)]: - for order in [(2, 0), (3, 0), (4, 0)]: + for order in [(2, 0), (3, 0)]: # Non singlet kernels gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf) ker = non_singlet.dispatcher( From dfd5e922e73a90919e0980becf0e6af82422606c Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Thu, 12 Jan 2023 17:52:25 +0100 Subject: [PATCH 09/17] Fix sign of NNLO kernel scheme B --- src/eko/scale_variations/expanded.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 00185e3db..5ad651f04 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -51,7 +51,7 @@ def variation_as2(gamma, L, beta0, g0e2): complex variation at |NNLO| """ - return -gamma[1] * L - 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 + return -gamma[1] * L + 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 # changing alpha will change the sign of beta*gamma0 From 910c5b4e0faab6b6f27253429ddfdc81513d5d23 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 19 Jan 2023 10:05:38 +0100 Subject: [PATCH 10/17] fix test sv expended inconsitently --- tests/eko/scale_variations/test_expanded.py | 58 +++++++++++++-------- 1 file changed, 36 insertions(+), 22 deletions(-) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index c1fb72bc6..84bde35a8 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,7 @@ from eko import basis_rotation as br from eko.anomalous_dimensions import gamma_ns, gamma_singlet -from eko.beta import beta_qcd_as2 +from eko.beta import beta_qcd_as2, beta_qcd_as3 from eko.kernels import non_singlet, singlet from eko.scale_variations import Modes, expanded, exponentiated @@ -43,13 +43,12 @@ def test_singlet_sv_dispacher(): 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}` depends - only on the accuracy in the :math:`\alpha_s` expansion - and not in the size of the `fact_to_ren` itself. - However in our implementation the exponentiated mode depends on - the actual value of a0 and a1, since the evolution integral in :math:`\alpha_s` - is evaluated. Thus this test ratio :math:`(ker_A - ker_B)/ker_{unv}` - still contains a dependency on `fact_to_ren`. + 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. """ nf = 5 n = 10 @@ -68,31 +67,46 @@ def scheme_diff(g, k, pto, is_singlet): if pto[0] >= 3: b0 = beta_qcd_as2(nf) g02 = g[0] @ g[0] if is_singlet else g[0] ** 2 - # This would be the difference if scheme B was expanded in terms - # of alpha(Q). However, it is expanded in terms of alpha(xi_F*Q) - # so one of the terms gets canceled by the change of alpha. - # 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 - # ) 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 ) return diff - # TODO: perform this test also at N3LO, once evolution kernels - # will be implemented for L in [np.log(0.5), np.log(2)]: - for order in [(2, 0), (3, 0)]: + for order in [(2, 0), (3, 0), (4,0)]: # Non singlet kernels gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf) ker = non_singlet.dispatcher( From b9d83711ab4770636c8c5903ec159958dadd710e Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 19 Jan 2023 10:05:59 +0100 Subject: [PATCH 11/17] fix test sv expended inconsitently --- 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 84bde35a8..c1ecc33fd 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -47,7 +47,7 @@ def test_scale_variation_a_vs_b(): 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-. + shifts in path length and :math:`\alpha_s` values. The real difference is always higher order. """ nf = 5 From 9dcf7e6f6bb71e49ae49cc42a60ebd2522d2b9eb Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 19 Jan 2023 10:11:00 +0100 Subject: [PATCH 12/17] fix name shift --- src/eko/evolution_operator/__init__.py | 10 +++++----- src/eko/matching_conditions/operator_matrix_element.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 74111859e..874ae25ac 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -286,7 +286,7 @@ def grid_size(self): """Return the grid size.""" return self.int_disp.xgrid.size - def mur2_shift(self, q2): + def scheme_exponentiated_shift(self, q2): """Compute shifted renormalization scale. Parameters @@ -308,9 +308,9 @@ def a_s(self): """Return the computed values for :math:`a_s`.""" sc = self.managers["strong_coupling"] a0 = sc.a_s( - self.mur2_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf + self.scheme_exponentiated_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf ) - a1 = sc.a_s(self.mur2_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf) + a1 = sc.a_s(self.scheme_exponentiated_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf) return (a0, a1) @property @@ -466,8 +466,8 @@ def compute(self): logger.info( "%s: µ_R^2 distance: %e -> %e", self.log_label, - self.mur2_shift(self.q2_from), - self.mur2_shift(self.q2_to), + self.scheme_exponentiated_shift(self.q2_from), + self.scheme_exponentiated_shift(self.q2_to), ) if self.sv_mode != sv.Modes.unvaried: logger.info( diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index a8ad60109..4ca0a6ddb 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -366,7 +366,7 @@ def a_s(self): Note that here you need to use :math:`a_s^{n_f+1}` """ sc = self.managers["strong_coupling"] - return sc.a_s(self.mur2_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1) + return sc.a_s(self.scheme_exponentiated_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1) def compute(self): """Compute the actual operators (i.e. run the integrations)""" From bd02ec23efbf38f395987505058d7ae315adbc71 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 19 Jan 2023 10:18:47 +0100 Subject: [PATCH 13/17] fix new func name in test --- tests/eko/evolution_operator/test_init.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index 1fc5a0fcd..294abd7a6 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -184,7 +184,7 @@ def test_exponentiated(self, theory_ffns, operator_card, tmp_path): g = r.op_grid # setup objs o = Operator(g.config, g.managers, 3, 2.0, 10.0) - np.testing.assert_allclose(o.mur2_shift(40.0), 10.0) + np.testing.assert_allclose(o.scheme_exponentiated_shift(40.0), 10.0) o.compute() self.check_lo(o) From 1336f45f22fa3605eb43cbde58fef4e5c8f6db81 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Thu, 19 Jan 2023 12:05:26 +0100 Subject: [PATCH 14/17] Change name of func --- src/eko/evolution_operator/__init__.py | 16 +++++++++++----- .../operator_matrix_element.py | 6 ++++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 874ae25ac..5818c90fc 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -286,7 +286,7 @@ def grid_size(self): """Return the grid size.""" return self.int_disp.xgrid.size - def scheme_exponentiated_shift(self, q2): + def sv_exponentiated_shift(self, q2): """Compute shifted renormalization scale. Parameters @@ -308,9 +308,15 @@ def a_s(self): """Return the computed values for :math:`a_s`.""" sc = self.managers["strong_coupling"] a0 = sc.a_s( - self.scheme_exponentiated_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf + self.sv_exponentiated_shift(self.q2_from), + fact_scale=self.q2_from, + nf_to=self.nf, + ) + a1 = sc.a_s( + self.sv_exponentiated_shift(self.q2_to), + fact_scale=self.q2_to, + nf_to=self.nf, ) - a1 = sc.a_s(self.scheme_exponentiated_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf) return (a0, a1) @property @@ -466,8 +472,8 @@ def compute(self): logger.info( "%s: µ_R^2 distance: %e -> %e", self.log_label, - self.scheme_exponentiated_shift(self.q2_from), - self.scheme_exponentiated_shift(self.q2_to), + self.sv_exponentiated_shift(self.q2_from), + self.sv_exponentiated_shift(self.q2_to), ) if self.sv_mode != sv.Modes.unvaried: logger.info( diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index 4ca0a6ddb..ed0299294 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -155,7 +155,7 @@ def build_ome(A, matching_order, a_s, backward_method): def quad_ker( u, order, mode0, mode1, is_log, logx, areas, a_s, nf, L, backward_method, is_msbar ): - """ + r""" Raw kernel inside quad Parameters @@ -366,7 +366,9 @@ def a_s(self): Note that here you need to use :math:`a_s^{n_f+1}` """ sc = self.managers["strong_coupling"] - return sc.a_s(self.scheme_exponentiated_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1) + return sc.a_s( + self.sv_exponentiated_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1 + ) def compute(self): """Compute the actual operators (i.e. run the integrations)""" From 5ee0d5149f331238b1e7d794d87f41520f49997d Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Thu, 19 Jan 2023 12:11:37 +0100 Subject: [PATCH 15/17] change name also in test --- tests/eko/evolution_operator/test_init.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index 294abd7a6..d40fbeb39 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -184,7 +184,7 @@ def test_exponentiated(self, theory_ffns, operator_card, tmp_path): g = r.op_grid # setup objs o = Operator(g.config, g.managers, 3, 2.0, 10.0) - np.testing.assert_allclose(o.scheme_exponentiated_shift(40.0), 10.0) + np.testing.assert_allclose(o.sv_exponentiated_shift(40.0), 10.0) o.compute() self.check_lo(o) From 007ba7240cc254780c4c833076dd2f4c61d08394 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 23 Jan 2023 17:28:42 +0100 Subject: [PATCH 16/17] Expand comment --- src/eko/scale_variations/expanded.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 5ad651f04..168de28c0 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -52,7 +52,9 @@ def variation_as2(gamma, L, beta0, g0e2): variation at |NNLO| """ return -gamma[1] * L + 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 - # changing alpha will change the sign of beta*gamma0 + # All this kernels are meant to be used with alpha_s evaluated at the + # factorization scale. If one expresses everything in terms of + # alpha_s evaluated at the process scale, the sign of beta*gamma0 flips. @nb.njit(cache=True) From 2afdbacb14d1559bdab084019ae7ba7d75b237c1 Mon Sep 17 00:00:00 2001 From: andreab1997 Date: Mon, 23 Jan 2023 17:52:37 +0100 Subject: [PATCH 17/17] Move comment to docs --- src/eko/scale_variations/expanded.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 168de28c0..3ceb5c6d4 100644 --- a/src/eko/scale_variations/expanded.py +++ b/src/eko/scale_variations/expanded.py @@ -35,6 +35,10 @@ def variation_as1(gamma, L): def variation_as2(gamma, L, beta0, g0e2): r"""Compute the |NNLO| anomalous dimension variation. + These kernels are meant to be used with alpha_s evaluated at the + factorization scale. If one expresses everything in terms of + alpha_s evaluated at the process scale, the sign of beta*gamma0 flips. + Parameters ---------- gamma : numpy.ndarray @@ -52,9 +56,6 @@ def variation_as2(gamma, L, beta0, g0e2): variation at |NNLO| """ return -gamma[1] * L + 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2 - # All this kernels are meant to be used with alpha_s evaluated at the - # factorization scale. If one expresses everything in terms of - # alpha_s evaluated at the process scale, the sign of beta*gamma0 flips. @nb.njit(cache=True)