Thresholds rewrite

benchmarks/ekomark/LHA.py

@@ -93,7 +93,7 @@ class LHABenchmarkPaper(Runner):
     def __init__(self, theory_path, operators_path, assets_dir, data_dir):
         super().__init__(theory_path, operators_path, assets_dir)
 
-        self.rotate_to_evolution_basis = False  # True
+        self.rotate_to_evolution_basis = not True
 
         if not np.isclose(self.theory["XIF"], 1.0):
             raise ValueError("XIF has to be 1")

benchmarks/input/LHA/theories/LO-EXA-FFNS.yaml

@@ -21,13 +21,13 @@ SxOrd: "LL"
 HQ: "POLE"
 mc: 1.4142135623730951
 Qmc: 1.4142135623730951
-kcThr: 1.0
+kcThr: 0.0
 mb: 4.5
 Qmb: 4.5
-kbThr: 1.0
+kbThr: 1.e+100
 mt: 175.0
 Qmt: 175.0
-ktThr: 1.0
+ktThr: 1.e+100
 CKM: "0.97428 0.22530 0.003470 0.22520 0.97345 0.041000 0.00862 0.04030 0.999152"
 MZ: 91.1876
 MW: 80.398

doc/source/code/Utilities.rst

@@ -3,7 +3,7 @@ Utility Classes
 
 Apart from the :doc:`operator <Operators>` classes, `eko` also provides some utility classes which are e.g. used in |yadism|
 
-- :class:`eko.thresholds.ThresholdsConfig`
+- :class:`eko.thresholds.ThresholdsAtlas`
 
   -  Implementation of the flavor number scheme and the quark thresholds both for
      the :class:`eko.strong_coupling.StrongCoupling` and the :doc:`operators <../theory/Matching>`

src/eko/operator/grid.py

@@ -1,7 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-This module contains the :class:`OperatorGrid`  and the
-:class:`OperatorMaster` class.
+This module contains the :class:`OperatorGrid` class.
 
 The first is the driver class of eko as it is the one that collects all the
 previously instantiated information and does the actual computation of the Q2s.
@@ -31,7 +30,7 @@ class OperatorGrid:
             Grid in Q2 on where to to compute the operators
         order: int
             order in perturbation theory
-        thresholds_config: eko.thresholds.ThresholdsConfig
+        thresholds_config: eko.thresholds.ThresholdsAtlas
             Instance of :class:`~eko.thresholds.Threshold` containing information about the
             thresholds
         strong_coupling: eko.strong_coupling.StrongCoupling
@@ -73,7 +72,6 @@ def __init__(
             strong_coupling=strong_coupling,
             interpol_dispatcher=interpol_dispatcher,
         )
-        self._op_grid = {}
         self._threshold_operators = {}
 
     @classmethod
