Scale Variations Expanded

benchmarks/apfel_bench.py

@@ -38,21 +38,24 @@ def skip_pdfs(_theory):
 
 
 class BenchmarkVFNS(ApfelBenchmark):
-    """Benckmark VFNS"""
+    """Benchmark VFNS"""
 
     vfns_theory = {
         "FNS": "ZM-VFNS",
         "ModEv": [
             "EXA",
-            # "EXP",
-            # "TRN",
+            "EXP",
+            "TRN",
         ],
         "kcThr": 1.0,
         "kbThr": 1.0,
         "ktThr": 1.0,
         "Qref": np.sqrt(2.0),
         "alphas": 0.35,
         "Q0": np.sqrt(2.0),
+        "nfref": 3,
+        "nf0": 3,
+        "mc": 1.51,
     }
     vfns_theory = tolist(vfns_theory)
 
@@ -65,11 +68,20 @@ def benchmark_plain(self, pto):
             cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"]
         )
 
-    def benchmark_sv(self, pto):
+    def benchmark_sv(self, pto, svmode):
         """Scale Variation"""
 
         th = self.vfns_theory.copy()
-        th.update({"PTO": [pto], "XIR": [0.7071067811865475, 1.4142135623730951]})
+        th.update(
+            {
+                "PTO": [pto],
+                "XIR": [1 / np.sqrt(2.0)],
+                "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
+                "nfref": [4],
+            }
+        )
         self.run(
             cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"]
         )
@@ -147,7 +159,7 @@ def benchmark_plain(self, pto):
             cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"]
         )
 
-    def benchmark_sv(self, pto):
+    def benchmark_sv(self, pto, svmode):
         """Scale Variation"""
 
         ts = []
@@ -157,6 +169,8 @@ def benchmark_sv(self, pto):
                 "PTO": [pto],
                 "XIR": [np.sqrt(0.5)],
                 "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
             }
         )
         ts.extend(cartesian_product(th))
@@ -166,6 +180,8 @@ def benchmark_sv(self, pto):
                 "PTO": [pto],
                 "XIR": [np.sqrt(2.0)],
                 "fact_to_ren_scale_ratio": [np.sqrt(0.5)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
             }
         )
         ts.extend(cartesian_product(th))
@@ -177,7 +193,7 @@ def benchmark_sv(self, pto):
     obj = BenchmarkVFNS()
     # obj = BenchmarkFFNS()
 
-    # obj.benchmark_plain(1)
-    # obj.benchmark_sv(1)
+    # obj.benchmark_plain(2)
+    obj.benchmark_sv(2, "exponentiated")
     # obj.benchmark_kthr(2)
-    obj.benchmark_msbar(2)
+    # obj.benchmark_msbar(2)

benchmarks/banana.yaml

@@ -1,5 +1,7 @@
-database_path: data/benchmark.db
-input_tables:
-  - theories
-  - operators
-  - cache
+paths:
+  database: data/benchmark.db
+input:
+  tables:
+    - theories
+    - operators
+    - cache

benchmarks/eko/benchmark_ad.py

@@ -14,7 +14,6 @@ def benchmark_melling_g3_pegasus():
         check_melling_g3_pegasus(N)
 
 
-@pytest.mark.isolated
 def check_melling_g3_pegasus(N):
     S1 = h.harmonic_S1(N)
     N1 = N + 1.0
@@ -40,25 +39,23 @@ def check_melling_g3_pegasus(N):
         - 0.3174 * (zeta2 - S15 / N5) / N5
         + 0.0699 * (zeta2 - S16 / N6) / N6
     )
-
     np.testing.assert_allclose(h.mellin_g3(N), SPMOM)
 
 
 @pytest.mark.isolated
 def benchmark_gamma_ns_1_pegasus():
-    for N in range(10):
+    # remember that singlet has pole at N=1
+    for N in [2, 3, 4, 1 + 1j, 1 - 1j, 2 + 1j, 3 + 1j]:
         for NF in [3, 4, 5]:
-            check_gamma_ns_1_pegasus(N, NF)
+            check_gamma_1_pegasus(N, NF)
 
 
-@pytest.mark.isolated
-def check_gamma_ns_1_pegasus(N, NF):
-    # pylint: disable=line-too-long,too-many-locals,too-many-statements
+def check_gamma_1_pegasus(N, NF):
     # Test against pegasus implementation
     ZETA2 = h.zeta2
     ZETA3 = h.zeta3
 
-    N = np.random.rand(1) + np.random.rand(1) * 1j
+    # N = np.random.rand(1) + np.random.rand(1) * 1j
     S1 = h.harmonic_S1(N)
     S2 = h.harmonic_S2(N)
 

