Fix LHA VFNS SV benchmark 2

benchmarks/apfel_bench.py

@@ -75,7 +75,7 @@ def benchmark_sv(self, pto, svmode):
             {
                 "PTO": [pto],
                 "XIR": [1 / np.sqrt(2.0)],
-                "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "XIF": [np.sqrt(2.0)],
                 "ModSV": [svmode],
                 "EScaleVar": [0],
                 "nfref": [4],
@@ -177,7 +177,7 @@ def benchmark_sv(self, pto, svmode):
             {
                 "PTO": [pto],
                 "XIR": [np.sqrt(0.5)],
-                "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "XIF": [np.sqrt(2.0)],
                 "ModSV": [svmode],
                 "EScaleVar": [0],
             }
@@ -188,7 +188,7 @@ def benchmark_sv(self, pto, svmode):
             {
                 "PTO": [pto],
                 "XIR": [np.sqrt(2.0)],
-                "fact_to_ren_scale_ratio": [np.sqrt(0.5)],
+                "XIF": [np.sqrt(0.5)],
                 "ModSV": [svmode],
                 "EScaleVar": [0],
             }
@@ -198,7 +198,6 @@ def benchmark_sv(self, pto, svmode):
 
 
 if __name__ == "__main__":
-
     # obj = BenchmarkVFNS()
     obj = BenchmarkFFNS()
 

benchmarks/lha_paper_bench.py

@@ -76,12 +76,10 @@ def sv_theories(self, pto):
         """
         low = self.theory.copy()
         low["PTO"] = pto
-        low["fact_to_ren_scale_ratio"] = np.sqrt(1.0 / 2.0)
         low["XIF"] = np.sqrt(1.0 / 2.0)
         low["ModSV"] = "exponentiated"
         high = self.theory.copy()
         high["PTO"] = pto
-        high["fact_to_ren_scale_ratio"] = np.sqrt(2.0)
         high["XIF"] = np.sqrt(2.0)
         high["ModSV"] = "exponentiated"
         return [high, low]
@@ -246,10 +244,11 @@ def run_lha(self, theory_updates):
 
 if __name__ == "__main__":
     # Benchmark to LHA
-    obj = BenchmarkFFNS_polarized()
+    # obj = BenchmarkFFNS_polarized()
     # obj = BenchmarkFFNS()
+    obj = BenchmarkVFNS()
     # obj.benchmark_plain(1)
-    obj.benchmark_sv(1)
+    obj.benchmark_sv(2)
 
     # # VFNS benchmarks with LHA settings
     # programs = ["LHA", "pegasus", "apfel"]

benchmarks/pegasus_bench.py

@@ -83,7 +83,7 @@ def benchmark_sv(self, pto, svmode):
         th.update(
             {
                 "PTO": [pto],
-                "fact_to_ren_scale_ratio": [np.sqrt(0.5), np.sqrt(2.0)],
+                "XIF": [np.sqrt(0.5), np.sqrt(2.0)],
                 "ModSV": [svmode],
             }
         )
@@ -141,7 +141,7 @@ def benchmark_sv(self, pto, svmode):
         th.update(
             {
                 "PTO": [pto],
-                "fact_to_ren_scale_ratio": [np.sqrt(0.5), np.sqrt(2.0)],
+                "XIF": [np.sqrt(0.5), np.sqrt(2.0)],
                 "ModSV": [svmode],
             }
         )
@@ -151,7 +151,6 @@ def benchmark_sv(self, pto, svmode):
 
 
 if __name__ == "__main__":
-
     # obj = BenchmarkVFNS()
     obj = BenchmarkFFNS()
     obj.benchmark_plain_pol(1)

benchmarks/sandbox.py

@@ -67,7 +67,7 @@ def doit(self):
             "QED": 0,
             # "ModEv": "EXA",
             # "XIR": 0.5,
-            # "fact_to_ren_scale_ratio": 2.0,
+            # "XIF": 2.0,
             "Q0": 1.65,  # np.sqrt(10),
             # "Qref": 1.5,
             # "alphas": 0.35,
@@ -116,7 +116,7 @@ def lha(self):
             "FNS": "FFNS",
             "NfFF": 4,
             "ModEv": "EXA",
-            "fact_to_ren_scale_ratio": np.sqrt(2),
+            "XIF": np.sqrt(2),
             "Q0": np.sqrt(2),
             "kcThr": 0.0,
             "kbThr": np.inf,

doc/source/code/IO.rst

@@ -32,12 +32,12 @@ and environment. The benchmark settings are available at :mod:`banana.data.theor
       ``iterate-exact`` abbreviated with ``EXA``, ``decompose-exact``, ``perturbative-exact``,
       ``iterate-expanded`` abbreviated with ``EXP``, ``decompose-expanded``, ``perturbative-expanded``,
       ``truncated`` abbreviated with ``TRN``, ``ordered-truncated``.
-  * - ``fact_to_ren_scale_ratio``
+  * - ``XIF``
     - :py:obj:`float`
     - Factorization to renormalization scale ratio. ``1`` means no scale variation.
   * - ``ModSV``
     - :py:obj:`str`
-    - Scale variation method, used only if ``fact_to_ren_scale_ratio!=1``. Possible options are:
+    - Scale variation method, used only if ``XIF!=1``. Possible options are:
       ``expanded`` or ``exponentiated``.
   * - ``Q0``
     - :py:obj:`float`

extras/matching/check-matching.py

@@ -26,7 +26,6 @@
     "IC": 0,
     "NfFF": 3,
     "Qref": 1.4142135623730951,
-    "fact_to_ren_scale_ratio": 1.0,
     "mb": 4.5,
     "CKM": "0.97428 0.22530 0.003470 0.22520 0.97345 0.041000 0.00862 0.04030 0.999152",
     "IB": 0,

src/eko/evolution_operator/__init__.py

@@ -313,35 +313,23 @@ def grid_size(self):
         """Return the grid size."""
         return self.int_disp.xgrid.size
 
-    def sv_exponentiated_shift(self, q2):
-        """Compute shifted renormalization scale.
-
-        Parameters
-        ----------
-        q2 : float
-            factorization scale
-
-        Returns
-        -------
-        float
-            renormalization scale
-        """
-        if self.sv_mode == sv.Modes.exponentiated:
-            return q2 / self.xif2
-        return q2
-
     @property
     def a_s(self):
         """Return the computed values for :math:`a_s`."""
         sc = self.managers["strong_coupling"]
         a0 = sc.a_s(
-            self.sv_exponentiated_shift(self.q2_from),
-            fact_scale=self.q2_from,
+            self.q2_from,
+            fact_scale=self.q2_from * self.xif2,
             nf_to=self.nf,
         )
         a1 = sc.a_s(
-            self.sv_exponentiated_shift(self.q2_to),
-            fact_scale=self.q2_to,
+            self.q2_to
+            * (
+                self.xif2
+                if not self.is_threshold and self.sv_mode == sv.Modes.expanded
+                else 1.0
+            ),
+            fact_scale=self.q2_to * self.xif2,
             nf_to=self.nf,
         )
         return (a0, a1)
@@ -483,7 +471,9 @@ def compute(self):
             # unless we have to do some scale variation
             # TODO remove if K is factored out of here
             if not (
-                self.sv_mode == sv.Modes.expanded and not np.isclose(self.xif2, 1.0)
+                self.sv_mode == sv.Modes.expanded
+                and not np.isclose(self.xif2, 1.0)
+                and not self.is_threshold
             ):
                 logger.info(
                     "%s: skipping unity operator at %e", self.log_label, self.q2_from
@@ -501,8 +491,8 @@ def compute(self):
         logger.info(
             "%s: µ_R^2 distance: %e -> %e",
             self.log_label,
-            self.sv_exponentiated_shift(self.q2_from),
-            self.sv_exponentiated_shift(self.q2_to),
+            self.q2_from,
+            self.q2_to,
         )
         if self.sv_mode != sv.Modes.unvaried:
             logger.info(

src/eko/evolution_operator/grid.py

@@ -201,11 +201,7 @@ def compute(self, q2grid=None):
         # And now return the grid
         grid_return = {}
         for q2 in q2grid:
-            # shift path for expanded scheme
-            q2_gen = (
-                q2 * self.config["xif2"] if self.sv_mode == sv.Modes.expanded else q2
-            )
-            grid_return[q2] = self.generate(q2_gen)
+            grid_return[q2] = self.generate(q2)
         return grid_return
 
     def generate(self, q2):

src/eko/evolution_operator/operator_matrix_element.py

@@ -322,9 +322,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.sv_exponentiated_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1
-        )
+        return sc.a_s(self.q2_from, self.q2_from * self.xif2, nf_to=self.nf + 1)
 
     def compute(self):
         """Compute the actual operators (i.e. run the integrations)."""

src/eko/runner/legacy.py

@@ -88,6 +88,7 @@ def __init__(
         )
 
         # strong coupling
+        new_theory.couplings.alphas.scale *= new_theory.xif
         sc = Couplings(
             couplings=new_theory.couplings,
             order=new_theory.order,

src/ekomark/benchmark/external/LHA_utils.py

@@ -50,10 +50,11 @@
     ]
 )
 
+
 # rotate basis
 def rotate_data(raw, is_ffns_nnlo=False, rotate_to_evolution_basis=False):
     """Rotate data in flavor space.
-    
+
     Rotate either to flavor basis  or evolution basis from the LHA basis, which is yet an other basis.
 
     Parameters
@@ -140,7 +141,7 @@ def compute_LHA_data(theory, operators, rotate_to_evolution_basis=False):
     with open(here / yaml_file, encoding="utf-8") as o:
         data = yaml.safe_load(o)
     fns = theory["FNS"]
-    xif = (theory["fact_to_ren_scale_ratio"]) ** 2
+    xif = (theory["XIF"]) ** 2
     if order == 0 and xif != 1.0:
         raise ValueError("LO LHA tables with scale variations are not available")
     table = None

src/ekomark/benchmark/external/pegasus_utils.py

@@ -48,9 +48,7 @@ def compute_pegasus_data(theory, operators, skip_pdfs, rotate_to_evolution_basis
         raise ValueError("Initial scale Q0 must be equal to Qref in Pegasus.")
 
     if operators["polarized"]:
-        pegasus.initpol(
-            imodev, theory["PTO"], ivfns, nf, theory["fact_to_ren_scale_ratio"] ** 2
-        )
+        pegasus.initpol(imodev, theory["PTO"], ivfns, nf, theory["XIF"] ** 2)
         pegasus.initpinp(
             theory["alphas"],
             theory["Qref"] ** 2,
@@ -59,9 +57,7 @@ def compute_pegasus_data(theory, operators, skip_pdfs, rotate_to_evolution_basis
             (theory["ktThr"] * theory["mt"]) ** 2,
         )
     else:
-        pegasus.initevol(
-            imodev, theory["PTO"], ivfns, nf, theory["fact_to_ren_scale_ratio"] ** 2
-        )
+        pegasus.initevol(imodev, theory["PTO"], ivfns, nf, theory["XIF"] ** 2)
         pegasus.initinp(
             theory["alphas"],
             theory["Qref"] ** 2,
@@ -87,7 +83,6 @@ def compute_pegasus_data(theory, operators, skip_pdfs, rotate_to_evolution_basis
     # run pegaus
     out_tabs = {}
     for q2 in operators["Q2grid"]:
-
         tab = {}
         for x in target_xgrid:
             # last two numbers are the min and max pid to calculate,

src/ekomark/navigator/navigator.py

@@ -1,3 +1,4 @@
+"""EKO Navigator."""
 import os
 import webbrowser
 
@@ -51,10 +52,9 @@ def fill_theories(self, theo, obj):
             "Qref",
             "alphas",
             "XIF",
-            "XIR",
         ]:
             obj[f] = theo[f]
-        obj["F/R"] = theo["fact_to_ren_scale_ratio"]
+        obj["XIF"] = theo["XIF"]
         obj["mcThr"] = theo["mc"] * theo["kcThr"]
         obj["mbThr"] = theo["mb"] * theo["kbThr"]
         obj["mtThr"] = theo["mt"] * theo["ktThr"]
@@ -133,7 +133,7 @@ def fill_logs(self, lg, obj):
 
     def check_log(self, doc_hash, perc_thr=1, abs_thr=1e-6):
         """
-        Check if the log passed the default assertions
+        Check if the log passed the default assertions.
 
         Parameters
         ----------
@@ -154,12 +154,13 @@ def check_log(self, doc_hash, perc_thr=1, abs_thr=1e-6):
 
     def plot_pdfs(self, doc_hash):
         """
-        Plots all PDFs at the final scale.
+        Plot all PDFs at the final scale.
 
         Parameters
         ----------
-            doc_hash : hash
-                log hash
+        doc_hash : hash
+            log hash
+
         """
         log = self.get(bnav.l, doc_hash)
         dfd = log["log"]
@@ -175,7 +176,6 @@ def plot_pdfs(self, doc_hash):
         print(f"Writing pdf plots to {path}")
 
         with PdfPages(path) as pp:
-
             # print setup
             theory = self.get(bnav.t, log["t_hash"][: self.hash_len])
             ops = self.get(bnav.o, log["o_hash"][: self.hash_len])
@@ -201,12 +201,13 @@ def plot_pdfs(self, doc_hash):
 
     def display_pdfs(self, doc_hash):
         """
-        Display PDF generated by ekomark.navigator.navigator.plot_pdfs
+        Display PDF generated by ekomark.navigator.navigator.plot_pdfs.
 
         Parameters
         ----------
-            doc_hash : hash
-                log hash
+        doc_hash : hash
+            log hash
+
         """
         log = self.get(bnav.l, doc_hash)
         directory = (
@@ -230,6 +231,21 @@ def display_pdfs(self, doc_hash):
 
     @staticmethod
     def is_valid_physical_object(name):
+        """Identify physical objects.
+
+        Used to test names, in order to distinguish physical quantities from metadata.
+
+        Parameters
+        ----------
+        name: str
+            name to test
+
+        Returns
+        -------
+        bool
+            test response
+
+        """
         return name in br.evol_basis or name in br.flavor_basis_names
 
     def compare_external(self, dfd1, dfd2):

tests/eko/evolution_operator/test_grid.py

@@ -81,12 +81,10 @@ def test_mod_expanded(theory_card, theory_ffns, operator_card, tmp_path: pathlib
         theory.order = (1, 0)
         theory.num_flavs_init = nf0
         theory.matching
-        theory.fact_to_ren = 1.0
         path.unlink(missing_ok=True)
         opgrid = legacy.Runner(theory, operator_card, path=path).op_grid
         opg = opgrid.compute(3)
-        theory.fact_to_ren = 1.0 + epsilon
-        theory_update["fact_to_ren_scale_ratio"] = 1.0 + epsilon
+        theory_update["XIF"] = 1.0 + epsilon
         path.unlink(missing_ok=True)
         sv_opgrid = legacy.Runner(theory, operator_card, path=path).op_grid
         sv_opg = sv_opgrid.compute(3)

tests/eko/evolution_operator/test_init.py

@@ -222,7 +222,6 @@ 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.sv_exponentiated_shift(40.0), 10.0)
         o.compute()
         self.check_lo(o)