@@ -88,18 +86,12 @@ def from_dict(
         """
         Create the object from the theory dictionary.
 
-        Read keys:
-
-            - Q2grid : required, target grid
-            - debug_skip_singlet : optional, avoid singlet integrations?
-            - debug_skip_non_singlet : optional, avoid non-singlet integrations?
-
         Parameters
         ----------
-            setup : dict
+            theory_card : dict
                 theory dictionary
-            thresholds_config : eko.thresholds.ThresholdsConfig
-                An instance of the ThresholdsConfig class
+            thresholds_config : eko.thresholds.ThresholdsAtlas
+                An instance of the ThresholdsAtlas class
             strong_coupling : eko.strong_coupling.StrongCoupling
                 An instance of the StrongCoupling class
             interpol_dispatcher : eko.interpolation.InterpolatorDispatcher
@@ -134,9 +126,9 @@ def from_dict(
             config, q2_grid, thresholds_config, strong_coupling, interpol_dispatcher
         )
 
-    def _generate_thresholds_op(self, to_q2):
+    def get_threshold_operators(self, path):
         """
-        Generate the threshold operators
+        Generate the threshold operators.
 
         This method is called everytime the OperatorGrid is asked for a grid on Q^2
         with a list of the relevant areas.
@@ -148,171 +140,86 @@ def _generate_thresholds_op(self, to_q2):
 
         Parameters
         ----------
-            to_q2: float
-                value of q2 for which the OperatorGrid will need to pass thresholds
+            path: list(PathSegment)
+                thresholds path
         """
-        # The lists of areas as produced by the thresholds
-        area_list = self.managers["thresholds_config"].get_path_from_q2_ref(to_q2)
         # The base area is always that of the reference q
-        q2_from = self.managers["thresholds_config"].q2_ref
-        nf = self.managers["thresholds_config"].nf_ref
-        for area in area_list:
-            q2_to = area.q2_ref
-            # TODO due to this continue we're actually seeing the area *one below*
-            if q2_to == q2_from:
-                continue
-            new_op_key = (q2_from, q2_to)
-            logger.info(
-                "Evolution: Compute threshold operator from %e to %e", q2_from, q2_to
-            )
+        thr_ops = []
+        for seg in path[:-1]:
+            new_op_key = seg.tuple
             if new_op_key not in self._threshold_operators:
-                # Compute the operator in place and store it
-                op_th = Operator(self.config, self.managers, nf, q2_from, q2_to)
+                # Compute the operator and store it
+                logger.info(
+                    "Threshold operator: %e -> %e, nf=%d",
+                    seg.q2_from,
+                    seg.q2_to,
+                    seg.nf,
+                )
+                op_th = Operator(
+                    self.config, self.managers, seg.nf, seg.q2_from, seg.q2_to
+                )
                 op_th.compute()
                 self._threshold_operators[new_op_key] = op_th
-            nf = area.nf
-            q2_from = q2_to
-
-    def _get_jumps(self, qsq):
-        """
-        Given a value of q^2, generates the list of operators that need to be
-        composed in order to get there from q0^2
-
-        Parameters
-        ----------
-            qsq: float
-                Target value of q^2
-
-        Returns
-        -------
-            op_list: list
-                List of threshold operators
-        """
-        # Get the list of areas to be crossed
-        full_area_path = self.managers["thresholds_config"].get_path_from_q2_ref(qsq)
-        # The last one is where q resides so it is not needed
-        area_path = full_area_path[:-1]
-        op_list = []
-        # Now loop over the areas to collect the necessary threshold operators
-        for area in area_path:
-            q2_from = area.q2_ref
-            q2_to = area.q2_towards(qsq)
-            op_list.append(self._threshold_operators[(q2_from, q2_to)])
-        return op_list
-
-    def set_q2_limits(self, q2min, q2max):
-        """
-        Sets up the limits of the grid in q^2 to be computed by the OperatorGrid
+            thr_ops.append(self._threshold_operators[new_op_key])
+        return thr_ops
 
-        This function is a wrapper to compute the necessary operators to go between areas
-
-        Parameters
-        ----------
-            q2min: float
-                Minimum value of q^2 that will be computed
-            q2max: float
-                Maximum value of q^2 that will be computed
-        """
-        # Sanity checks
-        if q2min <= 0.0 or q2max <= 0.0:
-            raise ValueError(
-                f"Values of q^2 below 0.0 are not accepted, received [{q2min},{q2max}]"
-            )
-        if q2min > q2max:
-            raise ValueError(
-                f"Minimum q^2 is above maximum q^2 (error: {q2max} < {q2min})"
-            )
-        # Ensure we have all the necessary operators to go from q2ref to q2min and qmax
-        self._generate_thresholds_op(q2min)
-        self._generate_thresholds_op(q2max)
-
-    def _compute_raw_grid(self, q2grid):
-        """
-        Receives a grid in q^2 and computes each opeator inside its
-        area with reference value the q_ref of its area
-
-        Parameters
-        ----------
-            q2grid: list
-                List of q^2
+    def compute(self, q2grid=None):
         """
-        area_list = self.managers["thresholds_config"].get_areas(q2grid)
-        for area, q2 in zip(area_list, q2grid):
-            q2_from = area.q2_ref
-            nf = area.nf
-            logger.info("Evolution: compute operators %e -> %e", q2_from, q2)
-            op = Operator(self.config, self.managers, nf, q2_from, q2)
-            op.compute()
-            self._op_grid[q2] = op
-
-    def compute_q2grid(self, q2grid=None):
-        """
-        Receives a grid in q^2 and computes all operations necessary
-        to return any operator at any given q for the evolution between q2ref and q2grid
+        Computes all ekos for the q2grid.
 
         Parameters
         ----------
-            q2grid: list
+            q2grid: list(float)
                 List of q^2
 
         Returns
         -------
-            grid_return: list
-                List of PhysicalOperator for each value of q^2
+            grid_return: list(dict)
+                List of ekos for each value of q^2
         """
         # use input?
         if q2grid is None:
             q2grid = self.q2_grid
         # normalize input
         if isinstance(q2grid, numbers.Number):
             q2grid = [q2grid]
-        # Check max and min of the grid and reset the limits if necessary
-        q2max = np.max(q2grid)
-        q2min = np.min(q2grid)
-        self.set_q2_limits(q2min, q2max)
-        # Now compute all raw operators
-        self._compute_raw_grid(q2grid)
         # And now return the grid
         grid_return = {}
         for q2 in q2grid:
-            grid_return[q2] = self.get_op_at_q2(q2)
+            grid_return[q2] = self.generate(q2)
         return grid_return
 
-    def get_op_at_q2(self, qsq):
+    def generate(self, q2):
         """
-        Given a value of q^2, returns the PhysicalOperator to get
-        to q^2 from q0^2
+        Computes an single EKO.
 
         Parameters
         ----------
-            qsq: float
+            q2: float
                 Target value of q^2
 
         Returns
         -------
-            final_op: eko.evolution_operator.PhysicalOperator
-                Op(q^2 <- q_0^2)
+            final_op: dict
+                eko E(q^2 <- q_0^2) in flavor basis as numpy array
         """
-        # Check the path to q0 for this operator
-        if qsq in self._op_grid:
-            operator = self._op_grid[qsq]
-        else:
-            logger.warning("Evolution: Q2=%e not found in the grid, computing...", qsq)
-            self.compute_q2grid(qsq)
-            operator = self._op_grid[qsq]
+        # The lists of areas as produced by the thresholds
+        path = self.managers["thresholds_config"].path(q2)
         # Prepare the path for the composition of the operator
-        operators_to_q2 = self._get_jumps(operator.q2_from)
-        is_vfns = self.managers["thresholds_config"].scheme != "FFNS"
+        thr_ops = self.get_threshold_operators(path)
+        operator = Operator(
+            self.config, self.managers, path[-1].nf, path[-1].q2_from, path[-1].q2_to
+        )
+        operator.compute()
         final_op = physical.PhysicalOperator.ad_to_evol_map(
             operator.op_members,
             operator.nf,
             operator.q2_to,
-            is_vfns,
             self.config["intrinsic_range"],
         )
-        for op in reversed(operators_to_q2):
+        for op in reversed(thr_ops):
             phys_op = physical.PhysicalOperator.ad_to_evol_map(
-                op.op_members, op.nf, op.q2_to, is_vfns, self.config["intrinsic_range"]
+                op.op_members, op.nf, op.q2_to, self.config["intrinsic_range"]
             )
             final_op = final_op @ phys_op
         values, errors = final_op.to_flavor_basis_tensor()

src/eko/operator/physical.py

@@ -27,7 +27,7 @@ def __init__(self, op_members, q2_final):
         self.q2_final = q2_final
 
     @classmethod
-    def ad_to_evol_map(cls, op_members, nf, q2_final, is_vfns, intrinsic_range=None):
+    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range=None):
         """
         Obtain map between the 3-dimensional anomalous dimension basis and the
         4-dimensional evolution basis.
@@ -41,8 +41,6 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, is_vfns, intrinsic_range=None)
                 operator members in anomalous dimension basis
             nf : int
                 number of active light flavors
