diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 67e901972..dfd0862f5 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -32,7 +32,7 @@ repos: - id: isort args: ["--profile", "black"] - repo: https://github.com/asottile/pyupgrade - rev: v3.16.0 + rev: v3.17.0 hooks: - id: pyupgrade - repo: https://github.com/pycqa/pydocstyle @@ -43,6 +43,13 @@ repos: args: ["--add-ignore=D107,D105"] additional_dependencies: - toml + - repo: https://github.com/pre-commit/mirrors-mypy + rev: v1.4.1 + hooks: + - id: mypy + additional_dependencies: [types-PyYAML] + pass_filenames: false + args: ["--ignore-missing-imports", "src/"] - repo: local hooks: - id: fmt @@ -53,6 +60,6 @@ repos: files: ^crates/.*\.rs$ args: [] - repo: https://github.com/pre-commit/pre-commit - rev: v3.7.1 + rev: v3.8.0 hooks: - id: validate_manifest diff --git a/crates/doc-header.html b/crates/doc-header.html new file mode 100644 index 000000000..3b2f9651b --- /dev/null +++ b/crates/doc-header.html @@ -0,0 +1,65 @@ + + + + diff --git a/crates/eko/Cargo.toml b/crates/eko/Cargo.toml index 0fc9b2ee7..0fffdea1a 100644 --- a/crates/eko/Cargo.toml +++ b/crates/eko/Cargo.toml @@ -13,7 +13,7 @@ rust-version.workspace = true version.workspace = true [package.metadata.docs.rs] -rustdoc-args = ["--html-in-header", "katex-header.html"] +rustdoc-args = ["--html-in-header", "doc-header.html"] [lib] name = "ekors" diff --git a/crates/eko/doc-header.html b/crates/eko/doc-header.html new file mode 120000 index 000000000..22282661d --- /dev/null +++ b/crates/eko/doc-header.html @@ -0,0 +1 @@ +../doc-header.html \ No newline at end of file diff --git a/crates/eko/katex-header.html b/crates/eko/katex-header.html deleted file mode 120000 index 810e68dcb..000000000 --- a/crates/eko/katex-header.html +++ /dev/null @@ -1 +0,0 @@ -../katex-header.html \ No newline at end of file diff --git a/crates/ekore/Cargo.toml b/crates/ekore/Cargo.toml index 86e9784eb..419bc126a 100644 --- a/crates/ekore/Cargo.toml +++ b/crates/ekore/Cargo.toml @@ -13,7 +13,7 @@ rust-version.workspace = true version.workspace = true [package.metadata.docs.rs] -rustdoc-args = ["--html-in-header", "katex-header.html"] +rustdoc-args = ["--html-in-header", "doc-header.html"] # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html diff --git a/crates/ekore/doc-header.html b/crates/ekore/doc-header.html new file mode 120000 index 000000000..22282661d --- /dev/null +++ b/crates/ekore/doc-header.html @@ -0,0 +1 @@ +../doc-header.html \ No newline at end of file diff --git a/crates/ekore/katex-header.html b/crates/ekore/katex-header.html deleted file mode 120000 index 810e68dcb..000000000 --- a/crates/ekore/katex-header.html +++ /dev/null @@ -1 +0,0 @@ -../katex-header.html \ No newline at end of file diff --git a/crates/katex-header.html b/crates/katex-header.html deleted file mode 100644 index ff018b9c3..000000000 --- a/crates/katex-header.html +++ /dev/null @@ -1,19 +0,0 @@ - - - - diff --git a/crates/parse-abbrev.py b/crates/parse-abbrev.py new file mode 100644 index 000000000..b928ed03c --- /dev/null +++ b/crates/parse-abbrev.py @@ -0,0 +1,20 @@ +"""Parse abbreviations from sphinx to Rust.""" + +import pathlib +import re + +SRC = ( + pathlib.Path(__file__).parents[1] + / "doc" + / "source" + / "shared" + / "abbreviations.rst" +) + +cnt = SRC.read_text("utf-8") +for el in cnt.split(".."): + test = re.match(r"\s*(.+)\s+replace::\n\s+:abbr:`(.+?)\((.+)\)`", el.strip()) + if test is None: + continue + # Print to terminal - the user can dump to the relevant file + print(f'"{test[2].strip()}": "{test[3].strip()}",') diff --git a/src/eko/couplings.py b/src/eko/couplings.py index 2b261afdb..a0826c09a 100644 --- a/src/eko/couplings.py +++ b/src/eko/couplings.py @@ -9,10 +9,11 @@ """ import logging -from typing import Iterable, List +from typing import Dict, Iterable, List, Tuple import numba as nb import numpy as np +import numpy.typing as npt import scipy from . import constants, matchings @@ -383,6 +384,10 @@ def couplings_expanded_fixed_alphaem(order, couplings_ref, nf, scale_from, scale return np.array([res_as, aem]) +_CouplingsCacheKey = Tuple[float, float, int, float, float] +"""Cache key containing (a0, a1, nf, scale_from, scale_to).""" + + class Couplings: r"""Compute the strong and electromagnetic coupling constants :math:`a_s, a_{em}`. @@ -480,7 +485,7 @@ def assert_positive(name, var): self.decoupled_running, ) # cache - self.cache = {} + self.cache: Dict[_CouplingsCacheKey, npt.NDArray] = {} @property def mu2_ref(self): diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 1c759c5c5..fe07ade9e 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -8,6 +8,7 @@ import os import time from multiprocessing import Pool +from typing import Dict, Tuple import numba as nb import numpy as np @@ -20,6 +21,8 @@ from .. import basis_rotation as br from .. import interpolation, mellin from .. import scale_variations as sv +from ..io.types import EvolutionMethod, OperatorLabel +from ..kernels import ev_method from ..kernels import non_singlet as ns from ..kernels import non_singlet_qed as qed_ns from ..kernels import singlet as s @@ -600,7 +603,11 @@ def quad_ker_qed( return ker -class Operator(sv.ModeMixin): +OpMembers = Dict[OperatorLabel, OpMember] +"""Map of all operators.""" + + +class Operator(sv.ScaleVariationModeMixin): """Internal representation of a single EKO. The actual matrices are computed upon calling :meth:`compute`. @@ -625,8 +632,8 @@ class Operator(sv.ModeMixin): log_label = "Evolution" # complete list of possible evolution operators labels - full_labels = br.full_labels - full_labels_qed = br.full_unified_labels + full_labels: Tuple[OperatorLabel, ...] = br.full_labels + full_labels_qed: Tuple[OperatorLabel, ...] = br.full_unified_labels def __init__( self, config, managers, segment: Segment, mellin_cut=5e-2, is_threshold=False @@ -639,9 +646,9 @@ def __init__( # TODO make 'cut' external parameter? self._mellin_cut = mellin_cut self.is_threshold = is_threshold - self.op_members = {} + self.op_members: OpMembers = {} self.order = tuple(config["order"]) - self.alphaem_running = self.managers["couplings"].alphaem_running + self.alphaem_running = self.managers.couplings.alphaem_running if self.log_label == "Evolution": self.a = self.compute_a() self.as_list, self.a_half_list = self.compute_aem_list() @@ -663,7 +670,7 @@ def xif2(self): @property def int_disp(self): """Return the interpolation dispatcher.""" - return self.managers["interpol_dispatcher"] + return self.managers.interpolator @property def grid_size(self): @@ -686,7 +693,7 @@ def mu2(self): def compute_a(self): """Return the computed values for :math:`a_s` and :math:`a_{em}`.""" - coupling = self.managers["couplings"] + coupling = self.managers.couplings a0 = coupling.a( self.mu2[0], nf_to=self.nf, @@ -722,7 +729,7 @@ def compute_aem_list(self): as_list = np.array([self.a_s[0], self.a_s[1]]) a_half = np.zeros((ev_op_iterations, 2)) else: - couplings = self.managers["couplings"] + couplings = self.managers.couplings mu2_steps = np.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations) mu2_l = mu2_steps[0] as_list = np.array( @@ -780,6 +787,11 @@ def labels(self): labels.extend(br.singlet_unified_labels) return labels + @property + def ev_method(self): + """Return the evolution method.""" + return ev_method(EvolutionMethod(self.config["method"])) + def quad_ker(self, label, logx, areas): """Return partially initialized integrand function. @@ -803,7 +815,7 @@ def quad_ker(self, label, logx, areas): order=self.order, mode0=label[0], mode1=label[1], - method=self.config["method"], + method=self.ev_method, is_log=self.int_disp.log, logx=logx, areas=areas, diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch index a3bc24ec1..995d77b6b 100644 --- a/src/eko/evolution_operator/__init__.py.patch +++ b/src/eko/evolution_operator/__init__.py.patch @@ -1,8 +1,8 @@ diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py -index 1c759c5c..5eb394d0 100644 +index fe07ade9..0f58c9e5 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py -@@ -3,15 +3,15 @@ r"""Contains the central operator classes. +@@ -3,16 +3,16 @@ r"""Contains the central operator classes. See :doc:`Operator overview `. """ @@ -11,6 +11,7 @@ index 1c759c5c..5eb394d0 100644 import os import time from multiprocessing import Pool + from typing import Dict, Tuple +import ekors import numba as nb @@ -20,7 +21,7 @@ index 1c759c5c..5eb394d0 100644 import ekore.anomalous_dimensions.polarized.space_like as ad_ps import ekore.anomalous_dimensions.unpolarized.space_like as ad_us -@@ -27,92 +27,10 @@ from ..kernels import singlet_qed as qed_s +@@ -30,92 +30,10 @@ from ..kernels import singlet_qed as qed_s from ..kernels import valence_qed as qed_v from ..matchings import Segment, lepton_number from ..member import OpMember @@ -114,7 +115,7 @@ index 1c759c5c..5eb394d0 100644 spec = [ ("is_singlet", nb.boolean), ("is_QEDsinglet", nb.boolean), -@@ -184,422 +102,6 @@ class QuadKerBase: +@@ -187,421 +105,6 @@ class QuadKerBase: return self.path.prefactor * pj * self.path.jac @@ -533,13 +534,12 @@ index 1c759c5c..5eb394d0 100644 - ) - return ker - -- - class Operator(sv.ModeMixin): - """Internal representation of a single EKO. -@@ -780,50 +282,6 @@ class Operator(sv.ModeMixin): - labels.extend(br.singlet_unified_labels) - return labels + OpMembers = Dict[OperatorLabel, OpMember] + """Map of all operators.""" +@@ -792,50 +295,6 @@ class Operator(sv.ScaleVariationModeMixin): + """Return the evolution method.""" + return ev_method(EvolutionMethod(self.config["method"])) - def quad_ker(self, label, logx, areas): - """Return partially initialized integrand function. @@ -564,7 +564,7 @@ index 1c759c5c..5eb394d0 100644 - order=self.order, - mode0=label[0], - mode1=label[1], -- method=self.config["method"], +- method=self.ev_method, - is_log=self.int_disp.log, - logx=logx, - areas=areas, @@ -588,7 +588,7 @@ index 1c759c5c..5eb394d0 100644 def initialize_op_members(self): """Init all operators with the identity or zeros.""" eye = OpMember( -@@ -846,10 +304,7 @@ class Operator(sv.ModeMixin): +@@ -858,10 +317,7 @@ class Operator(sv.ScaleVariationModeMixin): else: self.op_members[n] = zero.copy() @@ -600,7 +600,7 @@ index 1c759c5c..5eb394d0 100644 """Run the integration for each grid point. Parameters -@@ -864,18 +319,56 @@ class Operator(sv.ModeMixin): +@@ -876,18 +332,56 @@ class Operator(sv.ScaleVariationModeMixin): """ column = [] k, logx = log_grid diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index e64885418..6a9a3d9db 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -6,8 +6,8 @@ """ import logging -from dataclasses import astuple -from typing import Dict, List, Optional +from dataclasses import dataclass +from typing import Any, Dict, List, Optional import numpy as np import numpy.typing as npt @@ -18,9 +18,9 @@ from ..interpolation import InterpolatorDispatcher from ..io.runcards import Configs, Debug from ..io.types import EvolutionPoint as EPoint -from ..io.types import Order +from ..io.types import Order, SquaredScale from ..matchings import Atlas, Segment, flavor_shift, is_downward_path -from . import Operator, flavors, matching_condition, physical +from . import Operator, OpMembers, flavors, matching_condition, physical from .operator_matrix_element import OperatorMatrixElement logger = logging.getLogger(__name__) @@ -29,7 +29,16 @@ """In particular, only the ``operator`` and ``error`` fields are expected.""" -class OperatorGrid(sv.ModeMixin): +@dataclass(frozen=True) +class Managers: + """Set of steering objects.""" + + atlas: Atlas + couplings: Couplings + interpolator: InterpolatorDispatcher + + +class OperatorGrid(sv.ScaleVariationModeMixin): """Collection of evolution operators for several scales. The operator grid is the driver class of the evolution. @@ -64,7 +73,7 @@ def __init__( use_fhmruvv: bool, ): # check - config = {} + config: Dict[str, Any] = {} config["order"] = order config["intrinsic_range"] = intrinsic_flavors config["xif2"] = xif**2 @@ -95,13 +104,13 @@ def __init__( self.config = config self.q2_grid = mu2grid - self.managers = dict( - thresholds_config=atlas, + self.managers = Managers( + atlas=atlas, couplings=couplings, - interpol_dispatcher=interpol_dispatcher, + interpolator=interpol_dispatcher, ) - self._threshold_operators = {} - self._matching_operators = {} + self._threshold_operators: Dict[Segment, Operator] = {} + self._matching_operators: Dict[SquaredScale, OpMembers] = {} def get_threshold_operators(self, path: List[Segment]) -> List[Operator]: """Generate the threshold operators. @@ -123,7 +132,6 @@ def get_threshold_operators(self, path: List[Segment]) -> List[Operator]: is_downward = is_downward_path(path) shift = flavor_shift(is_downward) for seg in path[:-1]: - new_op_key = astuple(seg) kthr = self.config["thresholds_ratios"][seg.nf - shift] ome = OperatorMatrixElement( self.config, @@ -134,13 +142,13 @@ def get_threshold_operators(self, path: List[Segment]) -> List[Operator]: np.log(kthr), self.config["HQ"] == "MSBAR", ) - if new_op_key not in self._threshold_operators: + if seg not in self._threshold_operators: # Compute the operator and store it logger.info("Prepare threshold operator") op_th = Operator(self.config, self.managers, seg, is_threshold=True) op_th.compute() - self._threshold_operators[new_op_key] = op_th - thr_ops.append(self._threshold_operators[new_op_key]) + self._threshold_operators[seg] = op_th + thr_ops.append(self._threshold_operators[seg]) # Compute the matching conditions and store it if seg.target not in self._matching_operators: @@ -159,7 +167,7 @@ def generate(self, q2: EPoint) -> OpDict: """ # The lists of areas as produced by the thresholds - path = self.managers["thresholds_config"].path(q2) + path = self.managers.atlas.path(q2) # Prepare the path for the composition of the operator thr_ops = self.get_threshold_operators(path) # we start composing with the highest operator ... diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index 202fd18d9..fc29aafc6 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -1,6 +1,7 @@ """The |OME| for the non-trivial matching conditions in the |VFNS| evolution.""" import copy +import enum import functools import logging @@ -20,6 +21,33 @@ logger = logging.getLogger(__name__) +class MatchingMethods(enum.IntEnum): + """Enumerate matching methods.""" + + FORWARD = enum.auto() + BACKWARD_EXACT = enum.auto() + BACKWARD_EXPANDED = enum.auto() + + +def matching_method(s: InversionMethod) -> MatchingMethods: + """Return the matching method. + + Parameters + ---------- + s : + string representation + + Returns + ------- + i : + int representation + + """ + if s is not None: + return MatchingMethods["BACKWARD_" + s.value.upper()] + return MatchingMethods.FORWARD + + @nb.njit(cache=True) def build_ome(A, matching_order, a_s, backward_method): r"""Construct the matching expansion in :math:`a_s` with the appropriate method. @@ -32,7 +60,7 @@ def build_ome(A, matching_order, a_s, backward_method): perturbation matching order a_s : float strong coupling, needed only for the exact inverse - backward_method : InversionMethod or None + backward_method : MatchingMethods empty or method for inverting the matching condition (exact or expanded) Returns @@ -51,7 +79,7 @@ def build_ome(A, matching_order, a_s, backward_method): ome = np.eye(len(A[0]), dtype=np.complex_) A = A[:, :, :] A = np.ascontiguousarray(A) - if backward_method is InversionMethod.EXPANDED: + if backward_method is MatchingMethods.BACKWARD_EXPANDED: # expended inverse if matching_order[0] >= 1: ome -= a_s * A[0] @@ -68,7 +96,7 @@ def build_ome(A, matching_order, a_s, backward_method): if matching_order[0] >= 3: ome += a_s**3 * A[2] # need inverse exact ? so add the missing pieces - if backward_method is InversionMethod.EXACT: + if backward_method is MatchingMethods.BACKWARD_EXACT: ome = np.linalg.inv(ome) return ome @@ -199,7 +227,7 @@ class OperatorMatrixElement(Operator): log_label = "Matching" # complete list of possible matching operators labels - full_labels = [ + full_labels = ( *br.singlet_labels, (br.matching_hplus_pid, 21), (br.matching_hplus_pid, 100), @@ -210,13 +238,15 @@ class OperatorMatrixElement(Operator): (200, br.matching_hminus_pid), (br.matching_hminus_pid, 200), (br.matching_hminus_pid, br.matching_hminus_pid), - ] + ) # still valid in QED since Sdelta and Vdelta matchings are diagonal full_labels_qed = copy.deepcopy(full_labels) def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar): super().__init__(config, managers, Segment(q2, q2, nf)) - self.backward_method = config["backward_inversion"] if is_backward else None + self.backward_method = matching_method( + config["backward_inversion"] if is_backward else None + ) if is_backward: self.is_intrinsic = True else: @@ -241,7 +271,7 @@ def labels(self): logger.warning("%s: skipping non-singlet sector", self.log_label) else: labels.append((200, 200)) - if self.is_intrinsic or self.backward_method is not None: + if self.is_intrinsic or self.backward_method is not MatchingMethods.FORWARD: # intrinsic labels, which are not zero at NLO labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) # These contributions are always 0 for the moment @@ -257,7 +287,7 @@ def labels(self): (br.matching_hplus_pid, 100), ] ) - if self.is_intrinsic or self.backward_method is not None: + if self.is_intrinsic or self.backward_method is not MatchingMethods.FORWARD: labels.extend( [ (21, br.matching_hplus_pid), @@ -309,7 +339,7 @@ def a_s(self): Note that here you need to use :math:`a_s^{n_f+1}` """ - sc = self.managers["couplings"] + sc = self.managers.couplings return sc.a_s( self.q2_from * (self.xif2 if self.sv_mode == sv.Modes.exponentiated else 1.0), @@ -329,7 +359,7 @@ def compute(self): self.log_label, self.order[0], self.order[1], - self.backward_method, + MatchingMethods(self.backward_method).name, ) self.integrate() diff --git a/src/eko/io/inventory.py b/src/eko/io/inventory.py index 5bb98842f..9b25ad75d 100644 --- a/src/eko/io/inventory.py +++ b/src/eko/io/inventory.py @@ -3,7 +3,7 @@ import base64 from dataclasses import asdict, dataclass, field from pathlib import Path -from typing import Dict, Generic, Optional, Type, TypeVar +from typing import Dict, Generic, Literal, Optional, Type, TypeVar import yaml @@ -11,7 +11,7 @@ from .items import Header, Operator NBYTES = 8 -ENDIANNESS = "little" +ENDIANNESS: Literal["little", "big"] = "little" HEADER_EXT = ".yaml" ARRAY_EXT = [".npy", ".npz"] diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py index c05d150ec..0458b6a78 100644 --- a/src/eko/io/legacy.py +++ b/src/eko/io/legacy.py @@ -13,15 +13,23 @@ import numpy as np import yaml -from eko.interpolation import XGrid -from eko.io.runcards import flavored_mugrid -from eko.quantities.heavy_quarks import HeavyInfo, HeavyQuarkMasses, MatchingRatios - +from ..interpolation import XGrid +from ..io.runcards import flavored_mugrid +from ..quantities.heavy_quarks import ( + HeavyInfo, + HeavyQuarkMasses, + MatchingRatios, + QuarkMassScheme, +) from . import raw from .dictlike import DictLike from .struct import EKO, Operator from .types import EvolutionPoint as EPoint -from .types import RawCard +from .types import RawCard, ReferenceRunning + +_MC = 1.51 +_MB = 4.92 +_MT = 172.5 def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False): @@ -39,8 +47,8 @@ def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False): whether to load also errors (default ``False``) """ - with tempfile.TemporaryDirectory() as tmpdir: - tmpdir = pathlib.Path(tmpdir) + with tempfile.TemporaryDirectory() as tmpdirr: + tmpdir = pathlib.Path(tmpdirr) with tarfile.open(source, "r") as tar: raw.safe_extractall(tar, tmpdir) @@ -60,7 +68,7 @@ def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False): if op5 is None: op5 = metaold["mu2grid"] grid = op5to4( - flavored_mugrid(op5, theory.heavy.masses, theory.heavy.matching_ratios), arrays + flavored_mugrid(op5, [_MC, _MB, _MT], theory.heavy.matching_ratios), arrays ) with EKO.create(dest) as builder: @@ -93,10 +101,16 @@ def from_old(cls, old: RawCard): """Load from old metadata.""" heavy = HeavyInfo( num_flavs_init=4, - num_flavs_max_pdf=None, - intrinsic_flavors=None, - masses=HeavyQuarkMasses([1.51, 4.92, 172.5]), - masses_scheme=None, + num_flavs_max_pdf=5, + intrinsic_flavors=[], + masses=HeavyQuarkMasses( + [ + ReferenceRunning([_MC, np.inf]), + ReferenceRunning([_MB, np.inf]), + ReferenceRunning([_MT, np.inf]), + ] + ), + masses_scheme=QuarkMassScheme.POLE, matching_ratios=MatchingRatios([1.0, 1.0, 1.0]), ) return cls(heavy=heavy) @@ -125,7 +139,7 @@ def from_old(cls, old: RawCard): mu2list = old["mu2grid"] mu2grid = np.array(mu2list) evolgrid = flavored_mugrid( - np.sqrt(mu2grid).tolist(), [1.51, 4.92, 172.5], [1, 1, 1] + np.sqrt(mu2grid).tolist(), [_MC, _MB, _MT], [1, 1, 1] ) xgrid = XGrid(old["interpolation_xgrid"]) diff --git a/src/eko/io/struct.py b/src/eko/io/struct.py index ea0f4cadf..1520e1003 100644 --- a/src/eko/io/struct.py +++ b/src/eko/io/struct.py @@ -267,7 +267,7 @@ def approx( Raises ------ ValueError - if multiple values are find in the neighbourhood + if multiple values are found in the neighbourhood """ eps = np.array([ep_ for ep_ in self if ep_[1] == ep[1]]) @@ -275,7 +275,9 @@ def approx( close = eps[np.isclose(ep[0], mu2s, rtol=rtol, atol=atol)] if len(close) == 1: - return tuple(close[0]) + found = close[0] + assert isinstance(found[0], float) + return (found[0], int(found[1])) if len(close) == 0: return None raise ValueError(f"Multiple values of Q2 have been found close to {ep}") @@ -368,17 +370,17 @@ def open(cls, path: os.PathLike, mode="r"): raise ValueError(f"Unknown file mode: {mode}") tmpdir = pathlib.Path(tempfile.mkdtemp(prefix=TEMP_PREFIX)) - if load: - cls.load(path, tmpdir) - metadata = Metadata.load(tmpdir) - opened = cls( - **inventories(tmpdir, access), - metadata=metadata, - access=access, - ) - opened.operators.sync() - else: - opened = Builder(path=tmpdir, access=access) + if not load: + return Builder(path=tmpdir, access=access) + # load existing instead + cls.load(path, tmpdir) + metadata = Metadata.load(tmpdir) + opened: EKO = cls( + **inventories(tmpdir, access), + metadata=metadata, + access=access, + ) + opened.operators.sync() return opened diff --git a/src/eko/io/types.py b/src/eko/io/types.py index d8377a4e6..20ac52db5 100644 --- a/src/eko/io/types.py +++ b/src/eko/io/types.py @@ -1,8 +1,7 @@ """Common type definitions, only used for static analysis.""" import enum -import typing -from typing import Any, Dict, Generic, Tuple, TypeVar +from typing import Any, Dict, Generic, List, Tuple, TypeVar # Energy scales # ------------- @@ -23,8 +22,9 @@ Order = Tuple[int, int] FlavorsNumber = int FlavorIndex = int -IntrinsicFlavors = typing.List[FlavorIndex] -N3LOAdVariation = typing.Tuple[int, int, int, int] +IntrinsicFlavors = List[FlavorIndex] +N3LOAdVariation = Tuple[int, int, int, int] +OperatorLabel = Tuple[int, int] # Evolution coordinates # --------------------- diff --git a/src/eko/kernels/__init__.py b/src/eko/kernels/__init__.py index 7640727b5..594fbf1fa 100644 --- a/src/eko/kernels/__init__.py +++ b/src/eko/kernels/__init__.py @@ -1 +1,35 @@ """The solutions to the |DGLAP| equations.""" + +import enum + +from ..io.types import EvolutionMethod + + +class EvoMethods(enum.IntEnum): + """Enumerate evolution methods.""" + + ITERATE_EXACT = enum.auto() + ITERATE_EXPANDED = enum.auto() + PERTURBATIVE_EXACT = enum.auto() + PERTURBATIVE_EXPANDED = enum.auto() + TRUNCATED = enum.auto() + ORDERED_TRUNCATED = enum.auto() + DECOMPOSE_EXACT = enum.auto() + DECOMPOSE_EXPANDED = enum.auto() + + +def ev_method(s: EvolutionMethod) -> EvoMethods: + """Return the evolution method. + + Parameters + ---------- + s : + string representation + + Returns + ------- + i : + int representation + + """ + return EvoMethods[s.value.upper().replace("-", "_")] diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py index bd5bd10f0..15b913438 100644 --- a/src/eko/kernels/non_singlet.py +++ b/src/eko/kernels/non_singlet.py @@ -4,6 +4,7 @@ import numpy as np from .. import beta +from . import EvoMethods from . import as4_evolution_integrals as as4_ei from . import evolution_integrals as ei @@ -355,7 +356,7 @@ def dispatcher( ---------- order : tuple(int,int) perturbation order - method : str + method : int method gamma_ns : numpy.ndarray non-singlet anomalous dimensions @@ -378,38 +379,38 @@ def dispatcher( # use always exact in LO if order[0] == 1: return lo_exact(gamma_ns, a1, a0, betalist) - if method == "ordered-truncated": + if method is EvoMethods.ORDERED_TRUNCATED: return eko_ordered_truncated( gamma_ns, a1, a0, betalist, order, ev_op_iterations ) - if method == "truncated": + if method is EvoMethods.TRUNCATED: return eko_truncated(gamma_ns, a1, a0, betalist, order, ev_op_iterations) # NLO if order[0] == 2: if method in [ - "iterate-expanded", - "decompose-expanded", - "perturbative-expanded", + EvoMethods.ITERATE_EXPANDED, + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.PERTURBATIVE_EXPANDED, ]: return nlo_expanded(gamma_ns, a1, a0, betalist) - # if method in ["iterate-exact", "decompose-exact", "perturbative-exact"]: + # exact is left return nlo_exact(gamma_ns, a1, a0, betalist) # NNLO if order[0] == 3: if method in [ - "iterate-expanded", - "decompose-expanded", - "perturbative-expanded", + EvoMethods.ITERATE_EXPANDED, + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.PERTURBATIVE_EXPANDED, ]: return nnlo_expanded(gamma_ns, a1, a0, betalist) return nnlo_exact(gamma_ns, a1, a0, betalist) # N3LO if order[0] == 4: if method in [ - "iterate-expanded", - "decompose-expanded", - "perturbative-expanded", + EvoMethods.ITERATE_EXPANDED, + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.PERTURBATIVE_EXPANDED, ]: return n3lo_expanded(gamma_ns, a1, a0, betalist) return n3lo_exact(gamma_ns, a1, a0, betalist) diff --git a/src/eko/kernels/non_singlet_qed.py b/src/eko/kernels/non_singlet_qed.py index 7a2de816a..d81d5076b 100644 --- a/src/eko/kernels/non_singlet_qed.py +++ b/src/eko/kernels/non_singlet_qed.py @@ -146,7 +146,7 @@ def as4_exact(gamma_ns, a1, a0, beta): @nb.njit(cache=True) def dispatcher( order, - method, + _method, gamma_ns, as_list, aem_half, @@ -164,7 +164,7 @@ def dispatcher( ---------- order : tuple(int,int) perturbation order - method : str + method : int method gamma_ns : numpy.ndarray non-singlet anomalous dimensions diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index 5983c0c67..e637af046 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -6,6 +6,7 @@ from ekore import anomalous_dimensions as ad from .. import beta +from . import EvoMethods from . import as4_evolution_integrals as as4_ei from . import evolution_integrals as ei @@ -579,7 +580,7 @@ def dispatcher( # pylint: disable=too-many-return-statements ---------- order : tuple(int,int) perturbative order - method : str + method : int method gamma_singlet : numpy.ndarray singlet anomalous dimensions matrices @@ -609,9 +610,9 @@ def dispatcher( # pylint: disable=too-many-return-statements return lo_exact(gamma_singlet, a1, a0, betalist) # Common method for NLO and NNLO - if method in ["iterate-exact", "iterate-expanded"]: + if method in [EvoMethods.ITERATE_EXACT, EvoMethods.ITERATE_EXPANDED]: return eko_iterate(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) - if method == "perturbative-exact": + if method is EvoMethods.PERTURBATIVE_EXACT: return eko_perturbative( gamma_singlet, a1, @@ -622,7 +623,7 @@ def dispatcher( # pylint: disable=too-many-return-statements ev_op_max_order, True, ) - if method == "perturbative-expanded": + if method is EvoMethods.PERTURBATIVE_EXPANDED: return eko_perturbative( gamma_singlet, a1, @@ -633,19 +634,19 @@ def dispatcher( # pylint: disable=too-many-return-statements ev_op_max_order, False, ) - if method in ["truncated", "ordered-truncated"]: + if method in [EvoMethods.TRUNCATED, EvoMethods.ORDERED_TRUNCATED]: return eko_truncated(gamma_singlet, a1, a0, betalist, order, ev_op_iterations) # These methods are scattered for nlo and nnlo - if method == "decompose-exact": + if method is EvoMethods.DECOMPOSE_EXACT: if order[0] == 2: return nlo_decompose_exact(gamma_singlet, a1, a0, betalist) - elif order[0] == 3: + if order[0] == 3: return nnlo_decompose_exact(gamma_singlet, a1, a0, betalist) return n3lo_decompose_exact(gamma_singlet, a1, a0, nf) - if method == "decompose-expanded": + if method is EvoMethods.DECOMPOSE_EXPANDED: if order[0] == 2: return nlo_decompose_expanded(gamma_singlet, a1, a0, betalist) - elif order[0] == 3: + if order[0] == 3: return nnlo_decompose_expanded(gamma_singlet, a1, a0, betalist) return n3lo_decompose_expanded(gamma_singlet, a1, a0, nf) raise NotImplementedError("Selected method is not implemented") diff --git a/src/eko/kernels/singlet_qed.py b/src/eko/kernels/singlet_qed.py index c13dee95c..d63a2f1ac 100644 --- a/src/eko/kernels/singlet_qed.py +++ b/src/eko/kernels/singlet_qed.py @@ -6,6 +6,7 @@ from ekore import anomalous_dimensions as ad from .. import beta +from . import EvoMethods @nb.njit(cache=True) @@ -66,7 +67,7 @@ def dispatcher( a_half, nf, ev_op_iterations, - ev_op_max_order, + _ev_op_max_order, ): """Determine used kernel and call it. @@ -74,7 +75,7 @@ def dispatcher( ---------- order : tuple(int,int) perturbative order - method : str + method : int method gamma_singlet : numpy.ndarray singlet anomalous dimensions matrices @@ -96,7 +97,7 @@ def dispatcher( e_s : numpy.ndarray singlet EKO """ - if method in ["iterate-exact", "iterate-expanded"]: + if method is EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_singlet, as_list, a_half, nf, order, ev_op_iterations, 4 ) diff --git a/src/eko/kernels/valence_qed.py b/src/eko/kernels/valence_qed.py index 0d5b747c3..186c70219 100644 --- a/src/eko/kernels/valence_qed.py +++ b/src/eko/kernels/valence_qed.py @@ -2,6 +2,7 @@ import numba as nb +from . import EvoMethods from .singlet_qed import eko_iterate @@ -14,7 +15,7 @@ def dispatcher( a_half, nf, ev_op_iterations, - ev_op_max_order, + _ev_op_max_order, ): """ Determine used kernel and call it. @@ -23,7 +24,7 @@ def dispatcher( ---------- order : tuple(int,int) perturbative order - method : str + method : int method gamma_singlet : numpy.ndarray singlet anomalous dimensions matrices @@ -45,7 +46,7 @@ def dispatcher( e_v : numpy.ndarray singlet EKO """ - if method in ["iterate-exact", "iterate-expanded"]: + if method is EvoMethods.ITERATE_EXACT: return eko_iterate( gamma_valence, as_list, a_half, nf, order, ev_op_iterations, 2 ) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 26fb9874e..de605764b 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -396,7 +396,7 @@ def sc(thr_masses): heavy_quarks = quark_names[3:] hq_idxs = np.arange(0, 3) if nf_ref > 4: - heavy_quarks = reversed(heavy_quarks) + heavy_quarks = "".join(reversed(heavy_quarks)) hq_idxs = reversed(hq_idxs) # loop on heavy quarks and compute the msbar masses diff --git a/src/eko/quantities/couplings.py b/src/eko/quantities/couplings.py index aec1d083f..b223ce2a7 100644 --- a/src/eko/quantities/couplings.py +++ b/src/eko/quantities/couplings.py @@ -2,7 +2,6 @@ import dataclasses import enum -from typing import Optional from ..io.dictlike import DictLike from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, Scalar @@ -23,7 +22,7 @@ class CouplingsInfo(DictLike): alphaem: Coupling scale: LinearScale max_num_flavs: FlavorsNumber - num_flavs_ref: Optional[FlavorsNumber] + num_flavs_ref: FlavorsNumber r"""Number of active flavors at strong coupling reference scale. I.e. :math:`n_{f,\text{ref}}(\mu^2_{\text{ref}})`, formerly called diff --git a/src/eko/runner/operators.py b/src/eko/runner/operators.py index 07ac34093..1efa8bb07 100644 --- a/src/eko/runner/operators.py +++ b/src/eko/runner/operators.py @@ -17,7 +17,9 @@ def retrieve( elements = [] for head in headers: inv = parts if isinstance(head, Evolution) else parts_matching - elements.append(inv[head]) + op = inv[head] + assert op is not None + elements.append(op) return elements diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py index f9379bf76..798353f71 100644 --- a/src/eko/runner/parts.py +++ b/src/eko/runner/parts.py @@ -16,13 +16,14 @@ from ..evolution_operator import matching_condition from ..evolution_operator import operator_matrix_element as ome from ..evolution_operator import physical +from ..evolution_operator.grid import Managers from ..io import EKO from ..io.items import Evolution, Matching, Operator from ..quantities.heavy_quarks import QuarkMassScheme from . import commons -def managers(eko: EKO) -> dict: +def managers(eko: EKO) -> Managers: """Collect managers for operator computation. .. todo:: @@ -31,10 +32,10 @@ def managers(eko: EKO) -> dict: """ tcard = eko.theory_card ocard = eko.operator_card - return dict( - thresholds_config=commons.atlas(tcard, ocard), + return Managers( + atlas=commons.atlas(tcard, ocard), couplings=commons.couplings(tcard, ocard), - interpol_dispatcher=commons.interpolator(ocard), + interpolator=commons.interpolator(ocard), ) diff --git a/src/eko/scale_variations/__init__.py b/src/eko/scale_variations/__init__.py index 4af80592b..8d3f01556 100644 --- a/src/eko/scale_variations/__init__.py +++ b/src/eko/scale_variations/__init__.py @@ -5,30 +5,32 @@ """ import enum +from typing import Any, Dict +from ..io.types import ScaleVariationsMethod from . import expanded, exponentiated class Modes(enum.IntEnum): - """Enumerate scale Variation modes.""" + """Enumerate scale variation modes.""" unvaried = enum.auto() exponentiated = enum.auto() expanded = enum.auto() -def sv_mode(s): +def sv_mode(s: ScaleVariationsMethod) -> Modes: """Return the scale variation mode. Parameters ---------- - s : str + s : string representation Returns ------- - enum.IntEnum - enum representation + i : + int representation """ if s is not None: @@ -36,10 +38,12 @@ def sv_mode(s): return Modes.unvaried -class ModeMixin: +class ScaleVariationModeMixin: """Mixin to cast scale variation mode.""" + config: Dict[str, Any] + @property - def sv_mode(self): + def sv_mode(self) -> Modes: """Return the scale variation mode.""" return sv_mode(self.config["ModSV"]) diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py index 08be115c4..bf9f2083d 100644 --- a/src/ekobox/apply.py +++ b/src/ekobox/apply.py @@ -1,10 +1,14 @@ """Apply operator evolution to PDF set.""" +from dataclasses import dataclass +from typing import Dict, Optional, Union + import numpy as np from eko import basis_rotation as br from eko import interpolation from eko.io import EKO +from eko.io.types import EvolutionPoint def apply_pdf( @@ -49,6 +53,17 @@ def apply_pdf( CONTRACTION = "ajbk,bk" +_PdfLabel = Union[int, str] +"""PDF identifier: either PID or label.""" + + +@dataclass +class PdfResult: + """Helper class to collect PDF results.""" + + pdfs: Dict[_PdfLabel, float] + errors: Optional[Dict[_PdfLabel, float]] = None + def apply_pdf_flavor( eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None, labels=None @@ -85,35 +100,34 @@ def apply_pdf_flavor( ) # build output - out_grid = {} + out_grid: Dict[EvolutionPoint, PdfResult] = {} for ep, elem in eko.items(): pdf_final = np.einsum(CONTRACTION, elem.operator, pdfs, optimize="optimal") if elem.error is not None: error_final = np.einsum(CONTRACTION, elem.error, pdfs, optimize="optimal") else: error_final = None - out_grid[ep] = { - "pdfs": dict(zip(br.flavor_basis_pids, pdf_final)), - "errors": None, - } + out_grid[ep] = PdfResult(dict(zip(br.flavor_basis_pids, pdf_final))) if error_final is not None: - out_grid[ep]["errors"] = dict(zip(br.flavor_basis_pids, error_final)) + out_grid[ep].errors = dict(zip(br.flavor_basis_pids, error_final)) qed = eko.theory_card.order[1] > 0 # rotate to evolution basis if flavor_rotation is not None: for q2, op in out_grid.items(): pdf = flavor_rotation @ np.array( - [op["pdfs"][pid] for pid in br.flavor_basis_pids] + [op.pdfs[pid] for pid in br.flavor_basis_pids] ) - if labels is None: - labels = list(range(flavor_rotation.shape[0])) - op["pdfs"] = dict(zip(labels, pdf)) - if op["errors"] is not None: + if not qed: + evol_basis = br.evol_basis + else: + evol_basis = br.unified_evol_basis + op.pdfs = dict(zip(evol_basis, pdf)) + if op.errors is not None: errors = flavor_rotation @ np.array( - [op["errors"][pid] for pid in br.flavor_basis_pids] + [op.errors[pid] for pid in br.flavor_basis_pids] ) - op["errors"] = dict(zip(labels, errors)) + op.errors = dict(zip(evol_basis, errors)) # rotate/interpolate to target grid if targetgrid is not None: @@ -125,11 +139,12 @@ def apply_pdf_flavor( rot = b.get_interpolation(targetgrid) for evpdf in out_grid.values(): - for pdf_label in evpdf["pdfs"]: - evpdf["pdfs"][pdf_label] = np.matmul(rot, evpdf["pdfs"][pdf_label]) - if evpdf["errors"] is not None: - evpdf["errors"][pdf_label] = np.matmul( - rot, evpdf["errors"][pdf_label] - ) - - return out_grid + for pdf_label in evpdf.pdfs: + evpdf.pdfs[pdf_label] = np.matmul(rot, evpdf.pdfs[pdf_label]) + if evpdf.errors is not None: + evpdf.errors[pdf_label] = np.matmul(rot, evpdf.errors[pdf_label]) + # cast back to be backward compatible + real_out_grid = {} + for ep, res in out_grid.items(): + real_out_grid[ep] = {"pdfs": res.pdfs, "errors": res.errors} + return real_out_grid diff --git a/src/ekobox/cli/inspect.py b/src/ekobox/cli/inspect.py index 5787c5f1e..c355fce61 100644 --- a/src/ekobox/cli/inspect.py +++ b/src/ekobox/cli/inspect.py @@ -28,7 +28,7 @@ def subcommand(ctx, path: pathlib.Path): @click.pass_obj def sub_mu2(operator: EKO): """Check operator's mu2grid.""" - rich.print_json(data=operator.mu2grid.tolist()) + rich.print_json(data=operator.mu2grid) @subcommand.command("cards") diff --git a/src/ekobox/cli/run.py b/src/ekobox/cli/run.py index 9a91e36ab..44cfe216e 100644 --- a/src/ekobox/cli/run.py +++ b/src/ekobox/cli/run.py @@ -53,11 +53,11 @@ def subcommand(paths: Sequence[pathlib.Path]): else: output = operator.parent / OUTPUT - theory = yaml.safe_load(theory.read_text(encoding="utf-8")) - if "order" in theory: - theory = TheoryCard.from_dict(theory) - operator = yaml.safe_load(operator.read_text(encoding="utf-8")) - if "configs" in operator: - operator = OperatorCard.from_dict(operator) - - eko.solve(theory, operator, path=output) + tc = yaml.safe_load(theory.read_text(encoding="utf-8")) + if "order" in tc: + tc = TheoryCard.from_dict(tc) + oc = yaml.safe_load(operator.read_text(encoding="utf-8")) + if "configs" in oc: + oc = OperatorCard.from_dict(oc) + + eko.solve(tc, oc, path=output) diff --git a/src/ekobox/cli/runcards.py b/src/ekobox/cli/runcards.py index 020a58917..227e012f4 100644 --- a/src/ekobox/cli/runcards.py +++ b/src/ekobox/cli/runcards.py @@ -3,6 +3,8 @@ import logging import pathlib +import numpy as np + from .. import cards from . import library as lib from .base import command @@ -35,6 +37,6 @@ def sub_example(destination: pathlib.Path): cards.dump(theory.raw, path=destination / "theory.yaml") operator = cards.example.operator() operator.mu0 = 1.65 - operator.mu2grid = [1e5] + operator.mugrid = [(np.sqrt(1e5), 5)] cards.dump(operator.raw, path=destination / "operator.yaml") _logger.info(f"Runcards generated to '{destination}'") diff --git a/src/ekobox/evol_pdf.py b/src/ekobox/evol_pdf.py index ddd80f75a..9da6dda65 100644 --- a/src/ekobox/evol_pdf.py +++ b/src/ekobox/evol_pdf.py @@ -1,13 +1,15 @@ """Tools to evolve actual PDFs.""" import pathlib -from collections import defaultdict + +import numpy as np from eko import basis_rotation as br from eko.io import EKO from eko.runner import managed from . import apply, genpdf, info_file +from .utils import regroup_evolgrid DEFAULT_NAME = "eko.tar" @@ -46,6 +48,18 @@ def evolve_pdfs( info_update : dict dict of info to add or update to default info file """ + # separate by nf the evolgrid (and order per nf/q) + q2block_per_nf = regroup_evolgrid(operators_card.mugrid) + + # check we have disjoint scale ranges + nfs = list(q2block_per_nf.keys()) + for j in range(len(nfs) - 1): + # equal points are allowed by LHAPDF + if q2block_per_nf[nfs[j]][-1] > q2block_per_nf[nfs[j + 1]][0]: + raise ValueError( + f"Last scale point for nf={nfs[j]} is bigger than first in nf={nfs[j+1]}" + ) + # update op and th cards if path is not None: eko_path = pathlib.Path(path) @@ -59,12 +73,16 @@ def evolve_pdfs( # apply PDF to eko evolved_PDF_list = [] + q2block_per_nf = {} with EKO.read(eko_path) as eko_output: for initial_PDF in initial_PDF_list: evolved_PDF_list.append( apply.apply_pdf(eko_output, initial_PDF, targetgrid) ) + # separate by nf the evolgrid (and order per nf/q) + q2block_per_nf = regroup_evolgrid(eko_output.evolgrid) # pylint: disable=E1101 + # update info file if targetgrid is None: targetgrid = operators_card.xgrid @@ -79,9 +97,6 @@ def evolve_pdfs( info_update=info_update, ) - # separate by nf the evolgrid (and order per nf/q) - q2block_per_nf = regroup_evolgrid(eko_output.evolgrid) - # write all replicas all_member_blocks = [] for evolved_PDF in evolved_PDF_list: @@ -95,14 +110,6 @@ def evolve_pdfs( genpdf.install_pdf(name) -def regroup_evolgrid(evolgrid: list): - """Split evolution points by nf and sort by scale.""" - by_nf = defaultdict(list) - for q, nf in sorted(evolgrid, key=lambda ep: ep[1]): - by_nf[nf].append(q) - return {nf: sorted(qs) for nf, qs in by_nf.items()} - - def collect_blocks(evolved_PDF: dict, q2block_per_nf: dict, xgrid: list): """Collect all LHAPDF blocks for a given replica. @@ -129,7 +136,7 @@ def pdf_xq2(pid, x, Q2): pdf_xq2, xgrid=xgrid, sorted_q2grid=q2grid, - pids=br.flavor_basis_pids, + pids=np.array(br.flavor_basis_pids), ) all_blocks.append(block) return all_blocks diff --git a/src/ekobox/genpdf/__init__.py b/src/ekobox/genpdf/__init__.py index 4666891bc..d6c427989 100644 --- a/src/ekobox/genpdf/__init__.py +++ b/src/ekobox/genpdf/__init__.py @@ -17,7 +17,7 @@ def take_data( - parent_pdf_set: Optional[Union[str, dict]] = None, + parent_pdf_set=None, members: bool = False, xgrid: Optional[List[float]] = None, evolgrid: Optional[List[EPoint]] = None, @@ -63,7 +63,7 @@ def take_data( all_blocks.append( [ generate_block( - toylh.xfxQ2, xgrid, sorted_q2grid, br.flavor_basis_pids + toylh.xfxQ2, xgrid, sorted_q2grid, list(br.flavor_basis_pids) ) ] ) @@ -87,7 +87,7 @@ def take_data( ), xgrid, sorted_q2grid, - br.flavor_basis_pids, + list(br.flavor_basis_pids), ) ] ) diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index 26e3cb170..56439599d 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -2,6 +2,7 @@ import copy import math +from typing import Any, Dict import numpy as np @@ -9,6 +10,7 @@ from eko.io.runcards import OperatorCard, TheoryCard from .genpdf import load +from .utils import regroup_evolgrid def build( @@ -16,19 +18,19 @@ def build( operators_card: OperatorCard, num_members: int, info_update: dict, -): - """Generate a lhapdf info file from theory and operators card. +) -> dict: + """Generate a lhapdf info file. Parameters ---------- - theory_card : dict + theory_card : theory card - operators_card : dict - operators_card - num_members : int + operators_card : + operators card + num_members : number of pdf set members - info_update : dict - info to update + info_update : + additional info to update Returns ------- @@ -56,8 +58,35 @@ def build( template_info["MCharm"] = theory_card.heavy.masses.c.value template_info["MBottom"] = theory_card.heavy.masses.b.value template_info["MTop"] = theory_card.heavy.masses.t.value + # dump alphas + template_info.update(build_alphas(theory_card, operators_card)) + return template_info + + +def build_alphas( + theory_card: TheoryCard, + operators_card: OperatorCard, +) -> dict: + """Generate a couplings section of lhapdf info file. + + Parameters + ---------- + theory_card : dict + theory card + operators_card : dict + operators card + + Returns + ------- + dict + info file section in lhapdf format + """ + # start with meta stuff + template_info: Dict[str, Any] = {} template_info["AlphaS_MZ"] = theory_card.couplings.alphas template_info["AlphaS_OrderQCD"] = theory_card.order[0] - 1 + # prepare + evolgrid = regroup_evolgrid(operators_card.mugrid) evmod = couplings.couplings_mod_ev(operators_card.configs.evolution_method) quark_masses = [(x.value) ** 2 for x in theory_card.heavy.masses] sc = couplings.Couplings( @@ -68,19 +97,14 @@ def build( hqm_scheme=theory_card.heavy.masses_scheme, thresholds_ratios=np.power(list(iter(theory_card.heavy.matching_ratios)), 2.0), ) - alphas_values = np.array( - [ - 4.0 - * np.pi - * sc.a_s( - muf2, - ) - for muf2 in operators_card.mu2grid - ], - dtype=float, - ) - template_info["AlphaS_Vals"] = alphas_values.tolist() - template_info["AlphaS_Qs"] = np.array( - [mu for mu, _ in operators_card.mugrid] - ).tolist() + # add actual values + alphas_values = [] + alphas_qs = [] + for nf, mus in evolgrid.items(): + for mu in mus: + alphas_values.append(float(4.0 * np.pi * sc.a_s(mu * mu, nf_to=nf))) + alphas_qs.append(mu) + + template_info["AlphaS_Vals"] = alphas_values + template_info["AlphaS_Qs"] = alphas_qs return template_info diff --git a/src/ekobox/utils.py b/src/ekobox/utils.py index 0047d0872..091cde15e 100644 --- a/src/ekobox/utils.py +++ b/src/ekobox/utils.py @@ -1,6 +1,7 @@ """Generic utilities to work with EKOs.""" import os +from collections import defaultdict from typing import Optional import numpy as np @@ -78,3 +79,11 @@ def ekos_product( if path is not None: final_eko.close() + + +def regroup_evolgrid(evolgrid: list): + """Split evolution points by nf and sort by scale.""" + by_nf = defaultdict(list) + for q, nf in sorted(evolgrid, key=lambda ep: ep[1]): + by_nf[nf].append(q) + return {nf: sorted(qs) for nf, qs in by_nf.items()} diff --git a/src/ekomark/data/operators.py b/src/ekomark/data/operators.py index ad928574d..03ad4d349 100644 --- a/src/ekomark/data/operators.py +++ b/src/ekomark/data/operators.py @@ -15,7 +15,7 @@ ev_op_max_order=10, ev_op_iterations=10, backward_inversion="expanded", - n_integration_cores=0, + n_integration_cores=1, debug_skip_non_singlet=False, debug_skip_singlet=False, mugrid=[10], diff --git a/src/ekomark/plots.py b/src/ekomark/plots.py index a228cefcd..5e2bee420 100644 --- a/src/ekomark/plots.py +++ b/src/ekomark/plots.py @@ -2,6 +2,7 @@ import io import pprint +from typing import Dict, List import matplotlib.pyplot as plt import numpy as np @@ -251,7 +252,7 @@ def save_operators_to_pdf( # plot the operators # it's necessary to reshuffle the eko output - for mu2, op in me.items(): + for ep, op in me.items(): if rotate_to_evolution_basis: if not qed: op = manipulate.to_evol(op, source=True, target=True) @@ -264,8 +265,8 @@ def save_operators_to_pdf( for label_out, res, res_err in zip(ops_names, results, errors): if label_out in skip_pdfs: continue - new_op = {} - new_op_err = {} + new_op: Dict[int, List[float]] = {} + new_op_err: Dict[int, List[float]] = {} # loop on xgrid point for j in range(len(me.xgrid)): # loop on pid in @@ -281,17 +282,10 @@ def save_operators_to_pdf( for label_in in ops_names: if label_in in skip_pdfs: continue - lab_in = label_in - lab_out = label_out - if rotate_to_evolution_basis: - index_in = br.evol_basis_pids.index(label_in) - index_out = br.evol_basis_pids.index(label_out) - lab_in = br.evol_basis[index_in] - lab_out = br.evol_basis[index_out] fig = None try: fig = plot_operator( - f"Operator ({lab_in};{lab_out}) µ_F^2 = {mu2} GeV^2", + f"Operator ({label_in};{label_out}) µ_F^2 = {ep[0]} GeV^2, nf = {ep[1]}", new_op[label_in], new_op_err[label_in], ) diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index 597f645f8..ed003f925 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -1,4 +1,5 @@ import os +from dataclasses import dataclass import numpy as np import pytest @@ -11,6 +12,7 @@ from eko.evolution_operator import Operator, quad_ker from eko.interpolation import InterpolatorDispatcher from eko.io.runcards import OperatorCard, ScaleVariationsMethod, TheoryCard +from eko.kernels import EvoMethods from eko.kernels import non_singlet as ns from eko.kernels import non_singlet_qed as qed_ns from eko.kernels import singlet as s @@ -25,7 +27,7 @@ def test_quad_ker_errors(): order=(1, 0), mode0=mode0, mode1=0, - method="", + method="iterate-exact", is_log=True, logx=np.log(0.1), areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], @@ -60,15 +62,29 @@ def test_quad_ker(monkeypatch): monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0) monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2)) params = [ - ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0), - ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.123, 1.0), - ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0), - ((1, 0), 100, 100, "", 0.123, 1.0), - ((1, 0), 100, 21, "", 0.0, 0.0), - ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0), - ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0), - ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0), - ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0), + ((1, 0), br.non_singlet_pids_map["ns+"], 0, EvoMethods.ITERATE_EXACT, 0.0, 0.0), + ( + (1, 0), + br.non_singlet_pids_map["ns+"], + 0, + EvoMethods.ITERATE_EXACT, + 0.123, + 1.0, + ), + ( + (3, 1), + br.non_singlet_pids_map["ns+u"], + 0, + EvoMethods.ITERATE_EXACT, + 0.0, + 0.0, + ), + ((1, 0), 100, 100, EvoMethods.ITERATE_EXACT, 0.123, 1.0), + ((1, 0), 100, 21, EvoMethods.ITERATE_EXACT, 0.0, 0.0), + ((1, 1), 100, 100, EvoMethods.ITERATE_EXACT, 0.123, 1.0), + ((1, 1), 100, 21, EvoMethods.ITERATE_EXACT, 0.123, 0.0), + ((1, 1), 10200, 10200, EvoMethods.ITERATE_EXACT, 0.123, 1.0), + ((1, 1), 10200, 10204, EvoMethods.ITERATE_EXACT, 0.123, 0.0), ] for order, mode0, mode1, method, logx, res in params: for is_log in [True, False]: @@ -107,7 +123,7 @@ def test_quad_ker(monkeypatch): order=(1, 0), mode0=label[0], mode1=label[1], - method="", + method=EvoMethods.ITERATE_EXACT, is_log=True, logx=0.123, areas=np.zeros(3), @@ -143,7 +159,7 @@ def test_quad_ker(monkeypatch): order=(1, 1), mode0=label[0], mode1=label[1], - method="iterate-exact", + method=EvoMethods.ITERATE_EXACT, is_log=True, logx=0.123, areas=np.zeros(3), @@ -171,7 +187,7 @@ def test_quad_ker(monkeypatch): order=(1, 0), mode0=br.non_singlet_pids_map["ns+"], mode1=0, - method="", + method=EvoMethods.ITERATE_EXACT, is_log=True, logx=0.0, areas=np.zeros(3), @@ -212,7 +228,12 @@ def compute(self, a_ref, nf, scale_from, scale_to): return a_ref -fake_managers = {"couplings": FakeCoupling()} +@dataclass(frozen=True) +class FakeManagers: + couplings: FakeCoupling + + +fake_managers = FakeManagers(couplings=FakeCoupling()) class TestOperator: @@ -436,7 +457,7 @@ def quad_ker_pegasus( order = (2, 0) mode0 = br.non_singlet_pids_map["ns+"] mode1 = 0 - method = "" + method = EvoMethods.ITERATE_EXACT logxs = np.log(int_disp.xgrid.raw) a1 = 1 a0 = 2 diff --git a/tests/eko/evolution_operator/test_ome.py b/tests/eko/evolution_operator/test_ome.py index b844290c9..553e5823c 100644 --- a/tests/eko/evolution_operator/test_ome.py +++ b/tests/eko/evolution_operator/test_ome.py @@ -8,6 +8,7 @@ from eko import interpolation, mellin from eko import scale_variations as sv from eko.evolution_operator.operator_matrix_element import ( + MatchingMethods, OperatorMatrixElement, build_ome, quad_ker, @@ -33,7 +34,7 @@ def test_build_ome_as(): aS = A_singlet((o, 0), N, nf, L, is_msbar) for a in [aNS, aS]: - for method in [None, InversionMethod.EXPANDED, InversionMethod.EXACT]: + for method in MatchingMethods: dim = len(a[0]) if o != 1: assert len(a) == o @@ -53,7 +54,7 @@ def test_build_ome_nlo(): aNSi = A_non_singlet((1, 0), N, nf, L) aSi = A_singlet((1, 0), N, nf, L, is_msbar) for a in [aNSi, aSi]: - for method in [None, InversionMethod.EXPANDED, InversionMethod.EXACT]: + for method in MatchingMethods: dim = len(a[0]) # hh assert a[0, -1, -1] != 0.0 @@ -86,7 +87,7 @@ def test_quad_ker_errors(): is_log=True, logx=0.123, areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], - backward_method=None, + backward_method=MatchingMethods.FORWARD, a_s=0.0, nf=3, L=0.0, @@ -119,7 +120,7 @@ def test_quad_ker(monkeypatch): is_log=is_log, logx=0.123, areas=np.zeros(3), - backward_method=None, + backward_method=MatchingMethods.FORWARD, a_s=0.0, nf=3, L=0.0, @@ -138,7 +139,7 @@ def test_quad_ker(monkeypatch): is_log=is_log, logx=0.123, areas=np.zeros(3), - backward_method=None, + backward_method=MatchingMethods.FORWARD, a_s=0.0, nf=3, L=0.0, @@ -157,7 +158,7 @@ def test_quad_ker(monkeypatch): is_log=is_log, logx=0.0, areas=np.zeros(3), - backward_method=None, + backward_method=MatchingMethods.FORWARD, a_s=0.0, nf=3, L=0.0, @@ -191,7 +192,7 @@ def test_quad_ker(monkeypatch): is_log=True, logx=0.123, areas=np.zeros(3), - backward_method=InversionMethod.EXPANDED, + backward_method=MatchingMethods.BACKWARD_EXPANDED, a_s=0.0, nf=3, L=0.0, @@ -228,7 +229,7 @@ def test_quad_ker(monkeypatch): is_log=True, logx=0.123, areas=np.zeros(3), - backward_method=InversionMethod.EXACT, + backward_method=MatchingMethods.BACKWARD_EXACT, a_s=0.0, nf=3, L=0.0, @@ -252,7 +253,7 @@ def test_quad_ker(monkeypatch): is_log=True, logx=0.0, areas=np.array([0.01, 0.1, 1.0]), - backward_method=None, + backward_method=MatchingMethods.FORWARD, a_s=0.0, nf=3, L=0.0, diff --git a/tests/eko/kernels/test_init.py b/tests/eko/kernels/test_init.py new file mode 100644 index 000000000..3717de275 --- /dev/null +++ b/tests/eko/kernels/test_init.py @@ -0,0 +1,21 @@ +from eko.io.types import EvolutionMethod +from eko.kernels import EvoMethods, ev_method + + +def test_ev_method(): + methods = { + "iterate-expanded": EvoMethods.ITERATE_EXPANDED, + "decompose-expanded": EvoMethods.DECOMPOSE_EXPANDED, + "perturbative-expanded": EvoMethods.PERTURBATIVE_EXPANDED, + "truncated": EvoMethods.TRUNCATED, + "ordered-truncated": EvoMethods.ORDERED_TRUNCATED, + "iterate-exact": EvoMethods.ITERATE_EXACT, + "decompose-exact": EvoMethods.DECOMPOSE_EXACT, + "perturbative-exact": EvoMethods.PERTURBATIVE_EXACT, + } + assert len(methods.keys()) == len(EvolutionMethod) + assert len(methods.keys()) == len(EvoMethods) + for s, i in methods.items(): + j = ev_method(EvolutionMethod(s)) + assert j == i + assert isinstance(j, int) diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index bb18c3e17..9fce3097b 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -3,6 +3,7 @@ from eko import basis_rotation as br from eko.couplings import Couplings +from eko.kernels import EvoMethods from eko.kernels import non_singlet_qed as ns from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo from eko.quantities.heavy_quarks import QuarkMassScheme @@ -14,7 +15,7 @@ # "perturbative-expanded", # "truncated", # "ordered-truncated", - "iterate-exact", + EvoMethods.ITERATE_EXACT, # "decompose-exact", # "perturbative-exact", ] diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py index b953c053c..2fee4e17f 100644 --- a/tests/eko/kernels/test_kernels_QEDsinglet.py +++ b/tests/eko/kernels/test_kernels_QEDsinglet.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from eko.kernels import EvoMethods from eko.kernels import singlet_qed as s from ekore.anomalous_dimensions.unpolarized import space_like as ad @@ -10,7 +11,7 @@ # "perturbative-expanded", # "truncated", # "ordered-truncated", - "iterate-exact", + EvoMethods.ITERATE_EXACT, # "decompose-exact", # "perturbative-exact", ] @@ -125,5 +126,12 @@ def test_zero_true_gamma(monkeypatch): def test_error(): with pytest.raises(NotImplementedError): s.dispatcher( - (3, 2), "AAA", np.random.rand(4, 3, 2, 2), [0.2, 0.1], [0.01], 3, 10, 10 + (3, 2), + "iterate-exact", + np.random.rand(4, 3, 2, 2), + [0.2, 0.1], + [0.01], + 3, + 10, + 10, ) diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py index 29f0f3d7a..4b12a492b 100644 --- a/tests/eko/kernels/test_kernels_QEDvalence.py +++ b/tests/eko/kernels/test_kernels_QEDvalence.py @@ -1,6 +1,7 @@ import numpy as np import pytest +from eko.kernels import EvoMethods from eko.kernels import valence_qed as val from ekore import anomalous_dimensions @@ -10,13 +11,13 @@ # "perturbative-expanded", # "truncated", # "ordered-truncated", - "iterate-exact", + EvoMethods.ITERATE_EXACT, # "decompose-exact", # "perturbative-exact", ] -def test_zero(monkeypatch): +def test_zero(): """No evolution results in exp(0)""" nf = 3 ev_op_iterations = 2 @@ -57,7 +58,7 @@ def test_zero(monkeypatch): ) -def test_zero_true_gamma(monkeypatch): +def test_zero_true_gamma(): """No evolution results in exp(0)""" nf = 3 ev_op_iterations = 2 @@ -101,5 +102,12 @@ def test_zero_true_gamma(monkeypatch): def test_error(): with pytest.raises(NotImplementedError): val.dispatcher( - (3, 2), "AAA", np.random.rand(4, 3, 2, 2), [0.2, 0.1], [0.01], 3, 10, 10 + (3, 2), + "iterate-exact", + np.random.rand(4, 3, 2, 2), + [0.2, 0.1], + [0.01], + 3, + 10, + 10, ) diff --git a/tests/eko/kernels/test_ns.py b/tests/eko/kernels/test_ns.py index f78282142..70e7faf8d 100644 --- a/tests/eko/kernels/test_ns.py +++ b/tests/eko/kernels/test_ns.py @@ -5,19 +5,9 @@ import ekore.anomalous_dimensions.unpolarized.space_like as ad from eko import beta +from eko.kernels import EvoMethods from eko.kernels import non_singlet as ns -methods = [ - "iterate-expanded", - "decompose-expanded", - "perturbative-expanded", - "truncated", - "ordered-truncated", - "iterate-exact", - "decompose-exact", - "perturbative-exact", -] - def test_zero(): """No evolution results in exp(0)""" @@ -25,7 +15,7 @@ def test_zero(): ev_op_iterations = 2 gamma_ns = np.array([1 + 0.0j, 1 + 0j, 1 + 0j, 1 + 0j]) for order in [1, 2, 3, 4]: - for method in methods: + for method in EvoMethods: np.testing.assert_allclose( ns.dispatcher( (order, 0), method, gamma_ns, 1.0, 1.0, nf, ev_op_iterations @@ -54,7 +44,7 @@ def test_ode_lo(): a0 = 0.3 for a1 in [0.1, 0.2]: r = a1 * gamma_ns / (beta.beta_qcd((2, 0), nf) * a1**2) - for method in methods: + for method in EvoMethods: rhs = r * ns.dispatcher( (1, 0), method, gamma_ns, a1, a0, nf, ev_op_iterations ) @@ -91,7 +81,7 @@ def test_ode_nlo(): r = (a1 * gamma_ns[0] + a1**2 * gamma_ns[1]) / ( beta.beta_qcd((2, 0), nf) * a1**2 + beta.beta_qcd((3, 0), nf) * a1**3 ) - for method in ["iterate-exact"]: + for method in [EvoMethods.ITERATE_EXACT]: rhs = r * ns.dispatcher( (2, 0), method, gamma_ns, a1, a0, nf, ev_op_iterations ) @@ -130,7 +120,7 @@ def test_ode_nnlo(): + beta.beta_qcd((3, 0), nf) * a1**2 + beta.beta_qcd((4, 0), nf) * a1**3 ) - for method in ["iterate-exact"]: + for method in [EvoMethods.ITERATE_EXACT]: rhs = r * ns.dispatcher( (3, 0), method, gamma_ns, a1, a0, nf, ev_op_iterations ) @@ -172,7 +162,7 @@ def test_ode_n3lo(): + beta.beta_qcd((4, 0), nf) * a1**3 + beta.beta_qcd((5, 0), nf) * a1**4 ) - for method in ["iterate-exact"]: + for method in [EvoMethods.ITERATE_EXACT]: rhs = r * ns.dispatcher( (4, 0), method, gamma_ns, a1, a0, nf, ev_op_iterations ) @@ -202,7 +192,9 @@ def test_ode_n3lo(): def test_error(monkeypatch): monkeypatch.setattr("eko.beta.beta_qcd", lambda *_args: 1.0) with pytest.raises(NotImplementedError, match="order is not implemented"): - ns.dispatcher((5, 0), "iterate-exact", np.random.rand(3) + 0j, 0.2, 0.1, 3, 10) + ns.dispatcher( + (5, 0), EvoMethods.ITERATE_EXACT, np.random.rand(3) + 0j, 0.2, 0.1, 3, 10 + ) with pytest.raises(NotImplementedError): ad.gamma_ns((2, 0), 10202, 1, (0, 0, 0, 0, 0, 0, 0), 3) @@ -216,7 +208,7 @@ def test_gamma_usage(): gamma_ns = np.full(4, np.nan) for order in range(1, 5): gamma_ns[order - 1] = np.random.rand() - for method in methods: + for method in EvoMethods: r = ns.dispatcher( (order, 0), method, gamma_ns, a1, a0, nf, ev_op_iterations ) @@ -225,8 +217,8 @@ def test_gamma_usage(): for order in range(1, 5): gamma_ns = np.random.rand(order) gamma_ns[order - 1] = np.nan - for method in methods: - if method == "ordered-truncated": + for method in EvoMethods: + if method is EvoMethods.ORDERED_TRUNCATED: # we are actually dividing by a np.nan, # since the sum of U vec is nan warnings.simplefilter("ignore", RuntimeWarning) diff --git a/tests/eko/kernels/test_s.py b/tests/eko/kernels/test_s.py index 5380ba319..ff4fdfc0d 100644 --- a/tests/eko/kernels/test_s.py +++ b/tests/eko/kernels/test_s.py @@ -3,20 +3,10 @@ import numpy as np import pytest +from eko.kernels import EvoMethods from eko.kernels import singlet as s from ekore import anomalous_dimensions as ad -methods = [ - "iterate-expanded", - "decompose-expanded", - "perturbative-expanded", - "truncated", - "ordered-truncated", - "iterate-exact", - "decompose-exact", - "perturbative-exact", -] - def test_zero_lo(monkeypatch): """No evolution results in exp(0)""" @@ -35,7 +25,7 @@ def test_zero_lo(monkeypatch): np.array([[0, 0], [0, 1]]), ), ) - for method in methods: + for method in EvoMethods: np.testing.assert_allclose( s.dispatcher( (1, 0), method, gamma_s, 1, 1, nf, ev_op_iterations, ev_op_max_order @@ -75,8 +65,8 @@ def test_zero_nlo_decompose(monkeypatch): ), ) for method in [ - "decompose-expanded", - "decompose-exact", + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.DECOMPOSE_EXACT, ]: np.testing.assert_allclose( s.dispatcher( @@ -117,8 +107,8 @@ def test_zero_nnlo_decompose(monkeypatch): ), ) for method in [ - "decompose-expanded", - "decompose-exact", + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.DECOMPOSE_EXACT, ]: np.testing.assert_allclose( s.dispatcher( @@ -159,8 +149,8 @@ def test_zero_n3lo_decompose(monkeypatch): ), ) for method in [ - "decompose-expanded", - "decompose-exact", + EvoMethods.DECOMPOSE_EXPANDED, + EvoMethods.DECOMPOSE_EXACT, ]: np.testing.assert_allclose( s.dispatcher( @@ -200,7 +190,7 @@ def test_similarity(): ]: ref = s.dispatcher( order, - "decompose-exact", + EvoMethods.DECOMPOSE_EXACT, gamma_s, a1, a0, @@ -208,7 +198,7 @@ def test_similarity(): ev_op_iterations, ev_op_max_order, ) - for method in methods: + for method in EvoMethods: np.testing.assert_allclose( s.dispatcher( order, @@ -227,7 +217,9 @@ def test_similarity(): def test_error(): with pytest.raises(NotImplementedError): - s.dispatcher((4, 0), "AAA", np.random.rand(3, 2, 2), 0.2, 0.1, 3, 10, 10) + s.dispatcher( + (4, 0), "iterate-exact", np.random.rand(3, 2, 2), 0.2, 0.1, 3, 10, 10 + ) def mk_almost_diag_matrix(n, max_ang=np.pi / 8.0): @@ -247,7 +239,7 @@ def test_gamma_usage(): gamma_s = np.full((4, 2, 2), np.nan) for order in range(1, 5): gamma_s[order - 1] = mk_almost_diag_matrix(1) - for method in methods: + for method in EvoMethods: r = s.dispatcher( (order, 0), method, @@ -263,8 +255,8 @@ def test_gamma_usage(): for order in range(1, 5): gamma_s = mk_almost_diag_matrix(4) gamma_s[order - 1] = np.full((2, 2), np.nan) - for method in methods: - if "iterate" in method: + for method in EvoMethods: + if method in [EvoMethods.ITERATE_EXACT, EvoMethods.ITERATE_EXPANDED]: # we are actually dividing by the determinant of # matrix full of np.nan warnings.simplefilter("ignore", RuntimeWarning) @@ -287,8 +279,8 @@ def test_singlet_back(): nf = 4 a1 = 3.0 a0 = 4.0 - s10 = s.dispatcher(order, "iterate-exact", gamma_s, a1, a0, nf, 15, 1) + s10 = s.dispatcher(order, EvoMethods.ITERATE_EXACT, gamma_s, a1, a0, nf, 15, 1) np.testing.assert_allclose( np.linalg.inv(s10), - s.dispatcher(order, "iterate-exact", gamma_s, a0, a1, nf, 15, 1), + s.dispatcher(order, EvoMethods.ITERATE_EXACT, gamma_s, a0, a1, nf, 15, 1), ) diff --git a/tests/eko/quantities/test_couplings.py b/tests/eko/quantities/test_couplings.py index 0a3b97415..00ccae9d1 100644 --- a/tests/eko/quantities/test_couplings.py +++ b/tests/eko/quantities/test_couplings.py @@ -3,7 +3,7 @@ def test_couplings_ref(): scale = 90.0 - d = dict(alphas=0.1, alphaem=0.01, scale=scale, max_num_flavs=6, num_flavs_ref=None) + d = dict(alphas=0.1, alphaem=0.01, scale=scale, max_num_flavs=6, num_flavs_ref=5) couplings = CouplingsInfo.from_dict(d) assert couplings.scale == scale assert not couplings.em_running diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index c0dd67236..8ade8bcac 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -11,7 +11,7 @@ from eko import basis_rotation as br from eko.beta import beta_qcd_as2, beta_qcd_as3 from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo -from eko.kernels import non_singlet, singlet +from eko.kernels import EvoMethods, non_singlet, singlet from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import expanded, exponentiated from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet @@ -19,7 +19,7 @@ NF = 4 Q02 = 1.65**2 Q12 = 100**2 -EV_METHOD = "truncated" +EV_METHOD = EvoMethods.TRUNCATED def compute_a_s(q2, order): diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 6af22dd51..bfacc1c96 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,7 @@ from eko import basis_rotation as br from eko.couplings import CouplingEvolutionMethod, Couplings, CouplingsInfo -from eko.kernels import non_singlet, singlet +from eko.kernels import EvoMethods, non_singlet, singlet from eko.quantities.heavy_quarks import QuarkMassScheme from eko.scale_variations import Modes, expanded from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet @@ -10,7 +10,7 @@ NF = 4 Q02 = 1.65**2 Q12 = 100**2 -EV_METHOD = "truncated" +EV_METHOD = EvoMethods.TRUNCATED def compute_a_s(q2, order): diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py index ef32fb0d0..f8e051e4c 100644 --- a/tests/eko/test_couplings.py +++ b/tests/eko/test_couplings.py @@ -55,7 +55,7 @@ def test_init(self): alphas=alpharef[0], alphaem=alpharef[1], scale=muref, - num_flavs_ref=None, + num_flavs_ref=5, max_num_flavs=6, ) ) @@ -153,18 +153,19 @@ def test_ref(self): (0, np.inf, np.inf), (2, 4, 175), ] + nfrefs = (3, 4, 5) alpharef = (0.118, 0.00781) muref = 91.0 - couplings = CouplingsInfo.from_dict( - dict( - alphas=alpharef[0], - alphaem=alpharef[1], - scale=muref, - num_flavs_ref=None, - max_num_flavs=6, + for thresh_setup, nfref in zip(thresh_setups, nfrefs): + couplings = CouplingsInfo.from_dict( + dict( + alphas=alpharef[0], + alphaem=alpharef[1], + scale=muref, + num_flavs_ref=nfref, + max_num_flavs=6, + ) ) - ) - for thresh_setup in thresh_setups: for order_s in [1, 2, 3, 4]: for order_em in [0, 1, 2]: for evmod in CouplingEvolutionMethod: @@ -221,9 +222,10 @@ def test_exact(self): (0, np.inf, np.inf), (2, 4, 175), ] + nfrefs = (3, 4, 5) alpharef = np.array([0.118, 0.00781]) muref = 91.0 - for thresh_setup in thresh_setups: + for thresh_setup, nfref in zip(thresh_setups, nfrefs): for qcd in range(1, 4 + 1): for qed in range(2 + 1): for em_running in [ @@ -236,7 +238,7 @@ def test_exact(self): alphas=alpharef[0], alphaem=alpharef[1], scale=muref, - num_flavs_ref=None, + num_flavs_ref=nfref, max_num_flavs=6, em_running=em_running, ) @@ -299,7 +301,7 @@ def benchmark_expanded_n3lo(self): alphas=alpharef[0], alphaem=alpharef[1], scale=muref, - num_flavs_ref=None, + num_flavs_ref=5, max_num_flavs=6, ) m2c = 2 diff --git a/tests/ekobox/test_evol_pdf.py b/tests/ekobox/test_evol_pdf.py index 54bce6623..141787b49 100644 --- a/tests/ekobox/test_evol_pdf.py +++ b/tests/ekobox/test_evol_pdf.py @@ -1,4 +1,5 @@ import numpy as np +import pytest from banana import toy import eko @@ -50,6 +51,21 @@ def test_evolve_pdfs_dump_path(fake_lhapdf, cd): assert p.exists() +def test_evolve_pdfs_bad_scales(fake_lhapdf, cd): + """Bad scale configurations.""" + theory, op = init_cards() + n = "test_evolve_pdfs_bad_scales" + op = cards.example.operator() + op.mugrid = [(5.0, 3), (15.0, 4), (10.0, 5), (100.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + with cd(fake_lhapdf): + ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf) + op.mugrid = [(5.0, 3), (10.0, 3), (10.0, 4), (15.0, 4), (10.0, 5)] + with pytest.raises(ValueError, match="is bigger"): + with cd(fake_lhapdf): + ev_p.evolve_pdfs([toy.mkPDF("", 0)], theory, op, name=n, path=fake_lhapdf) + + def test_evolve_pdfs_dump_file(fake_lhapdf, cd): theory, op = init_cards() n = "test_evolve_pdfs_dump_file" diff --git a/tests/ekobox/test_info_file.py b/tests/ekobox/test_info_file.py index 64af8eb93..acf3f7272 100644 --- a/tests/ekobox/test_info_file.py +++ b/tests/ekobox/test_info_file.py @@ -23,3 +23,27 @@ def test_build(): np.testing.assert_allclose(info["QMin"], math.sqrt(op.mu2grid[0]), rtol=1e-5) assert info["XMin"] == op.xgrid.raw[0] assert info["XMax"] == op.xgrid.raw[-1] == 1.0 + + +def test_build_alphas_good(): + """Good configurations.""" + theory = cards.example.theory() + theory.order = (2, 0) + theory.couplings.alphas = 0.2 + op = cards.example.operator() + op.mu0 = 1.0 + # base case + op.mugrid = [(100.0, 5), (10.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 2 + np.testing.assert_allclose(info["AlphaS_Qs"], [10.0, 100.0]) + # several nf + op.mugrid = [(5.0, 4), (10.0, 5), (1.0, 3), (5.0, 3), (10.0, 5), (100.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 6 + np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 5.0, 5.0, 10.0, 10.0, 100.0]) + # several nf with gap + op.mugrid = [(1.0, 3), (10.0, 3), (10.0, 5), (100.0, 5)] + info = info_file.build_alphas(theory, op) + assert len(info["AlphaS_Vals"]) == 4 + np.testing.assert_allclose(info["AlphaS_Qs"], [1.0, 10.0, 10.0, 100.0])