benchmarks/eko/benchmark_evol_to_unity.py

@@ -34,6 +34,7 @@ class BenchmarkBackwardForward:
         "MaxNfPdf": 6,
         "MaxNfAs": 6,
         "HQ": "POLE",
+        "ModSV": None,
     }
     operators_card = {
         "Q2grid": [10],

benchmarks/eko/benchmark_msbar_evolution.py

@@ -1,9 +1,12 @@
 # -*- coding: utf-8 -*-
 """This module benchmarks MSbar mass evolution against APFEL."""
+import copy
+
 import numpy as np
 import pytest
 
 from eko import msbar_masses
+from eko.basis_rotation import quark_names
 from eko.strong_coupling import StrongCoupling
 
 # try to load APFEL - if not available, we'll use the cached values
@@ -40,12 +43,18 @@
 @pytest.mark.isolated
 class BenchmarkMSbar:
     def benchmark_APFEL_msbar_evolution(self):
-        Q2s = np.power([1, 96, 150], 2)
-        alphas_ref = 0.118
-        scale_ref = 91.0**2
-        thresholds_ratios = np.power((1.0, 1.0, 1.0), 2)
-        Q2m = np.power([2.0, 4.5, 175], 2)
-        m2 = np.power((1.4, 4.5, 175), 2)
+        theory = copy.deepcopy(theory_dict)
+        bench_values = dict(zip(np.power([1, 96, 150], 2), [3, 5, 5]))
+        theory.update(
+            {
+                "mc": 1.4,
+                "mb": 4.5,
+                "mt": 175.0,
+                "Qmc": 2.0,
+                "Qmb": 4.5,
+                "Qmt": 175.0,
+            }
+        )
         apfel_vals_dict = {
             "exact": {
                 0: np.array(
@@ -96,27 +105,32 @@ def benchmark_APFEL_msbar_evolution(self):
         }
         # collect my values
         for method in ["exact", "expanded"]:
+            theory["ModEv"] = "EXP" if method == "expanded" else "EXA"
             for order in [0, 1, 2]:
+                theory["PTO"] = order
                 as_VFNS = StrongCoupling(
-                    alphas_ref,
-                    scale_ref,
-                    m2,
-                    thresholds_ratios,
+                    theory["alphas"],
+                    theory["Qref"] ** 2,
+                    msbar_masses.compute(theory),
+                    np.power([theory["kcThr"], theory["kbThr"], theory["ktThr"]], 2),
                     order=order,
                     method=method,
-                    hqm_scheme="MSBAR",
+                    nf_ref=theory["nfref"],
+                    hqm_scheme=theory["HQ"],
                 )
                 my_vals = []
-                for Q2 in Q2s:
+                for Q2, nf_to in bench_values.items():
                     my_masses = []
                     for n in [3, 4, 5]:
                         my_masses.append(
                             msbar_masses.evolve(
-                                m2[n - 3],
-                                Q2m[n - 3],
+                                theory[f"m{quark_names[n]}"] ** 2,
+                                theory[f"Qm{quark_names[n]}"] ** 2,
                                 strong_coupling=as_VFNS,
                                 fact_to_ren=1.0,
                                 q2_to=Q2,
+                                nf_ref=n,
+                                nf_to=nf_to,
                             )
                         )
                     my_vals.append(my_masses)
@@ -128,15 +142,17 @@ def benchmark_APFEL_msbar_evolution(self):
                     apfel.SetTheory("QCD")
                     apfel.SetPerturbativeOrder(order)
                     apfel.SetAlphaEvolution(method)
-                    apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref))
+                    apfel.SetAlphaQCDRef(theory["alphas"], theory["Qref"])
                     apfel.SetVFNS()
-                    apfel.SetMSbarMasses(*np.sqrt(m2))
-                    apfel.SetMassScaleReference(*np.sqrt(Q2m))
+                    apfel.SetMSbarMasses(theory["mc"], theory["mb"], theory["mt"])
+                    apfel.SetMassScaleReference(
+                        theory["Qmc"], theory["Qmb"], theory["Qmt"]
+                    )
                     apfel.SetRenFacRatio(1.0)
                     apfel.InitializeAPFEL()
                     # collect apfel masses
                     apfel_vals_cur = []
-                    for Q2 in Q2s:
+                    for Q2 in bench_values:
                         masses = []
                         for n in [4, 5, 6]:
                             masses.append(apfel.HeavyQuarkMass(n, np.sqrt(Q2)))
@@ -147,7 +163,7 @@ def benchmark_APFEL_msbar_evolution(self):
                     )
                 # check myself to APFEL
                 np.testing.assert_allclose(
-                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3
+                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.3e-3
                 )
 
     def benchmark_APFEL_msbar_solution(self):