-            is_vfns : bool
-                is |VFNS|?
             intrinsic_range : sequence
                 intrinsic heavy flavours
 
@@ -66,39 +64,29 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, is_vfns, intrinsic_range=None)
             m[f"T{n}.T{n}"] = op_members["NS_p"]
         # activate one higher element, i.e. where the next heavy quark could participate,
         # but actually it is trivial
-        if is_vfns:
-            n = (nf + 1) ** 2 - 1
-            # without this new heavy quark Vn = V and Tn = S
-            m[f"V{n}.V"] = op_members["NS_v"]
-            m[f"T{n}.S"] = op_members["S_qq"]
-            m[f"T{n}.g"] = op_members["S_qg"]
+        n = (nf + 1) ** 2 - 1
+        # without this new heavy quark Vn = V and Tn = S
+        m[f"V{n}.V"] = op_members["NS_v"]
+        m[f"T{n}.S"] = op_members["S_qq"]
+        m[f"T{n}.g"] = op_members["S_qg"]
         # deal with intrinsic heavy quark pdfs
         if intrinsic_range is not None:
             hqfl = "cbt"
-            for f in intrinsic_range:
-                hq = hqfl[f - 4]  # find name
+            for intr_fl in intrinsic_range:
+                hq = hqfl[intr_fl - 4]  # find name
                 # intrinsic means no evolution, i.e. they are evolving with the identity
                 len_xgrid = op_members["NS_v"].value.shape[0]
                 op_id = member.OpMember(
                     np.eye(len_xgrid), np.zeros((len_xgrid, len_xgrid))
                 )
-                if f > nf + 1:  # keep the higher quarks as they are
+                if intr_fl > nf + 1:  # keep the higher quarks as they are
                     m[f"{hq}+.{hq}+"] = op_id
                     m[f"{hq}-.{hq}-"] = op_id
-                elif f == nf + 1:  # next is comming hq?
-                    n = f ** 2 - 1
-                    if is_vfns:
-                        # e.g. T15 = (u+ + d+ + s+) - 3c+
-                        m[f"V{n}.{hq}-"] = -(f - 1) * op_id
-                        m[f"T{n}.{hq}+"] = -(f - 1) * op_id
-                    else:
-                        m[f"{hq}+.{hq}+"] = op_id
-                        m[f"{hq}-.{hq}-"] = op_id
-                else:  # f <= nf
-                    if not is_vfns:
-                        raise ValueError(
-                            f"{hq} is perturbative inside FFNS{nf} so can NOT be intrinsic"
-                        )
+                elif intr_fl == nf + 1:  # next is comming hq?
+                    n = intr_fl ** 2 - 1
+                    # e.g. T15 = (u+ + d+ + s+) - 3c+
+                    m[f"V{n}.{hq}-"] = -(intr_fl - 1) * op_id
+                    m[f"T{n}.{hq}+"] = -(intr_fl - 1) * op_id
         opms = {}
         for k, v in m.items():
             opms[flavors.MemberName(k)] = v.copy()