From 89d2de747beb214eedcb946d3940b9813cbfc6b1 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Thu, 26 Nov 2020 18:08:53 +0100 Subject: [PATCH 01/12] Init thresholds rewrite --- src/eko/thresholds.py | 156 ++++++++---------------------------------- 1 file changed, 28 insertions(+), 128 deletions(-) diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index d999927ae..428168a79 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -16,7 +16,7 @@ class Area: """ - Sets up a single threhold area with a fixed configuration. + Sets up a single threshold area with a fixed configuration. Parameters ---------- @@ -24,73 +24,32 @@ class Area: lower bound of the area q2_max : float upper bound of the area - q2_0 : float - reference point of the area (can be anywhere in the area) nf : float number of flavours in the area """ - def __init__(self, q2_min, q2_max, q2_0, nf): + def __init__(self, q2_min, q2_max, nf): self.q2_min = q2_min self.q2_max = q2_max self.nf = nf - self.has_q2_0 = False - # Now check which is the q2_ref for this area - if q2_0 > q2_max: - self.q2_ref = q2_max - elif q2_0 < q2_min: - self.q2_ref = q2_min - else: - self.has_q2_0 = True - self.q2_ref = q2_0 - - def q2_towards(self, q2): - """ - Return q2_min or q2_max depending on whether - we are going towards the max or the min or q2 - if we are alreay in the correct area - - Parameters - ---------- - q2 : float - reference point - - Returns - ------- - q2_next : float - the closest point to q2 that is within the area - """ - if q2 > self.q2_max: - return self.q2_max - elif q2 < self.q2_min: - return self.q2_min - else: - return q2 - def __gt__(self, target_area): - """Compares q2 of areas""" - return self.q2_min >= target_area.q2_max - - def __lt__(self, target_area): - """Compares q2 of areas""" - return self.q2_max <= target_area.q2_min - - def __call__(self, q2): - """ - Checks whether q2 is contained in the area - - Parameters - ---------- - q2 : float - testing point - - Returns - ------- - contained : bool - is point contained? - """ - return self.q2_min <= q2 <= self.q2_max +class PathSegment: + """ + Oriented path in the threshold area landscape. + Parameters + ---------- + q2_from : float + starting point + q2_to : float + final point + area : Area + containing area + """ + def __init__(self, q2_from, q2_to, area): + self.q2_from = q2_from + self.q2_to = q2_to + self.area = area class ThresholdsConfig: """ @@ -109,39 +68,18 @@ class ThresholdsConfig: Number of flavors for the FFNS """ - def __init__(self, q2_ref, scheme, *, threshold_list=None, nf=None): + def __init__(self, area_walls, q2_ref = None): # Initial values self.q2_ref = q2_ref - self.scheme = scheme - self._areas = [] - self._area_walls = [] - self._area_ref = 0 - self.max_nf = None - self.min_nf = None - - if scheme == "FFNS": - if nf is None: - raise ValueError("No value for nf in the FFNS was received") - if threshold_list is not None: - raise ValueError("The FFNS does not accept any thresholds") - self._areas = [Area(0, np.inf, self.q2_ref, nf)] - self.max_nf = nf - self.min_nf = nf - logger.info("Fixed flavor number scheme with %d flavors", nf) - elif scheme in ["ZM-VFNS", "FONLL-A", "FONLL-A'"]: - if scheme not in ["FONLL-A"] and nf is not None: - logger.warning( - "The VFNS configures its own value for nf, ignoring input nf=%d", - nf, - ) - if threshold_list is None: - raise ValueError( - "The VFN scheme was selected but no thresholds were given" - ) - self._setup_vfns(threshold_list, nf) - logger.info("Variable flavor number scheme (%s)", scheme) - else: - raise NotImplementedError(f"The scheme {scheme} is not implemented") + if area_walls != sorted(area_walls): + raise ValueError("area walls need to be sorted") + self.areas = [] + q2_min = 0 + for q2_max in area_walls + [np.inf]: + nf = nf + 1 + new_area = Area(q2_min, q2_max, nf) + self.areas.append(new_area) + q2_min = q2_max @classmethod def from_dict(cls, theory_card): @@ -182,44 +120,6 @@ def nf_ref(self): """ Number of flavours in the reference area """ return self._areas[self._area_ref].nf - def nf_range(self): - """ - Iterate number of flavours, including min_nf *and* max_nf - - Yields - ------ - nf : int - number of flavour - """ - return range(self.min_nf, self.max_nf + 1) - - def _setup_vfns(self, threshold_list, nf=None): - """ - Receives a list of thresholds and sets up the vfns scheme - - Parameters - ---------- - threshold_list: list - List of q^2 thresholds - """ - if nf is None: - nf = 3 - # Force sorting - self._area_walls = sorted(threshold_list) - # Generate areas - self._areas = [] - q2_min = 0 - q2_ref = self.q2_ref - self.min_nf = nf - for i, q2_max in enumerate(self._area_walls + [np.inf]): - nf = self.min_nf + i - new_area = Area(q2_min, q2_max, q2_ref, nf) - if new_area.has_q2_0: - self._area_ref = i - self._areas.append(new_area) - q2_min = q2_max - self.max_nf = nf - def get_path_from_q2_ref(self, q2): """ Get the Area path from q2_ref to q2. From 3d337994888fec795a0b06c53517dc046aa82370 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 27 Nov 2020 18:20:07 +0100 Subject: [PATCH 02/12] Rework threshold and strong coupling --- benchmarks/ekomark/LHA.py | 2 +- .../input/LHA/theories/LO-EXA-FFNS.yaml | 6 +- doc/source/code/Utilities.rst | 2 +- src/eko/operator/grid.py | 49 ++-- src/eko/runner.py | 4 +- src/eko/strong_coupling.py | 37 ++- src/eko/thresholds.py | 163 ++++++------ tests/benchmark_strong_coupling.py | 24 +- tests/test_op_grid.py | 6 +- tests/test_operator.py | 4 +- tests/test_strong_coupling.py | 38 +-- tests/test_thresholds.py | 251 +++++++++++------- 12 files changed, 300 insertions(+), 286 deletions(-) diff --git a/benchmarks/ekomark/LHA.py b/benchmarks/ekomark/LHA.py index 84c5f3a16..cd1e1a979 100644 --- a/benchmarks/ekomark/LHA.py +++ b/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") diff --git a/benchmarks/input/LHA/theories/LO-EXA-FFNS.yaml b/benchmarks/input/LHA/theories/LO-EXA-FFNS.yaml index c142f7964..09a9957e9 100644 --- a/benchmarks/input/LHA/theories/LO-EXA-FFNS.yaml +++ b/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 diff --git a/doc/source/code/Utilities.rst b/doc/source/code/Utilities.rst index 3d8cbc49e..083ace8c3 100644 --- a/doc/source/code/Utilities.rst +++ b/doc/source/code/Utilities.rst @@ -3,7 +3,7 @@ Utility Classes Apart from the :doc:`operator ` 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>` diff --git a/src/eko/operator/grid.py b/src/eko/operator/grid.py index 6d3cea05f..5fef1d58b 100644 --- a/src/eko/operator/grid.py +++ b/src/eko/operator/grid.py @@ -31,7 +31,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 @@ -98,8 +98,8 @@ def from_dict( ---------- setup : 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 @@ -152,26 +152,28 @@ def _generate_thresholds_op(self, to_q2): value of q2 for which the OperatorGrid will need to pass thresholds """ # The lists of areas as produced by the thresholds - area_list = self.managers["thresholds_config"].get_path_from_q2_ref(to_q2) + path = self.managers["thresholds_config"].path(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 + # q2_from = self.managers["thresholds_config"].q2_ref + # nf = self.managers["thresholds_config"].nf_ref + for seg in path[:-1]: + # q2_to = area.q2_ref # TODO due to this continue we're actually seeing the area *one below* - if q2_to == q2_from: + if seg.q2_to == seg.q2_from: continue - new_op_key = (q2_from, q2_to) + new_op_key = (seg.q2_from, seg.q2_to) logger.info( - "Evolution: Compute threshold operator from %e to %e", q2_from, q2_to + "Evolution: Compute threshold operator from %e to %e", + seg.q2_from, + seg.q2_to, ) 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) + 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): """ @@ -189,15 +191,12 @@ def _get_jumps(self, qsq): 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) + path = self.managers["thresholds_config"].path(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)]) + for seg in path[:-1]: + op_list.append(self._threshold_operators[(seg.q2_from, seg.q2_to)]) return op_list def set_q2_limits(self, q2min, q2max): @@ -236,10 +235,10 @@ def _compute_raw_grid(self, q2grid): q2grid: list List of q^2 """ - 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 + for q2 in q2grid: + path = self.managers["thresholds_config"].path(q2) + q2_from = path[-1].q2_from + nf = path[-1].nf logger.info("Evolution: compute operators %e -> %e", q2_from, q2) op = Operator(self.config, self.managers, nf, q2_from, q2) op.compute() @@ -302,7 +301,7 @@ def get_op_at_q2(self, qsq): operator = self._op_grid[qsq] # 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" + is_vfns = True # TODO #self.managers["thresholds_config"].scheme != "FFNS" final_op = physical.PhysicalOperator.ad_to_evol_map( operator.op_members, operator.nf, diff --git a/src/eko/runner.py b/src/eko/runner.py index 3e3303b4a..39a540678 100644 --- a/src/eko/runner.py +++ b/src/eko/runner.py @@ -8,7 +8,7 @@ from . import interpolation from .output import Output from .strong_coupling import StrongCoupling -from .thresholds import ThresholdsConfig +from .thresholds import ThresholdsAtlas from .operator.grid import OperatorGrid from . import basis_rotation as br @@ -79,7 +79,7 @@ def __init__(self, theory_card, operators_card): bfd = interpolation.InterpolatorDispatcher.from_dict(operators_card) self.out.update(bfd.to_dict()) # FNS - tc = ThresholdsConfig.from_dict(theory_card) + tc = ThresholdsAtlas.from_dict(theory_card) self.out["q2_ref"] = float(tc.q2_ref) # strong coupling sc = StrongCoupling.from_dict(theory_card, tc) diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index eb41b8e05..1eb87bff9 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -99,8 +99,8 @@ class StrongCoupling: alpha_s(!) at the reference scale :math:`\alpha_s(\mu_0^2)` scale_ref : float reference scale :math:`\mu_0^2` - threshold_holder : eko.thresholds.ThresholdsConfig - An instance of :class:`~eko.thresholds.ThresholdsConfig` + threshold_holder : eko.thresholds.ThresholdsAtlas + An instance of :class:`~eko.thresholds.ThresholdsAtlas` order: int Evaluated order of the beta function: ``0`` = LO, ... method : ["expanded", "exact"] @@ -110,7 +110,7 @@ class StrongCoupling: -------- >>> alpha_ref = 0.35 >>> scale_ref = 2 - >>> threshold_holder = ThresholdsConfig( ... ) + >>> threshold_holder = ThresholdsAtlas( ... ) >>> sc = StrongCoupling(alpha_ref, scale_ref, threshold_holder) >>> q2 = 91.1**2 >>> sc.a_s(q2) @@ -130,7 +130,7 @@ def __init__( raise ValueError(f"alpha_s_ref has to be positive - got {alpha_s_ref}") if scale_ref <= 0: raise ValueError(f"scale_ref has to be positive - got {scale_ref}") - if not isinstance(thresh, thresholds.ThresholdsConfig): + if not isinstance(thresh, thresholds.ThresholdsAtlas): raise ValueError("Needs a Threshold instance") if order not in [0, 1, 2]: raise NotImplementedError("a_s beyond NNLO is not implemented") @@ -141,14 +141,9 @@ def __init__( # create new threshold object self.as_ref = alpha_s_ref / 4.0 / np.pi # convert to a_s - if thresh.scheme == "FFNS": - self._threshold_holder = thresholds.ThresholdsConfig( - scale_ref, thresh.scheme, nf=thresh.nf_ref - ) - else: - self._threshold_holder = thresholds.ThresholdsConfig( - scale_ref, thresh.scheme, threshold_list=thresh._area_walls - ) + self._threshold_holder = thresholds.ThresholdsAtlas( + thresh.area_walls[1:-1], scale_ref + ) logger.info( "Strong Coupling: Reference a_s(Q^2=%f)=%f", self.q2_ref, self.as_ref ) @@ -174,7 +169,7 @@ def from_dict(cls, theory_card, thresholds_config=None): ---------- theory_card : dict theory dictionary - thresholds_config : eko.thresholds.ThresholdsConfig + thresholds_config : eko.thresholds.ThresholdsAtlas threshold configuration Returns @@ -204,7 +199,7 @@ def from_dict(cls, theory_card, thresholds_config=None): raise ValueError(f"Unknown evolution mode {mod_ev}") # eventually read my dependents if thresholds_config is None: - thresholds_config = thresholds.ThresholdsConfig.from_dict(theory_card) + thresholds_config = thresholds.ThresholdsAtlas.from_dict(theory_card) return cls(alpha_ref, q2_alpha, thresholds_config, order, method) def _compute_exact(self, as_ref, nf, scale_from, scale_to): @@ -298,20 +293,20 @@ def a_s(self, scale_to, fact_scale=None): """ # Set up the path to follow in order to go from q2_0 to q2_ref final_as = self.as_ref - area_path = self._threshold_holder.get_path_from_q2_ref(scale_to) + path = self._threshold_holder.path(scale_to) # as a default assume mu_F^2 = mu_R^2 if fact_scale is None: fact_scale = scale_to - for k, area in enumerate(area_path): - q2_from = area.q2_ref - q2_to = area.q2_towards(scale_to) + for k, seg in enumerate(path): + q2_from = seg.q2_from + q2_to = seg.q2_to if np.isclose(q2_from, q2_to): continue - new_as = self._compute(final_as, area.nf, q2_from, q2_to) + new_as = self._compute(final_as, seg.nf, q2_from, q2_to) # apply matching conditions: see hep-ph/9706430 # - if there is yet a step to go - if k < len(area_path) - 1: - next_nf_is_down = area_path[k + 1].nf < area.nf + if k < len(path) - 1: + next_nf_is_down = path[k + 1].nf < seg.nf # q2_to is the threshold value L = np.log(scale_to / fact_scale) if next_nf_is_down: diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index 428168a79..844dca8dd 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -7,7 +7,6 @@ .. include:: /code/IO-tabs/ThresholdConfig.rst """ import logging -import numbers import numpy as np @@ -33,6 +32,10 @@ def __init__(self, q2_min, q2_max, nf): self.q2_max = q2_max self.nf = nf + def __repr__(self): + return f"Area([{self.q2_min},{self.q2_max}], nf={self.nf})" + + class PathSegment: """ Oriented path in the threshold area landscape. @@ -46,43 +49,72 @@ class PathSegment: area : Area containing area """ + def __init__(self, q2_from, q2_to, area): self.q2_from = q2_from self.q2_to = q2_to - self.area = area + self._area = area + + @property + def nf(self): + return self._area.nf + + def __repr__(self): + return f"PathSegment({self.q2_from} -> {self.q2_to}, nf={self.nf})" -class ThresholdsConfig: + @classmethod + def intersect(cls, q2_from, q2_to, area): + if q2_from < q2_to: + return cls(max(q2_from, area.q2_min), min(q2_to, area.q2_max), area) + else: + return cls(min(q2_from, area.q2_max), max(q2_to, area.q2_min), area) + + +class ThresholdsAtlas: """ The threshold class holds information about the thresholds any Q2 has to pass in order to get there from a given q2_ref and scheme. Parameters ---------- - q2_ref : float - Reference q^2 - scheme: str - Choice of scheme - threshold_list: list - List of q^2 thresholds should the scheme accept it - nf: int - Number of flavors for the FFNS + thresholds: list(float) + List of q^2 thresholds + q2_ref: float + reference scale """ - def __init__(self, area_walls, q2_ref = None): + def __init__(self, thresholds, q2_ref=None): # Initial values self.q2_ref = q2_ref - if area_walls != sorted(area_walls): - raise ValueError("area walls need to be sorted") + thresholds = list(thresholds) + if thresholds != sorted(thresholds): + raise ValueError("thresholds need to be sorted") self.areas = [] + self.area_walls = [0] + thresholds + [np.inf] q2_min = 0 - for q2_max in area_walls + [np.inf]: - nf = nf + 1 + nf = 3 + for q2_max in self.area_walls[1:]: new_area = Area(q2_min, q2_max, nf) self.areas.append(new_area) + nf = nf + 1 q2_min = q2_max @classmethod - def from_dict(cls, theory_card): + def ffns(cls, nf, q2_ref=None): + """ + Create a |FFNS| setup. + + Parameters + ---------- + nf : int + number of light flavors + q2_ref : float + reference scale + """ + return cls([0] * (nf - 3) + [np.inf] * (6 - nf), q2_ref) + + @classmethod + def from_dict(cls, theory_card, prefix="k"): """ Create the object from the run card. @@ -96,91 +128,46 @@ def from_dict(cls, theory_card): cls : ThresholdConfig created object """ - FNS = theory_card["FNS"] - q2_ref = pow(theory_card["Q0"], 2) - if FNS == "FFNS": - return cls(q2_ref, FNS, threshold_list=None, nf=theory_card["NfFF"]) - # the threshold value does not necessarily coincide with the mass + def thres(pid): heavy_flavors = "cbt" flavor = heavy_flavors[pid - 4] - return pow(theory_card[f"m{flavor}"] * theory_card[f"k{flavor}Thr"], 2) - if FNS in ["FONLL-A"]: - nf = theory_card["NfFF"] - if nf <= 3: - raise ValueError("NfFF should point to the heavy quark! i.e. NfFF>3") - threshold_list = [thres(nf)] - return cls(q2_ref, FNS, threshold_list=threshold_list, nf=nf-1) - # setup VFNS - threshold_list = [thres(q) for q in range(4, 6 + 1)] - return cls(q2_ref, FNS, threshold_list=threshold_list) + return pow( + theory_card[f"m{flavor}"] * theory_card[f"{prefix}{flavor}Thr"], 2 + ) - @property - def nf_ref(self): - """ Number of flavours in the reference area """ - return self._areas[self._area_ref].nf + thresholds = [thres(q) for q in range(4, 6 + 1)] + # preset ref scale + q2_ref = pow(theory_card["Q0"], 2) + return cls(thresholds, q2_ref) - def get_path_from_q2_ref(self, q2): + def path(self, q2_to, q2_from=None): """ - Get the Area path from q2_ref to q2. + Get path from q2_from to q2_to. Parameters ---------- - q2: float - Target value of q2 + q2_to: float + target value of q2 + q2_from: float + starting value of q2 Returns ------- - area_path: list - List of Areas to go through in order to get from q2_ref - to q2. The first one is the one containg q2_ref while the - last one contains q2 + path: list(PathSegment) + List of PathSegment to go through in order to get from q2_from + to q2_to. """ - current_area = self.get_areas_idx(q2)[0] - if current_area < self._area_ref: + if q2_from is None: + q2_from = self.q2_ref + + ref_idx = np.digitize(q2_from, self.area_walls) + target_idx = np.digitize(q2_to, self.area_walls) + if q2_to < q2_from: rc = -1 else: rc = 1 - area_path = [ - self._areas[i] for i in range(self._area_ref, current_area + rc, rc) + return [ + PathSegment.intersect(q2_from, q2_to, self.areas[i - 1]) + for i in range(ref_idx, target_idx + rc, rc) ] - return area_path - - def get_areas_idx(self, q2arr): - """ - Returns the index of the area in which each value of q2arr falls. - - Parameters - ---------- - q2arr: np.array - array of values of q2 - - Returns - ------- - areas_idx: list - list with the indices of the corresponding areas for q2arr - """ - # Ensure q2arr is an array - if isinstance(q2arr, numbers.Number): - q2arr = np.array([q2arr]) - # Check in which area is every q2 - areas_idx = np.digitize(q2arr, self._area_walls) - return areas_idx - - def get_areas(self, q2arr): - """ - Returns the Areas in which each value of q2arr falls - - Parameters - ---------- - q2arr: np.array - array of values of q2 - - Returns - ------- - areas: list - list with the areas for q2arr - """ - idx = self.get_areas_idx(q2arr) - area_list = np.array(self._areas)[idx] - return list(area_list) diff --git a/tests/benchmark_strong_coupling.py b/tests/benchmark_strong_coupling.py index c2ff43e23..0f72f4dd7 100644 --- a/tests/benchmark_strong_coupling.py +++ b/tests/benchmark_strong_coupling.py @@ -33,7 +33,7 @@ def test_a_s(self): ref_alpha_s = 0.1181 ref_mu2 = 90 ask_q2 = 125 - threshold_holder = thresholds.ThresholdsConfig(ref_mu2, "FFNS", nf=5) + threshold_holder = thresholds.ThresholdsAtlas([0, 0, np.inf], ref_mu2) as_FFNS_LO = StrongCoupling(ref_alpha_s, ref_mu2, threshold_holder, order=0) # check local np.testing.assert_approx_equal( @@ -48,16 +48,14 @@ def benchmark_LHA_paper(self): # LO - FFNS # note that the LO-FFNS value reported in :cite:`Giele:2002hx` # was corrected in :cite:`Dittmar:2005ed` - threshold_holder = thresholds.ThresholdsConfig(2, "FFNS", nf=4) + threshold_holder = thresholds.ThresholdsAtlas((0, np.inf, np.inf)) as_FFNS_LO = StrongCoupling(0.35, 2, threshold_holder, order=0) me = as_FFNS_LO.a_s(1e4) * 4 * np.pi ref = 0.117574 np.testing.assert_approx_equal(me, ref, significant=6) # LO - VFNS threshold_list = [2, pow(4.5, 2), pow(175, 2)] - threshold_holder = thresholds.ThresholdsConfig( - 2, "ZM-VFNS", threshold_list=threshold_list - ) + threshold_holder = thresholds.ThresholdsAtlas(threshold_list) as_VFNS_LO = StrongCoupling(0.35, 2, threshold_holder, order=0) me = as_VFNS_LO.a_s(1e4) * 4 * np.pi ref = 0.122306 @@ -95,7 +93,7 @@ def benchmark_APFEL_ffns(self): ), } # collect my values - threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=nf) + threshold_holder = thresholds.ThresholdsAtlas.ffns(nf) for order in [0, 1, 2]: as_FFNS = StrongCoupling( alphas_ref, @@ -162,9 +160,7 @@ def benchmark_APFEL_vfns(self): ), } # collect my values - threshold_holder = thresholds.ThresholdsConfig( - scale_ref, "ZM-VFNS", threshold_list=threshold_list - ) + threshold_holder = thresholds.ThresholdsAtlas(threshold_list) for order in [0, 1, 2]: as_VFNS = StrongCoupling( alphas_ref, @@ -210,7 +206,7 @@ def benchmark_lhapdf_ffns_lo(self): scale_ref = 91.0 ** 2 nf = 4 # collect my values - threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=nf) + threshold_holder = thresholds.ThresholdsAtlas.ffns(nf) as_FFNS_LO = StrongCoupling(alphas_ref, scale_ref, threshold_holder, order=0) my_vals = [] for Q2 in Q2s: @@ -247,7 +243,7 @@ def benchmark_apfel_exact(self): alphas_ref = 0.118 scale_ref = 90 ** 2 # collect my values - threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=3) + threshold_holder = thresholds.ThresholdsAtlas.ffns(3) # LHAPDF cache apfel_vals_dict = { 0: np.array( @@ -313,7 +309,7 @@ def benchmark_lhapdf_exact(self): alphas_ref = 0.118 scale_ref = 90 ** 2 # collect my values - threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=3) + threshold_holder = thresholds.ThresholdsAtlas.ffns(3) # LHAPDF cache lhapdf_vals_dict = { 0: np.array( @@ -391,9 +387,7 @@ def benchmark_lhapdf_zmvfns_lo(self): # Lambda2_3 = self._get_Lambda2_LO(as_FFNS_LO_4.a_s(m2c), m2c, 3) # collect my values - threshold_holder = thresholds.ThresholdsConfig( - scale_ref, "ZM-VFNS", threshold_list=threshold_list - ) + threshold_holder = thresholds.ThresholdsAtlas(threshold_list) as_VFNS_LO = StrongCoupling(alphas_ref, scale_ref, threshold_holder, order=0) my_vals = [] for Q2 in Q2s: diff --git a/tests/test_op_grid.py b/tests/test_op_grid.py index 4255a6727..cd2565d3a 100644 --- a/tests/test_op_grid.py +++ b/tests/test_op_grid.py @@ -10,7 +10,7 @@ import numpy as np import eko.interpolation as interpolation from eko.strong_coupling import StrongCoupling -from eko.thresholds import ThresholdsConfig +from eko.thresholds import ThresholdsAtlas from eko.operator.grid import OperatorGrid @@ -57,7 +57,7 @@ def _get_operator_grid(self, use_FFNS=True): basis_function_dispatcher = interpolation.InterpolatorDispatcher.from_dict( operators_card ) - threshold_holder = ThresholdsConfig.from_dict(theory_card) + threshold_holder = ThresholdsAtlas.from_dict(theory_card) a_s = StrongCoupling.from_dict(theory_card, threshold_holder) return OperatorGrid.from_dict( theory_card, @@ -83,7 +83,7 @@ def test_sanity(self): basis_function_dispatcher = interpolation.InterpolatorDispatcher.from_dict( operators_card ) - threshold_holder = ThresholdsConfig.from_dict(theory_card) + threshold_holder = ThresholdsAtlas.from_dict(theory_card) a_s = StrongCoupling.from_dict(theory_card, threshold_holder) theory_card.update({"ModEv": "wrong"}) OperatorGrid.from_dict( diff --git a/tests/test_operator.py b/tests/test_operator.py index aa6877021..541ce6066 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -5,7 +5,7 @@ from eko.operator import Operator, gamma_ns_fact, gamma_singlet_fact, quad_ker from eko.operator.grid import OperatorGrid -from eko.thresholds import ThresholdsConfig +from eko.thresholds import ThresholdsAtlas from eko.strong_coupling import StrongCoupling from eko.interpolation import InterpolatorDispatcher from eko import anomalous_dimensions as ad @@ -164,7 +164,7 @@ def test_compute(self, monkeypatch): g = OperatorGrid.from_dict( theory_card, operators_card, - ThresholdsConfig.from_dict(theory_card), + ThresholdsAtlas.from_dict(theory_card), StrongCoupling.from_dict(theory_card), InterpolatorDispatcher.from_dict(operators_card), ) diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 8467ac98d..ac73427e1 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -16,7 +16,7 @@ def test_init(self): alphas_ref = 0.118 scale_ref = 91.0 ** 2 nf = 4 - threshold_holder = thresholds.ThresholdsConfig(scale_ref, "FFNS", nf=nf) + threshold_holder = thresholds.ThresholdsAtlas.ffns(nf) # create sc = StrongCoupling(alphas_ref, scale_ref, threshold_holder) assert sc.q2_ref == scale_ref @@ -33,7 +33,7 @@ def test_init(self): NfFF=nf, Q0=2, ) - sc2 = StrongCoupling.from_dict(setup) + sc2 = StrongCoupling.from_dict(setup, threshold_holder) assert sc2.q2_ref == scale_ref assert sc2.as_ref == alphas_ref / 4.0 / np.pi @@ -57,23 +57,14 @@ def test_init(self): def test_ref(self): # prepare thresh_setups = [ - {"FNS": "FFNS", "NfFF": 3}, - {"FNS": "FFNS", "NfFF": 4}, - { - "FNS": "ZM-VFNS", - "mc": 2, - "mb": 4, - "mt": 175, - "kcThr": 1, - "kbThr": 1, - "ktThr": 1, - }, + (np.inf, np.inf, np.inf), + (0, np.inf, np.inf), + (2, 4, 175), ] alphas_ref = 0.118 scale_ref = 91.0 ** 2 for thresh_setup in thresh_setups: - thresh_setup["Q0"] = 1 - thresholds_conf = thresholds.ThresholdsConfig.from_dict(thresh_setup) + thresholds_conf = thresholds.ThresholdsAtlas(thresh_setup, 1) for order in [0, 1, 2]: for method in ["exact", "expanded"]: # create @@ -87,23 +78,14 @@ def test_ref(self): def test_exact_LO(self): # prepare thresh_setups = [ - {"FNS": "FFNS", "NfFF": 3}, - {"FNS": "FFNS", "NfFF": 4}, - { - "FNS": "ZM-VFNS", - "mc": 2, - "mb": 4, - "mt": 175, - "kcThr": 1, - "kbThr": 1, - "ktThr": 1, - }, + (np.inf, np.inf, np.inf), + (0, np.inf, np.inf), + (2, 4, 175), ] alphas_ref = 0.118 scale_ref = 91.0 ** 2 for thresh_setup in thresh_setups: - thresh_setup["Q0"] = 1 - thresholds_conf = thresholds.ThresholdsConfig.from_dict(thresh_setup) + thresholds_conf = thresholds.ThresholdsAtlas(thresh_setup, 1) # in LO expanded = exact sc_expanded = StrongCoupling( alphas_ref, scale_ref, thresholds_conf, 0, "expanded" diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py index 4c33b57a9..7ace35b33 100644 --- a/tests/test_thresholds.py +++ b/tests/test_thresholds.py @@ -5,104 +5,161 @@ import pytest import numpy as np -from eko.thresholds import ThresholdsConfig, Area +from eko.thresholds import ThresholdsAtlas, Area, PathSegment + + +class TestPathSegment: + def test_intersect_fully_inside(self): + a = Area(0, 1, 3) + fwd = PathSegment.intersect(0.2, 0.8, a) + assert fwd.q2_from == 0.2 + assert fwd.q2_to == 0.8 + assert fwd._area == a + bwd = PathSegment.intersect(0.8, 0.2, a) + assert bwd.q2_from == 0.8 + assert bwd.q2_to == 0.2 + assert bwd._area == a + + def test_intersect_too_low(self): + a = Area(0, 1, 3) + fwd = PathSegment.intersect(-0.2, 0.8, a) + assert fwd.q2_from == 0 + assert fwd.q2_to == 0.8 + assert fwd._area == a + bwd = PathSegment.intersect(0.8, -0.2, a) + assert bwd.q2_from == 0.8 + assert bwd.q2_to == 0.0 + assert bwd._area == a + + def test_intersect_no_overlap(self): + a = Area(0, 1, 3) + fwd = PathSegment.intersect(-0.2, 1.8, a) + assert fwd.q2_from == 0 + assert fwd.q2_to == 1 + assert fwd._area == a + bwd = PathSegment.intersect(1.8, -0.2, a) + assert bwd.q2_from == 1 + assert bwd.q2_to == 0.0 + assert bwd._area == a class TestThresholdsConfig: def test_init(self): - # at the moment it is fine to have ZM-VFNS with NfFF set - _tc = ThresholdsConfig(2, "ZM-VFNS", threshold_list=[1, 2], nf=2) - # errors - with pytest.raises(ValueError): - ThresholdsConfig(2, "FFNS") - with pytest.raises(ValueError): - ThresholdsConfig(2, "FFNS", nf=2, threshold_list=[1.0, 2.0]) - with pytest.raises(NotImplementedError): - ThresholdsConfig( - 2, "VFNS" - ) # I enforced ZM-VFNS -> "explicit is better then implicit" - with pytest.raises(ValueError): - ThresholdsConfig(2, "ZM-VFNS") - - def test_from_dict(self): - # here it is ok to have *all* keys set - tc_vfns = ThresholdsConfig.from_dict( - { - "FNS": "ZM-VFNS", - "NfFF": 3, - "Q0": np.sqrt(2), - "mc": 2, - "mb": 4, - "mt": 175, - "kcThr": 1, - "kbThr": 1, - "ktThr": 1, - } - ) - assert tc_vfns.scheme == "ZM-VFNS" - tc_ffns = ThresholdsConfig.from_dict( - {"FNS": "FFNS", "NfFF": 3, "Q0": np.sqrt(2)} - ) - assert tc_ffns.scheme == "FFNS" - - def test_nf_range_nf_ref(self): - # FFNS - for nf in [3, 4, 5, 6]: - tc_ffns = ThresholdsConfig(2, "FFNS", nf=nf) - assert list(tc_ffns.nf_range()) == [nf] - assert tc_ffns.nf_ref == nf - # VFNS - 1 threshold - for k, q2_ref in enumerate([2, 4]): - tc_vfns1 = ThresholdsConfig(q2_ref, "ZM-VFNS", threshold_list=[3]) - assert list(tc_vfns1.nf_range()) == [3, 4] - assert ( - tc_vfns1.nf_ref == 3 + k - ) # = 3 if below and 4 if above the one threshold - - def test_ffns(self): - # Check the setup for FFNS produces the correct thing - tholder = ThresholdsConfig(4.5, "FFNS", nf=4) - assert len(tholder._areas) == 1 # pylint: disable=protected-access - area_path = tholder.get_path_from_q2_ref(87.4) - assert len(area_path) == 1 - area = area_path[0] - assert area.q2_min == 0.0 - assert area.q2_max == np.inf - - def test_vnfs(self): - # Tests the setup for VFNS - th_list = [5, 50] - q2fin = 42 - for q2_ref in [2, 8, 61]: - tholder = ThresholdsConfig(q2_ref, "ZM-VFNS", threshold_list=th_list) - assert ( - len(tholder._areas) # pylint: disable=protected-access - == len(th_list) + 1 - ) - area_path = tholder.get_path_from_q2_ref(q2fin) - # The first area should contain q2_ref and the last q2fin - assert area_path[0](q2_ref) - assert area_path[-1](q2fin) - - -class TestAreas: - def _get_areas(self, q2_ref=42.0, nf=4): - q2_change = q2_ref * 0.99 - area_left = Area(0.0, q2_change, q2_ref, nf) - area_right = Area(q2_change, q2_ref * 2, q2_ref, nf) - area_far = Area(q2_ref * 4, q2_ref * 5, q2_ref, nf) - return area_left, area_right, area_far - - def test_order(self): - # Checks that areas work - area_left, area_right, area_far = self._get_areas() - assert area_left < area_right - assert area_right > area_left - assert area_far > area_left - - def test_q2_towards(self): - q2_ref = 42 - area_left, area_right, area_far = self._get_areas(q2_ref) - assert area_left.q2_towards(q2_ref) == area_left.q2_max - assert area_far.q2_towards(q2_ref) == area_far.q2_min - assert area_right.q2_towards(q2_ref) == q2_ref + # 3 thr + tc3 = ThresholdsAtlas([1, 2, 3]) + assert tc3.area_walls == [0, 1, 2, 3, np.inf] + assert tc3.areas[0].nf == 3 + assert tc3.areas[1].nf == 4 + # 2 thr + tc2 = ThresholdsAtlas([0, 2, 3]) + assert tc2.area_walls == [0, 0, 2, 3, np.inf] + assert tc2.areas[0].nf == 3 + assert tc2.areas[1].nf == 4 + + def test_path_3thr(self): + tc = ThresholdsAtlas([1, 2, 3], 0.5) + p1 = tc.path(0.7) + assert len(p1) == 1 + assert p1[0].q2_from == 0.5 + assert p1[0].q2_to == 0.7 + assert p1[0].nf == 3 + + +# class TestThresholdsConfig: +# def test_init(self): +# # at the moment it is fine to have ZM-VFNS with NfFF set +# _tc = ThresholdsAtlas(2, "ZM-VFNS", threshold_list=[1, 2], nf=2) +# # errors +# with pytest.raises(ValueError): +# ThresholdsAtlas(2, "FFNS") +# with pytest.raises(ValueError): +# ThresholdsAtlas(2, "FFNS", nf=2, threshold_list=[1.0, 2.0]) +# with pytest.raises(NotImplementedError): +# ThresholdsAtlas( +# 2, "VFNS" +# ) # I enforced ZM-VFNS -> "explicit is better then implicit" +# with pytest.raises(ValueError): +# ThresholdsAtlas(2, "ZM-VFNS") + +# def test_from_dict(self): +# # here it is ok to have *all* keys set +# tc_vfns = ThresholdsAtlas.from_dict( +# { +# "FNS": "ZM-VFNS", +# "NfFF": 3, +# "Q0": np.sqrt(2), +# "mc": 2, +# "mb": 4, +# "mt": 175, +# "kcThr": 1, +# "kbThr": 1, +# "ktThr": 1, +# } +# ) +# assert tc_vfns.scheme == "ZM-VFNS" +# tc_ffns = ThresholdsAtlas.from_dict( +# {"FNS": "FFNS", "NfFF": 3, "Q0": np.sqrt(2)} +# ) +# assert tc_ffns.scheme == "FFNS" + +# def test_nf_range_nf_ref(self): +# # FFNS +# for nf in [3, 4, 5, 6]: +# tc_ffns = ThresholdsAtlas(2, "FFNS", nf=nf) +# assert list(tc_ffns.nf_range()) == [nf] +# assert tc_ffns.nf_ref == nf +# # VFNS - 1 threshold +# for k, q2_ref in enumerate([2, 4]): +# tc_vfns1 = ThresholdsAtlas(q2_ref, "ZM-VFNS", threshold_list=[3]) +# assert list(tc_vfns1.nf_range()) == [3, 4] +# assert ( +# tc_vfns1.nf_ref == 3 + k +# ) # = 3 if below and 4 if above the one threshold + +# def test_ffns(self): +# # Check the setup for FFNS produces the correct thing +# tholder = ThresholdsAtlas(4.5, "FFNS", nf=4) +# assert len(tholder._areas) == 1 # pylint: disable=protected-access +# area_path = tholder.get_path_from_q2_ref(87.4) +# assert len(area_path) == 1 +# area = area_path[0] +# assert area.q2_min == 0.0 +# assert area.q2_max == np.inf + +# def test_vnfs(self): +# # Tests the setup for VFNS +# th_list = [5, 50] +# q2fin = 42 +# for q2_ref in [2, 8, 61]: +# tholder = ThresholdsAtlas(q2_ref, "ZM-VFNS", threshold_list=th_list) +# assert ( +# len(tholder._areas) # pylint: disable=protected-access +# == len(th_list) + 1 +# ) +# area_path = tholder.get_path_from_q2_ref(q2fin) +# # The first area should contain q2_ref and the last q2fin +# assert area_path[0](q2_ref) +# assert area_path[-1](q2fin) + + +# class TestAreas: +# def _get_areas(self, q2_ref=42.0, nf=4): +# q2_change = q2_ref * 0.99 +# area_left = Area(0.0, q2_change, q2_ref, nf) +# area_right = Area(q2_change, q2_ref * 2, q2_ref, nf) +# area_far = Area(q2_ref * 4, q2_ref * 5, q2_ref, nf) +# return area_left, area_right, area_far + +# def test_order(self): +# # Checks that areas work +# area_left, area_right, area_far = self._get_areas() +# assert area_left < area_right +# assert area_right > area_left +# assert area_far > area_left + +# def test_q2_towards(self): +# q2_ref = 42 +# area_left, area_right, area_far = self._get_areas(q2_ref) +# assert area_left.q2_towards(q2_ref) == area_left.q2_max +# assert area_far.q2_towards(q2_ref) == area_far.q2_min +# assert area_right.q2_towards(q2_ref) == q2_ref From 922333925e7bd54f3633d97f9a570530c4aa9958 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 30 Nov 2020 14:40:35 +0100 Subject: [PATCH 03/12] Fix tests --- tests/test_op_grid.py | 6 ++++++ tests/test_operator.py | 6 ++++++ tests/test_runner.py | 6 ++++++ 3 files changed, 18 insertions(+) diff --git a/tests/test_op_grid.py b/tests/test_op_grid.py index cd2565d3a..6094ad121 100644 --- a/tests/test_op_grid.py +++ b/tests/test_op_grid.py @@ -27,6 +27,12 @@ def _get_setup(self, use_FFNS): "FNS": "FFNS", "NfFF": 3, "IC": 1, + "mc": 1.0, + "mb": 4.75, + "mt": 173.0, + "kcThr": 0, + "kbThr": np.inf, + "ktThr": np.inf, } operators_card = { "Q2grid": [1, 10], diff --git a/tests/test_operator.py b/tests/test_operator.py index 541ce6066..4570c27bb 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -150,6 +150,12 @@ def test_compute(self, monkeypatch): "FNS": "FFNS", "NfFF": 3, "IC": 0, + "mc": 1.0, + "mb": 4.75, + "mt": 173.0, + "kcThr": 0, + "kbThr": np.inf, + "ktThr": np.inf, } operators_card = { "Q2grid": [1, 10], diff --git a/tests/test_runner.py b/tests/test_runner.py index e3f29233a..2d05139b9 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -18,6 +18,12 @@ def test_runner(): "NfFF": 3, "ModEv": "EXA", "IC": 0, + "mc": 1.0, + "mb": 4.75, + "mt": 173.0, + "kcThr": 0, + "kbThr": np.inf, + "ktThr": np.inf, } operators_card = { "Q2grid": [10, 100], From 03ef73ac5a869dd0baad52a4e545f0b41e312e85 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 30 Nov 2020 14:58:14 +0100 Subject: [PATCH 04/12] Improve and increase ThresholdAtlas tests --- tests/test_thresholds.py | 147 +++++++++++---------------------------- 1 file changed, 42 insertions(+), 105 deletions(-) diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py index 7ace35b33..0d8a6a9fa 100644 --- a/tests/test_thresholds.py +++ b/tests/test_thresholds.py @@ -14,33 +14,40 @@ def test_intersect_fully_inside(self): fwd = PathSegment.intersect(0.2, 0.8, a) assert fwd.q2_from == 0.2 assert fwd.q2_to == 0.8 - assert fwd._area == a + assert fwd._area == a # pylint: disable=protected-access bwd = PathSegment.intersect(0.8, 0.2, a) assert bwd.q2_from == 0.8 assert bwd.q2_to == 0.2 - assert bwd._area == a + assert bwd._area == a # pylint: disable=protected-access def test_intersect_too_low(self): a = Area(0, 1, 3) fwd = PathSegment.intersect(-0.2, 0.8, a) assert fwd.q2_from == 0 assert fwd.q2_to == 0.8 - assert fwd._area == a + assert fwd._area == a # pylint: disable=protected-access bwd = PathSegment.intersect(0.8, -0.2, a) assert bwd.q2_from == 0.8 assert bwd.q2_to == 0.0 - assert bwd._area == a + assert bwd._area == a # pylint: disable=protected-access def test_intersect_no_overlap(self): a = Area(0, 1, 3) fwd = PathSegment.intersect(-0.2, 1.8, a) assert fwd.q2_from == 0 assert fwd.q2_to == 1 - assert fwd._area == a + assert fwd._area == a # pylint: disable=protected-access bwd = PathSegment.intersect(1.8, -0.2, a) assert bwd.q2_from == 1 assert bwd.q2_to == 0.0 - assert bwd._area == a + assert bwd._area == a # pylint: disable=protected-access + + def test_print(self): + a = Area(0, 1, 3) + assert "nf=3" in str(a) + + p = PathSegment(0, 1, a) + assert "nf=3" in str(p) class TestThresholdsConfig: @@ -56,6 +63,31 @@ def test_init(self): assert tc2.areas[0].nf == 3 assert tc2.areas[1].nf == 4 + # errors + with pytest.raises(ValueError): + ThresholdsAtlas([1.0, 0.0]) + + def test_from_dict(self): + tc = ThresholdsAtlas.from_dict( + { + "mc": 1.0, + "mb": 4.0, + "mt": 100.0, + "kcThr": 1, + "kbThr": 2.0, + "ktThr": np.inf, + "Q0": 1.0, + } + ) + assert tc.area_walls[1:-1] == [1.0, 64.0, np.inf] + assert tc.q2_ref == 1.0 + + def test_ffns(self): + tc3 = ThresholdsAtlas.ffns(3) + assert tc3.area_walls == [0] + [np.inf] * 4 + tc4 = ThresholdsAtlas.ffns(4) + assert tc4.area_walls == [0] * 2 + [np.inf] * 3 + def test_path_3thr(self): tc = ThresholdsAtlas([1, 2, 3], 0.5) p1 = tc.path(0.7) @@ -64,102 +96,7 @@ def test_path_3thr(self): assert p1[0].q2_to == 0.7 assert p1[0].nf == 3 - -# class TestThresholdsConfig: -# def test_init(self): -# # at the moment it is fine to have ZM-VFNS with NfFF set -# _tc = ThresholdsAtlas(2, "ZM-VFNS", threshold_list=[1, 2], nf=2) -# # errors -# with pytest.raises(ValueError): -# ThresholdsAtlas(2, "FFNS") -# with pytest.raises(ValueError): -# ThresholdsAtlas(2, "FFNS", nf=2, threshold_list=[1.0, 2.0]) -# with pytest.raises(NotImplementedError): -# ThresholdsAtlas( -# 2, "VFNS" -# ) # I enforced ZM-VFNS -> "explicit is better then implicit" -# with pytest.raises(ValueError): -# ThresholdsAtlas(2, "ZM-VFNS") - -# def test_from_dict(self): -# # here it is ok to have *all* keys set -# tc_vfns = ThresholdsAtlas.from_dict( -# { -# "FNS": "ZM-VFNS", -# "NfFF": 3, -# "Q0": np.sqrt(2), -# "mc": 2, -# "mb": 4, -# "mt": 175, -# "kcThr": 1, -# "kbThr": 1, -# "ktThr": 1, -# } -# ) -# assert tc_vfns.scheme == "ZM-VFNS" -# tc_ffns = ThresholdsAtlas.from_dict( -# {"FNS": "FFNS", "NfFF": 3, "Q0": np.sqrt(2)} -# ) -# assert tc_ffns.scheme == "FFNS" - -# def test_nf_range_nf_ref(self): -# # FFNS -# for nf in [3, 4, 5, 6]: -# tc_ffns = ThresholdsAtlas(2, "FFNS", nf=nf) -# assert list(tc_ffns.nf_range()) == [nf] -# assert tc_ffns.nf_ref == nf -# # VFNS - 1 threshold -# for k, q2_ref in enumerate([2, 4]): -# tc_vfns1 = ThresholdsAtlas(q2_ref, "ZM-VFNS", threshold_list=[3]) -# assert list(tc_vfns1.nf_range()) == [3, 4] -# assert ( -# tc_vfns1.nf_ref == 3 + k -# ) # = 3 if below and 4 if above the one threshold - -# def test_ffns(self): -# # Check the setup for FFNS produces the correct thing -# tholder = ThresholdsAtlas(4.5, "FFNS", nf=4) -# assert len(tholder._areas) == 1 # pylint: disable=protected-access -# area_path = tholder.get_path_from_q2_ref(87.4) -# assert len(area_path) == 1 -# area = area_path[0] -# assert area.q2_min == 0.0 -# assert area.q2_max == np.inf - -# def test_vnfs(self): -# # Tests the setup for VFNS -# th_list = [5, 50] -# q2fin = 42 -# for q2_ref in [2, 8, 61]: -# tholder = ThresholdsAtlas(q2_ref, "ZM-VFNS", threshold_list=th_list) -# assert ( -# len(tholder._areas) # pylint: disable=protected-access -# == len(th_list) + 1 -# ) -# area_path = tholder.get_path_from_q2_ref(q2fin) -# # The first area should contain q2_ref and the last q2fin -# assert area_path[0](q2_ref) -# assert area_path[-1](q2fin) - - -# class TestAreas: -# def _get_areas(self, q2_ref=42.0, nf=4): -# q2_change = q2_ref * 0.99 -# area_left = Area(0.0, q2_change, q2_ref, nf) -# area_right = Area(q2_change, q2_ref * 2, q2_ref, nf) -# area_far = Area(q2_ref * 4, q2_ref * 5, q2_ref, nf) -# return area_left, area_right, area_far - -# def test_order(self): -# # Checks that areas work -# area_left, area_right, area_far = self._get_areas() -# assert area_left < area_right -# assert area_right > area_left -# assert area_far > area_left - -# def test_q2_towards(self): -# q2_ref = 42 -# area_left, area_right, area_far = self._get_areas(q2_ref) -# assert area_left.q2_towards(q2_ref) == area_left.q2_max -# assert area_far.q2_towards(q2_ref) == area_far.q2_min -# assert area_right.q2_towards(q2_ref) == q2_ref + p2 = tc.path(1.5, 2.5) + assert len(p2) == 2 + assert p2[0].nf == 5 + assert p2[1].nf == 4 From 7e1adba36bb08859970db2090f532de87927d492 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 30 Nov 2020 15:35:57 +0100 Subject: [PATCH 05/12] Filter empty PathSeg --- src/eko/operator/grid.py | 24 ++++++++---------------- src/eko/strong_coupling.py | 6 +----- src/eko/thresholds.py | 13 +++++++++---- tests/test_thresholds.py | 6 ++++++ 4 files changed, 24 insertions(+), 25 deletions(-) diff --git a/src/eko/operator/grid.py b/src/eko/operator/grid.py index 5fef1d58b..c8aa3e162 100644 --- a/src/eko/operator/grid.py +++ b/src/eko/operator/grid.py @@ -136,7 +136,7 @@ def from_dict( def _generate_thresholds_op(self, to_q2): """ - 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. @@ -154,13 +154,7 @@ def _generate_thresholds_op(self, to_q2): # The lists of areas as produced by the thresholds path = self.managers["thresholds_config"].path(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 seg in path[:-1]: - # q2_to = area.q2_ref - # TODO due to this continue we're actually seeing the area *one below* - if seg.q2_to == seg.q2_from: - continue new_op_key = (seg.q2_from, seg.q2_to) logger.info( "Evolution: Compute threshold operator from %e to %e", @@ -246,18 +240,17 @@ def _compute_raw_grid(self, q2grid): 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: @@ -279,8 +272,7 @@ def compute_q2grid(self, q2grid=None): def get_op_at_q2(self, qsq): """ - Given a value of q^2, returns the PhysicalOperator to get - to q^2 from q0^2 + Computes an single EKO. Parameters ---------- @@ -289,8 +281,8 @@ def get_op_at_q2(self, qsq): 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 """ # Check the path to q0 for this operator if qsq in self._op_grid: diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 1eb87bff9..f3ecf43b3 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -298,11 +298,7 @@ def a_s(self, scale_to, fact_scale=None): if fact_scale is None: fact_scale = scale_to for k, seg in enumerate(path): - q2_from = seg.q2_from - q2_to = seg.q2_to - if np.isclose(q2_from, q2_to): - continue - new_as = self._compute(final_as, seg.nf, q2_from, q2_to) + new_as = self._compute(final_as, seg.nf, seg.q2_from, seg.q2_to) # apply matching conditions: see hep-ph/9706430 # - if there is yet a step to go if k < len(path) - 1: diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index 844dca8dd..d653e95f1 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -72,8 +72,8 @@ def intersect(cls, q2_from, q2_to, area): class ThresholdsAtlas: """ - The threshold class holds information about the thresholds any - Q2 has to pass in order to get there from a given q2_ref and scheme. + Holds information about the thresholds any Q2 has to pass in order to get + there from a given q2_ref. Parameters ---------- @@ -115,13 +115,17 @@ def ffns(cls, nf, q2_ref=None): @classmethod def from_dict(cls, theory_card, prefix="k"): - """ + r""" Create the object from the run card. + The thresholds are computed by :math:`(m_q \cdot k_q^{Thr})`. + Parameters ---------- theory_card : dict run card with the keys given at the head of the :mod:`module ` + prefix : str + prefix for the ratio parameters Returns ------- @@ -167,7 +171,8 @@ def path(self, q2_to, q2_from=None): rc = -1 else: rc = 1 - return [ + path = [ PathSegment.intersect(q2_from, q2_to, self.areas[i - 1]) for i in range(ref_idx, target_idx + rc, rc) ] + return list(filter(lambda s: not np.allclose(s.q2_from, q2_to), path)) diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py index 0d8a6a9fa..24ecc20a4 100644 --- a/tests/test_thresholds.py +++ b/tests/test_thresholds.py @@ -100,3 +100,9 @@ def test_path_3thr(self): assert len(p2) == 2 assert p2[0].nf == 5 assert p2[1].nf == 4 + + def test_path_filter(self): + ta1 = ThresholdsAtlas([0, 2, 3], 0.5) + assert len(ta1.path(2.5)) == 2 + ta2 = ThresholdsAtlas([1, 2, 3], 0.5) + assert len(ta2.path(2.5)) == 3 From ee7d4a46a9bc1263ee3a65eeb43672130e81012a Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 30 Nov 2020 18:33:38 +0100 Subject: [PATCH 06/12] Apply PathSeg scheme in op.grid --- src/eko/operator/grid.py | 125 +++++++---------------------------- src/eko/operator/physical.py | 34 +++------- src/eko/runner.py | 2 +- src/eko/thresholds.py | 6 +- 4 files changed, 41 insertions(+), 126 deletions(-) diff --git a/src/eko/operator/grid.py b/src/eko/operator/grid.py index c8aa3e162..5fb72732a 100644 --- a/src/eko/operator/grid.py +++ b/src/eko/operator/grid.py @@ -134,7 +134,7 @@ 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. @@ -151,94 +151,26 @@ def _generate_thresholds_op(self, to_q2): to_q2: float value of q2 for which the OperatorGrid will need to pass thresholds """ - # The lists of areas as produced by the thresholds - path = self.managers["thresholds_config"].path(to_q2) # The base area is always that of the reference q + thr_ops = [] for seg in path[:-1]: - new_op_key = (seg.q2_from, seg.q2_to) - logger.info( - "Evolution: Compute threshold operator from %e to %e", - seg.q2_from, - seg.q2_to, - ) + new_op_key = seg.tuple if new_op_key not in self._threshold_operators: - # Compute the operator in place and store it + # Compute the operator and store it + logger.info( + "Evolution: Compute threshold operator from %e to %e", + seg.q2_from, + seg.q2_to, + ) 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 + thr_ops.append(self._threshold_operators[new_op_key]) + return thr_ops - 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 - path = self.managers["thresholds_config"].path(qsq) - # The last one is where q resides so it is not needed - op_list = [] - # Now loop over the areas to collect the necessary threshold operators - for seg in path[:-1]: - op_list.append(self._threshold_operators[(seg.q2_from, seg.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 - - 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 - """ - for q2 in q2grid: - path = self.managers["thresholds_config"].path(q2) - q2_from = path[-1].q2_from - nf = path[-1].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): + def compute(self, q2grid=None): """ Computes all ekos for the q2grid. @@ -258,25 +190,19 @@ def compute_q2grid(self, q2grid=None): # 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): """ Computes an single EKO. Parameters ---------- - qsq: float + q2: float Target value of q^2 Returns @@ -284,26 +210,23 @@ def get_op_at_q2(self, qsq): final_op: dict eko E(q^2 <- q_0^2) in flavor basis """ - # 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 = True # TODO #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() diff --git a/src/eko/operator/physical.py b/src/eko/operator/physical.py index 8c58c3231..a11ad9c3a 100644 --- a/src/eko/operator/physical.py +++ b/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,17 +64,16 @@ 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( @@ -87,18 +84,9 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, is_vfns, intrinsic_range=None) 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" - ) + # 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() diff --git a/src/eko/runner.py b/src/eko/runner.py index 39a540678..66c3f2030 100644 --- a/src/eko/runner.py +++ b/src/eko/runner.py @@ -104,7 +104,7 @@ def get_output(self): # add all operators Q2grid = {} self.out["pids"] = br.flavor_basis_pids - for final_scale, op in self.op_grid.compute_q2grid().items(): + for final_scale, op in self.op_grid.compute().items(): Q2grid[float(final_scale)] = op self.out["Q2grid"] = Q2grid return copy.deepcopy(self.out) diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index d653e95f1..9a36fa52d 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -59,6 +59,10 @@ def __init__(self, q2_from, q2_to, area): def nf(self): return self._area.nf + @property + def tuple(self): + return (self.q2_from, self.q2_to) + def __repr__(self): return f"PathSegment({self.q2_from} -> {self.q2_to}, nf={self.nf})" @@ -175,4 +179,4 @@ def path(self, q2_to, q2_from=None): PathSegment.intersect(q2_from, q2_to, self.areas[i - 1]) for i in range(ref_idx, target_idx + rc, rc) ] - return list(filter(lambda s: not np.allclose(s.q2_from, q2_to), path)) + return list(filter(lambda s: not np.allclose(s.q2_from, s.q2_to), path)) From 93c9996982c9fa7c5c3a007ad603e1b5784aaf0f Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 30 Nov 2020 18:42:47 +0100 Subject: [PATCH 07/12] Update docstrings --- src/eko/operator/grid.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/src/eko/operator/grid.py b/src/eko/operator/grid.py index 5fb72732a..40109c394 100644 --- a/src/eko/operator/grid.py +++ b/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. @@ -73,7 +72,6 @@ def __init__( strong_coupling=strong_coupling, interpol_dispatcher=interpol_dispatcher, ) - self._op_grid = {} self._threshold_operators = {} @classmethod @@ -88,15 +86,9 @@ 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.ThresholdsAtlas An instance of the ThresholdsAtlas class @@ -148,8 +140,8 @@ def get_threshold_operators(self, path): Parameters ---------- - to_q2: float - value of q2 for which the OperatorGrid will need to pass thresholds + path: list(PathSegment) + thresholds path """ # The base area is always that of the reference q thr_ops = [] @@ -158,9 +150,10 @@ def get_threshold_operators(self, path): if new_op_key not in self._threshold_operators: # Compute the operator and store it logger.info( - "Evolution: Compute threshold operator from %e to %e", + "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 @@ -208,7 +201,7 @@ def generate(self, q2): Returns ------- final_op: dict - eko E(q^2 <- q_0^2) in flavor basis + eko E(q^2 <- q_0^2) in flavor basis as numpy array """ # The lists of areas as produced by the thresholds path = self.managers["thresholds_config"].path(q2) From cacffb1c68bcfbbbf030f93882130ace2f3c1889 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Tue, 1 Dec 2020 17:17:53 +0100 Subject: [PATCH 08/12] Add ta.nf method --- src/eko/thresholds.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/eko/thresholds.py b/src/eko/thresholds.py index 9a36fa52d..ced879004 100644 --- a/src/eko/thresholds.py +++ b/src/eko/thresholds.py @@ -180,3 +180,20 @@ def path(self, q2_to, q2_from=None): for i in range(ref_idx, target_idx + rc, rc) ] return list(filter(lambda s: not np.allclose(s.q2_from, s.q2_to), path)) + + def nf(self, q2): + """ + Finds the number of flavor active at the given scale. + + Parameters + ---------- + q2 : float + reference scale + + Returns + ------- + nf : int + number of active flavors + """ + ref_idx = np.digitize(q2, self.area_walls) + return self.areas[ref_idx - 1].nf From 961b5dc7164674224a812f759326b7f49d011898 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 2 Dec 2020 10:54:22 +0100 Subject: [PATCH 09/12] Recover tests, fix small bug --- src/eko/operator/grid.py | 2 +- src/eko/operator/physical.py | 6 +- tests/test_op_grid.py | 59 ++++++++---------- tests/test_op_physical.py | 48 ++++---------- tests/test_operator.py | 118 +++++++++++++++++------------------ tests/test_thresholds.py | 8 +++ 6 files changed, 109 insertions(+), 132 deletions(-) diff --git a/src/eko/operator/grid.py b/src/eko/operator/grid.py index 40109c394..5b77947e9 100644 --- a/src/eko/operator/grid.py +++ b/src/eko/operator/grid.py @@ -153,7 +153,7 @@ def get_threshold_operators(self, path): "Threshold operator: %e -> %e, nf=%d", seg.q2_from, seg.q2_to, - seg.nf + seg.nf, ) op_th = Operator( self.config, self.managers, seg.nf, seg.q2_from, seg.q2_to diff --git a/src/eko/operator/physical.py b/src/eko/operator/physical.py index a11ad9c3a..edc45f146 100644 --- a/src/eko/operator/physical.py +++ b/src/eko/operator/physical.py @@ -79,11 +79,11 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range=None): 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 + 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 diff --git a/tests/test_op_grid.py b/tests/test_op_grid.py index 6094ad121..297292df6 100644 --- a/tests/test_op_grid.py +++ b/tests/test_op_grid.py @@ -27,12 +27,9 @@ def _get_setup(self, use_FFNS): "FNS": "FFNS", "NfFF": 3, "IC": 1, - "mc": 1.0, - "mb": 4.75, - "mt": 173.0, - "kcThr": 0, - "kbThr": np.inf, - "ktThr": np.inf, + "mc": 2.0, + "mb": 4.0, + "mt": 100.0, } operators_card = { "Q2grid": [1, 10], @@ -47,11 +44,11 @@ def _get_setup(self, use_FFNS): if use_FFNS: theory_card["FNS"] = "FFNS" theory_card["NfFF"] = 3 + theory_card["kcThr"] = np.inf + theory_card["kbThr"] = np.inf + theory_card["ktThr"] = np.inf else: theory_card["FNS"] = "ZM-VFNS" - theory_card["mc"] = 2 - theory_card["mb"] = 4 - theory_card["mt"] = 100 theory_card["kcThr"] = 1 theory_card["kbThr"] = 1 theory_card["ktThr"] = 1 @@ -78,12 +75,6 @@ def test_sanity(self): opgrid = self._get_operator_grid(False) # errors - with pytest.raises(ValueError): - opgrid.set_q2_limits(-1, 4) - with pytest.raises(ValueError): - opgrid.set_q2_limits(-1, -4) - with pytest.raises(ValueError): - opgrid.set_q2_limits(4, 1) with pytest.raises(ValueError): theory_card, operators_card = self._get_setup(True) basis_function_dispatcher = interpolation.InterpolatorDispatcher.from_dict( @@ -100,23 +91,23 @@ def test_sanity(self): basis_function_dispatcher, ) - def test_compute_q2grid(self): - opgrid = self._get_operator_grid() - # q2 has not be precomputed - but should work nevertheless - opgrid.get_op_at_q2(3) - # we can also pass a single number - opg = opgrid.compute_q2grid() - assert len(opg) == 2 - opg = opgrid.compute_q2grid(3) - assert len(opg) == 1 - # errors - with pytest.raises(ValueError): - bad_grid = [100, -6, 3] - _ = opgrid.compute_q2grid(bad_grid) + # def test_compute_q2grid(self): + # opgrid = self._get_operator_grid() + # # q2 has not be precomputed - but should work nevertheless + # opgrid.get_op_at_q2(3) + # # we can also pass a single number + # opg = opgrid.compute_q2grid() + # assert len(opg) == 2 + # opg = opgrid.compute_q2grid(3) + # assert len(opg) == 1 + # # errors + # with pytest.raises(ValueError): + # bad_grid = [100, -6, 3] + # _ = opgrid.compute_q2grid(bad_grid) - def test_grid_computation_VFNS(self): - """ Checks that the grid can be computed """ - opgrid = self._get_operator_grid(False) - qgrid_check = [3, 5] - operators = opgrid.compute_q2grid(qgrid_check) - assert len(operators) == len(qgrid_check) + # def test_grid_computation_VFNS(self): + # """ Checks that the grid can be computed """ + # opgrid = self._get_operator_grid(False) + # qgrid_check = [3, 5] + # operators = opgrid.compute_q2grid(qgrid_check) + # assert len(operators) == len(qgrid_check) diff --git a/tests/test_op_physical.py b/tests/test_op_physical.py index 5c3b1b2b8..cd677cdc0 100644 --- a/tests/test_op_physical.py +++ b/tests/test_op_physical.py @@ -106,50 +106,28 @@ def mk_op_members(shape=(2, 2)): return om -def get_ad_to_evol_map(nf, is_vfns, intrinsic_range=None): +def get_ad_to_evol_map(nf, intrinsic_range=None): oms = mk_op_members() - m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, is_vfns, intrinsic_range) + m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, intrinsic_range) return sorted(map(str, m.op_members.keys())) -def test_ad_to_evol_map_ffns(): +def test_ad_to_evol_map(): triv_ops = ("S.S", "S.g", "g.S", "g.g", "V.V", "V3.V3", "T3.T3", "V8.V8", "T8.T8") - # FFNS3 - assert sorted(triv_ops) == get_ad_to_evol_map(3, False) - # FFNS3 + IC - assert sorted([*triv_ops, "c+.c+", "c-.c-"]) == get_ad_to_evol_map(3, False, [4]) - # FFNS3 + IC + IB - assert sorted( - [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, False, [4, 5]) - # FFNS4 + IC + IB - with pytest.raises(ValueError): - get_ad_to_evol_map(4, False, [4, 5]) - # FFNS4 + IB - assert sorted( - [*triv_ops, "V15.V15", "T15.T15", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(4, False, [5]) - # FFNS6 - ks = get_ad_to_evol_map(6, False) - assert len(ks) == 4 + 1 + 2 * 5 - - -def test_ad_to_evol_map_vfns(): - triv_ops = ("S.S", "S.g", "g.S", "g.g", "V.V", "V3.V3", "T3.T3", "V8.V8", "T8.T8") - # VFNS, nf=3 patch - assert sorted([*triv_ops, "V15.V", "T15.S", "T15.g"]) == get_ad_to_evol_map(3, True) - # VFNS, nf=3 patch + IC + # nf=3 + assert sorted([*triv_ops, "V15.V", "T15.S", "T15.g"]) == get_ad_to_evol_map(3) + # nf=3 + IC assert sorted( [*triv_ops, "V15.V", "T15.S", "T15.g", "T15.c+", "V15.c-"] - ) == get_ad_to_evol_map(3, True, [4]) - # VFNS, nf=3 patch + IC + IB + ) == get_ad_to_evol_map(3, [4]) + # nf=3 + IC + IB assert sorted( [*triv_ops, "V15.V", "T15.S", "T15.g", "T15.c+", "V15.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, True, [4, 5]) - # VFNS, nf=4 patch + IC + IB + ) == get_ad_to_evol_map(3, [4, 5]) + # nf=4 + IC + IB ks = sorted( [*triv_ops, "V15.V15", "T15.T15", "V24.V", "T24.S", "T24.g", "T24.b+", "V24.b-"] ) - assert ks == get_ad_to_evol_map(4, True, [4, 5]) - # VFNS, nf=4 patch + IB - assert ks == get_ad_to_evol_map(4, True, [5]) + assert ks == get_ad_to_evol_map(4, [4, 5]) + # nf=4 + IB + assert ks == get_ad_to_evol_map(4, [5]) diff --git a/tests/test_operator.py b/tests/test_operator.py index 4570c27bb..f754cfe0e 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -137,62 +137,62 @@ def test_labels(self): ) assert sorted(o.labels()) == [] - def test_compute(self, monkeypatch): - # setup objs - theory_card = { - "alphas": 0.35, - "PTO": 0, - "ModEv": "TRN", - "XIF": 1.0, - "XIR": 1.0, - "Qref": np.sqrt(2), - "Q0": np.sqrt(2), - "FNS": "FFNS", - "NfFF": 3, - "IC": 0, - "mc": 1.0, - "mb": 4.75, - "mt": 173.0, - "kcThr": 0, - "kbThr": np.inf, - "ktThr": np.inf, - } - operators_card = { - "Q2grid": [1, 10], - "interpolation_xgrid": [0.1, 1.0], - "interpolation_polynomial_degree": 1, - "interpolation_is_log": True, - "debug_skip_singlet": True, - "debug_skip_non_singlet": False, - "ev_op_max_order": 1, - "ev_op_iterations": 1, - } - g = OperatorGrid.from_dict( - theory_card, - operators_card, - ThresholdsAtlas.from_dict(theory_card), - StrongCoupling.from_dict(theory_card), - InterpolatorDispatcher.from_dict(operators_card), - ) - g._compute_raw_grid([10]) # pylint: disable=protected-access - o = g._op_grid[10] # pylint: disable=protected-access - # fake quad - monkeypatch.setattr( - scipy.integrate, "quad", lambda *args, **kwargs: np.random.rand(2) - ) - # LO - o.compute() - assert "NS_m" in o.op_members - np.testing.assert_allclose( - o.op_members["NS_m"].value, o.op_members["NS_p"].value - ) - np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_p"].value - ) - # NLO - o.config["order"] = 1 - o.compute() - assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) - np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_m"].value - ) + # def test_compute(self, monkeypatch): + # # setup objs + # theory_card = { + # "alphas": 0.35, + # "PTO": 0, + # "ModEv": "TRN", + # "XIF": 1.0, + # "XIR": 1.0, + # "Qref": np.sqrt(2), + # "Q0": np.sqrt(2), + # "FNS": "FFNS", + # "NfFF": 3, + # "IC": 0, + # "mc": 1.0, + # "mb": 4.75, + # "mt": 173.0, + # "kcThr": 0, + # "kbThr": np.inf, + # "ktThr": np.inf, + # } + # operators_card = { + # "Q2grid": [1, 10], + # "interpolation_xgrid": [0.1, 1.0], + # "interpolation_polynomial_degree": 1, + # "interpolation_is_log": True, + # "debug_skip_singlet": True, + # "debug_skip_non_singlet": False, + # "ev_op_max_order": 1, + # "ev_op_iterations": 1, + # } + # g = OperatorGrid.from_dict( + # theory_card, + # operators_card, + # ThresholdsAtlas.from_dict(theory_card), + # StrongCoupling.from_dict(theory_card), + # InterpolatorDispatcher.from_dict(operators_card), + # ) + # ops = g.compute([10]) + # o = ops[10] + # # fake quad + # monkeypatch.setattr( + # scipy.integrate, "quad", lambda *args, **kwargs: np.random.rand(2) + # ) + # # LO + # o.compute() + # assert "NS_m" in o.op_members + # np.testing.assert_allclose( + # o.op_members["NS_m"].value, o.op_members["NS_p"].value + # ) + # np.testing.assert_allclose( + # o.op_members["NS_v"].value, o.op_members["NS_p"].value + # ) + # # NLO + # o.config["order"] = 1 + # o.compute() + # assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) + # np.testing.assert_allclose( + # o.op_members["NS_v"].value, o.op_members["NS_m"].value + # ) diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py index 24ecc20a4..ec6881c93 100644 --- a/tests/test_thresholds.py +++ b/tests/test_thresholds.py @@ -106,3 +106,11 @@ def test_path_filter(self): assert len(ta1.path(2.5)) == 2 ta2 = ThresholdsAtlas([1, 2, 3], 0.5) assert len(ta2.path(2.5)) == 3 + + def test_nf(self): + nf4 = ThresholdsAtlas.ffns(4) + for q2 in [1.0, 1e1, 1e2, 1e3, 1e4]: + assert nf4.nf(q2) == 4 + ta = ThresholdsAtlas([1, 2, 3], 0.5) + assert ta.nf(0.9) == 3 + assert ta.nf(1.1) == 4 From e32e4d304eaeafa8a6ef2b92dc9ce018a65d081f Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 2 Dec 2020 11:06:55 +0100 Subject: [PATCH 10/12] Increase test cov --- tests/test_thresholds.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/test_thresholds.py b/tests/test_thresholds.py index ec6881c93..7fab4feca 100644 --- a/tests/test_thresholds.py +++ b/tests/test_thresholds.py @@ -49,6 +49,15 @@ def test_print(self): p = PathSegment(0, 1, a) assert "nf=3" in str(p) + def test_tuple(self): + a = Area(0, 1, 3) + p = PathSegment(0, 1, a) + assert p.tuple == (0, 1) + # is hashable? + d = dict() + d[p.tuple] = 1 + assert d[p.tuple] == 1 + class TestThresholdsConfig: def test_init(self): From fc4dd934ce12e2f4b3417aa2cb7fc00985095f4b Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Wed, 2 Dec 2020 14:06:23 +0100 Subject: [PATCH 11/12] Increase a_s test cov --- tests/test_strong_coupling.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index ac73427e1..c018cccef 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -11,6 +11,23 @@ class TestStrongCoupling: + def test_from_dict(self): + d = { + "alphas": 0.118, + "Qref": 91., + "Q0": 1, + "PTO": 0, + "ModEv": "EXA", + "mc": 2., + "mb": 4., + "mt": 175., + "kcThr": 1., + "kbThr": 1., + "ktThr": 1., + } + sc = StrongCoupling.from_dict(d) + assert sc.a_s(d["Qref"]**2) == d["alphas"]/(4.*np.pi) + def test_init(self): # prepare alphas_ref = 0.118 From b576a16a993ef9e3e48333ccc32d105639551525 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 4 Dec 2020 11:24:11 +0100 Subject: [PATCH 12/12] Max test cov --- tests/test_op_grid.py | 36 +++++------ tests/test_operator.py | 117 +++++++++++++++++----------------- tests/test_strong_coupling.py | 16 ++--- 3 files changed, 81 insertions(+), 88 deletions(-) diff --git a/tests/test_op_grid.py b/tests/test_op_grid.py index 297292df6..18ec01f63 100644 --- a/tests/test_op_grid.py +++ b/tests/test_op_grid.py @@ -72,8 +72,6 @@ def _get_operator_grid(self, use_FFNS=True): def test_sanity(self): """ Sanity checks for the input""" - opgrid = self._get_operator_grid(False) - # errors with pytest.raises(ValueError): theory_card, operators_card = self._get_setup(True) @@ -91,23 +89,19 @@ def test_sanity(self): basis_function_dispatcher, ) - # def test_compute_q2grid(self): - # opgrid = self._get_operator_grid() - # # q2 has not be precomputed - but should work nevertheless - # opgrid.get_op_at_q2(3) - # # we can also pass a single number - # opg = opgrid.compute_q2grid() - # assert len(opg) == 2 - # opg = opgrid.compute_q2grid(3) - # assert len(opg) == 1 - # # errors - # with pytest.raises(ValueError): - # bad_grid = [100, -6, 3] - # _ = opgrid.compute_q2grid(bad_grid) + def test_compute_q2grid(self): + opgrid = self._get_operator_grid() + # q2 has not be precomputed - but should work nevertheless + opgrid.compute(3) + # we can also pass a single number + opg = opgrid.compute() + assert len(opg) == 2 + opg = opgrid.compute(3) + assert len(opg) == 1 - # def test_grid_computation_VFNS(self): - # """ Checks that the grid can be computed """ - # opgrid = self._get_operator_grid(False) - # qgrid_check = [3, 5] - # operators = opgrid.compute_q2grid(qgrid_check) - # assert len(operators) == len(qgrid_check) + def test_grid_computation_VFNS(self): + """ Checks that the grid can be computed """ + opgrid = self._get_operator_grid(False) + qgrid_check = [3, 5] + operators = opgrid.compute(qgrid_check) + assert len(operators) == len(qgrid_check) diff --git a/tests/test_operator.py b/tests/test_operator.py index f754cfe0e..0baaa54ae 100644 --- a/tests/test_operator.py +++ b/tests/test_operator.py @@ -137,62 +137,61 @@ def test_labels(self): ) assert sorted(o.labels()) == [] - # def test_compute(self, monkeypatch): - # # setup objs - # theory_card = { - # "alphas": 0.35, - # "PTO": 0, - # "ModEv": "TRN", - # "XIF": 1.0, - # "XIR": 1.0, - # "Qref": np.sqrt(2), - # "Q0": np.sqrt(2), - # "FNS": "FFNS", - # "NfFF": 3, - # "IC": 0, - # "mc": 1.0, - # "mb": 4.75, - # "mt": 173.0, - # "kcThr": 0, - # "kbThr": np.inf, - # "ktThr": np.inf, - # } - # operators_card = { - # "Q2grid": [1, 10], - # "interpolation_xgrid": [0.1, 1.0], - # "interpolation_polynomial_degree": 1, - # "interpolation_is_log": True, - # "debug_skip_singlet": True, - # "debug_skip_non_singlet": False, - # "ev_op_max_order": 1, - # "ev_op_iterations": 1, - # } - # g = OperatorGrid.from_dict( - # theory_card, - # operators_card, - # ThresholdsAtlas.from_dict(theory_card), - # StrongCoupling.from_dict(theory_card), - # InterpolatorDispatcher.from_dict(operators_card), - # ) - # ops = g.compute([10]) - # o = ops[10] - # # fake quad - # monkeypatch.setattr( - # scipy.integrate, "quad", lambda *args, **kwargs: np.random.rand(2) - # ) - # # LO - # o.compute() - # assert "NS_m" in o.op_members - # np.testing.assert_allclose( - # o.op_members["NS_m"].value, o.op_members["NS_p"].value - # ) - # np.testing.assert_allclose( - # o.op_members["NS_v"].value, o.op_members["NS_p"].value - # ) - # # NLO - # o.config["order"] = 1 - # o.compute() - # assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) - # np.testing.assert_allclose( - # o.op_members["NS_v"].value, o.op_members["NS_m"].value - # ) + def test_compute(self, monkeypatch): + # setup objs + theory_card = { + "alphas": 0.35, + "PTO": 0, + "ModEv": "TRN", + "XIF": 1.0, + "XIR": 1.0, + "Qref": np.sqrt(2), + "Q0": np.sqrt(2), + "FNS": "FFNS", + "NfFF": 3, + "IC": 0, + "mc": 1.0, + "mb": 4.75, + "mt": 173.0, + "kcThr": np.inf, + "kbThr": np.inf, + "ktThr": np.inf, + } + operators_card = { + "Q2grid": [1, 10], + "interpolation_xgrid": [0.1, 1.0], + "interpolation_polynomial_degree": 1, + "interpolation_is_log": True, + "debug_skip_singlet": True, + "debug_skip_non_singlet": False, + "ev_op_max_order": 1, + "ev_op_iterations": 1, + } + g = OperatorGrid.from_dict( + theory_card, + operators_card, + ThresholdsAtlas.from_dict(theory_card), + StrongCoupling.from_dict(theory_card), + InterpolatorDispatcher.from_dict(operators_card), + ) + o = Operator(g.config, g.managers, 3, 2, 10) + # fake quad + monkeypatch.setattr( + scipy.integrate, "quad", lambda *args, **kwargs: np.random.rand(2) + ) + # LO + o.compute() + assert "NS_m" in o.op_members + np.testing.assert_allclose( + o.op_members["NS_m"].value, o.op_members["NS_p"].value + ) + np.testing.assert_allclose( + o.op_members["NS_v"].value, o.op_members["NS_p"].value + ) + # NLO + o.config["order"] = 1 + o.compute() + assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) + np.testing.assert_allclose( + o.op_members["NS_v"].value, o.op_members["NS_m"].value + ) diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index c018cccef..8fd31533d 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -14,19 +14,19 @@ class TestStrongCoupling: def test_from_dict(self): d = { "alphas": 0.118, - "Qref": 91., + "Qref": 91.0, "Q0": 1, "PTO": 0, "ModEv": "EXA", - "mc": 2., - "mb": 4., - "mt": 175., - "kcThr": 1., - "kbThr": 1., - "ktThr": 1., + "mc": 2.0, + "mb": 4.0, + "mt": 175.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, } sc = StrongCoupling.from_dict(d) - assert sc.a_s(d["Qref"]**2) == d["alphas"]/(4.*np.pi) + assert sc.a_s(d["Qref"] ** 2) == d["alphas"] / (4.0 * np.pi) def test_init(self): # prepare