@@ -184,10 +200,9 @@ def benchmark_APFEL_msbar_solution(self):
                 )
                 apfel.SetRenFacRatio(theory_dict["fact_to_ren_scale_ratio"])
                 apfel.InitializeAPFEL()
-                apfel.EnableWelcomeMessage(1)
             # check myself to APFEL
             np.testing.assert_allclose(
-                apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5
+                apfel_vals, np.sqrt(np.array(my_masses)), rtol=4e-4
             )
 
     def benchmark_msbar_solution_kthr(self):

benchmarks/eko/benchmark_strong_coupling.py

@@ -787,3 +787,81 @@ def benchmark_lhapdf_zmvfns_lo(self):
         # Max absolute difference: 2.58611282e-06
         # Max relative difference: 0.00013379
         np.testing.assert_allclose(lhapdf_vals, np.array(my_vals), rtol=1.5e-4)
+
+    def benchmark_APFEL_fact_to_ren_lha_settings(self):
+
+        theory_dict = {
+            "alphas": 0.35,
+            "Qref": np.sqrt(2.0),
+            "nfref": 4,
+            "nf0": 3,
+            "MaxNfPdf": 6,
+            "MaxNfAs": 6,
+            "Q0": np.sqrt(2.0),
+            "fact_to_ren_scale_ratio": np.sqrt(2.0),
+            "mc": np.sqrt(2.0),
+            "mb": 4.5,
+            "mt": 175.0,
+            "kcThr": 1.0,
+            "kbThr": 1.0,
+            "ktThr": 1.0,
+            "HQ": "POLE",
+            "Qmc": np.sqrt(2.0),
+            "Qmb": 4.5,
+            "Qmt": 175.0,
+            "PTO": 2,
+            "ModEv": "EXA",
+        }
+        Q2s = [2.0]
+        sc = StrongCoupling.from_dict(theory_dict)
+        fact_to_ren = theory_dict["fact_to_ren_scale_ratio"] ** 2
+        for Q2 in Q2s:
+
+            my_val = sc.a_s(Q2 / fact_to_ren, Q2)
+            path = sc.thresholds.path(Q2 / fact_to_ren)
+            my_val_4 = sc.a_s(Q2 / fact_to_ren, Q2, nf_to=4)
+            path_4 = sc.thresholds.path(Q2 / fact_to_ren, 4)
+            my_val_3 = sc.a_s(Q2 / fact_to_ren, Q2, nf_to=3)
+            path_3 = sc.thresholds.path(Q2 / fact_to_ren, 3)
+
+            # path_4 it's not matched
+            assert len(path_4) == 1
+
+            # path_3 is the same as path backward in nf and in q2.
+            assert len(path_3) == 2
+            assert len(path) == 2
+            assert path_3[1].nf < path_3[0].nf
+            assert path_3[1].q2_from < path_3[0].q2_from
+
+            apfel_val_ref = 0.03478112968976964
+            if use_APFEL:
+                # run apfel
+                apfel.CleanUp()
+                apfel.SetTheory("QCD")
+                apfel.SetPerturbativeOrder(theory_dict["PTO"])
+                apfel.SetAlphaEvolution("exact")
+                apfel.SetAlphaQCDRef(theory_dict["alphas"], theory_dict["Qref"])
+                apfel.SetVFNS()
+                apfel.SetPoleMasses(
+                    theory_dict["mc"],
+                    theory_dict["mb"],
+                    theory_dict["mt"],
+                )
+                apfel.SetMassMatchingScales(
+                    theory_dict["kcThr"],
+                    theory_dict["kbThr"],
+                    theory_dict["ktThr"],
+                )
+                apfel.SetRenFacRatio(1.0 / theory_dict["fact_to_ren_scale_ratio"])
+                apfel.InitializeAPFEL()
+                # collect a_s
+                apfel_val = apfel.AlphaQCD(
+                    np.sqrt(Q2) / theory_dict["fact_to_ren_scale_ratio"]
+                ) / (4.0 * np.pi)
+                # check APFEL cached value
+                np.testing.assert_allclose(apfel_val_ref, apfel_val)
+
+            # check myself to APFEL
+            np.testing.assert_allclose(apfel_val_ref, my_val, rtol=0.03)
+            np.testing.assert_allclose(apfel_val_ref, my_val_4)
+            np.testing.assert_allclose(my_val, my_val_3)