Merge branch 'master' into nix

.pre-commit-config.yaml

@@ -13,8 +13,8 @@ repos:
       - id: check-merge-conflict
       - id: debug-statements
         exclude: "debug.py"
-  - repo: https://github.com/psf/black
-    rev: 23.10.1
+  - repo: https://github.com/psf/black-pre-commit-mirror
+    rev: 23.11.0
     hooks:
       - id: black
   - repo: https://github.com/asottile/blacken-docs
@@ -31,7 +31,7 @@ repos:
     hooks:
       - id: pyupgrade
   - repo: https://github.com/hadialqattan/pycln
-    rev: v2.3.0
+    rev: v2.4.0
     hooks:
       - id: pycln
         args: [--config=pyproject.toml]

extras/heavy_cc_nnlo/.gitignore

@@ -0,0 +1,3 @@
+data/
+demoFelix/
+plots/

extras/heavy_cc_nnlo/reformat/__init__.py

@@ -0,0 +1,40 @@
+# -*- coding: utf-8 -*-
+import argparse
+import logging
+import pathlib
+
+from . import generate as gen
+from . import parse
+
+_logger = logging.getLogger(__name__)
+
+
+def parse_args():
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument("paths", nargs="*", type=pathlib.Path)
+    parser.add_argument("-a", "--actions", nargs="*", choices=["dump", "plot"])
+
+    return parser.parse_args()
+
+
+def main():
+    args = parse_args()
+
+    for path in args.paths:
+        ch, blocks = parse.parse(path)
+        ds = gen.to_xarray(blocks)
+
+        if "dump" in args.actions:
+            dumpdest = pathlib.Path.cwd() / "data" / f"{path.stem}.nc"
+            dumpdest.parent.mkdir(exist_ok=True)
+            ds.to_netcdf(dumpdest)
+
+            _logger.info(
+                f"Saved xgrid array to {dumpdest.relative_to(pathlib.Path.cwd())}"
+            )
+
+        if "plot" in args.actions:
+            dest = pathlib.Path.cwd() / "plots" / f"{path.stem}.png"
+            dest.parent.mkdir(exist_ok=True)
+            gen.plot(ds, dest)

extras/heavy_cc_nnlo/reformat/__main__.py

@@ -0,0 +1,12 @@
+# -*- coding: utf-8 -*-
+import logging
+
+from rich.logging import RichHandler
+
+from . import main
+
+logging.basicConfig(
+    level=logging.INFO, format="%(message)s", datefmt="[%X]", handlers=[RichHandler()]
+)
+
+main()

extras/heavy_cc_nnlo/reformat/generate.py

@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+import logging
+import os
+import pathlib
+
+import matplotlib.colors as clr
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+
+from . import parse
+
+_logger = logging.getLogger(__name__)
+
+
+def to_xarray(blocks: list[parse.Block]) -> xr.Dataset:
+    catdict = {
+        k: (["bl"], [getattr(dic, k) for dic in blocks]) for k in blocks[0].__dict__
+    }
+    catdict["grid"][0].extend(["x", "y"])
+    catdict["xgrid"][0].append("z")
+
+    return xr.Dataset(catdict)
+
+
+def plot(ds: xr.Dataset, path: os.PathLike):
+    path = pathlib.Path(path)
+
+    plt.figure()
+    ds["xgrid"].plot(norm=clr.LogNorm())
+    figpath = path.with_stem(f"{path.stem}-xgrid")
+    plt.savefig(figpath)
+
+    _logger.info(f"Saved xgrid plot to {figpath.relative_to(pathlib.Path.cwd())}")
+
+    plt.figure()
+    xbj = ds["xbj"]
+    nxbj = np.abs(np.log(xbj))
+    (nxbj / nxbj.max()).plot(label="xBj")
+    qc = ds["qcd_scale"]
+    (qc / qc.max()).plot(label="qcd")
+    en = ds["energy"]
+    nen = np.log(en)
+    (nen / nen.max()).plot(label="E")
+    plt.ylabel("logs and ratios")
+    plt.legend()
+    plt.title(
+        f"xBj = {xbj.min():0.3g} - {xbj.max():0.3g}\n"
+        f"qcd = {qc.min()} - {qc.max()}\n"
+        f"E = {en.min():0.3g} - {en.max():0.3g}"
+    )
+    plt.tight_layout()
+    figpath = path.with_stem(f"{path.stem}-xbj-scales")
+    plt.savefig(figpath)
+
+    _logger.info(f"Saved xgrid plot to {figpath.relative_to(pathlib.Path.cwd())}")

extras/heavy_cc_nnlo/reformat/parse.py

@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+import os
+import pathlib
+import re
+from dataclasses import dataclass
+
+import numpy as np
+
+DOUBLE = r"-?\d\.\d+D[-+]\d+"
+INT = r"-?\d+"
+
+BLOCK = r"^<Block([ \d]+)>(.*?)^<End([ \d]+)>"
+KINS = re.compile(
+    rf"^ Bjorken.*\n ({DOUBLE}) +({DOUBLE}) +({DOUBLE}) +({INT})", flags=re.MULTILINE
+)
+SCALES = re.compile(rf"^ Charm.*\n ({DOUBLE}) +({DOUBLE})", flags=re.MULTILINE)
+XGRID = re.compile(r"^ x-grid +(\d+)", flags=re.MULTILINE)
+INTERPOLATION = re.compile(
+    rf"^ Interpolation coefficients: +({INT}) +({INT}) +({INT}) +({INT})",
+    flags=re.MULTILINE,
+)
+
+
+@dataclass
+class Block:
+    xbj: float
+    ybj: float
+    energy: float
+    cpc: int
+    mc: float
+    qcd_scale: float
+    xgrid: np.ndarray
+    grid: np.ndarray
+
+
+def parse(path: os.PathLike) -> tuple[int, list[Block]]:
+    path = pathlib.Path(path)
+    content = path.read_text(encoding="utf-8")
+
+    charge = 1 if path.stem.startswith("nu") else -1
+
+    blocks = []
+    for bl in re.finditer(BLOCK, content, flags=re.MULTILINE | re.DOTALL):
+        blocks.append(block(bl))
+
+    return charge, blocks
+
+
+def dtof(double: str):
+    return float(double.replace("D", "e"))
+
+
+def block(matched: re.Match) -> Block:
+    if int(matched[1]) != int(matched[3]):
+        raise ValueError
+
+    content = matched[2]
+
+    parsed = {}
+
+    kins = re.search(KINS, content)
+    if kins is None:
+        raise ValueError()
+    parsed["xbj"] = dtof(kins[1])
+    parsed["ybj"] = dtof(kins[2])
+    parsed["energy"] = dtof(kins[3])
+    parsed["cpc"] = int(kins[4])
+
+    scales = re.search(SCALES, content)
+    if scales is None:
+        raise ValueError()
+    parsed["mc"] = dtof(scales[1])
+    parsed["qcd_scale"] = dtof(scales[2])
+
+    nx = int(re.search(XGRID, content)[1])
+    xgrid = [
+        dtof(x)
+        for x in content.split("x-grid")[1].split("Interpolation")[0].splitlines()[1:-1]
+    ]
+    parsed["xgrid"] = xgrid
+
+    if len(xgrid) != nx:
+        raise ValueError(
+            f"xgrid actual length, {len(xgrid)}, different from what declared, {nx}"
+        )
+
+    coeffs = re.search(INTERPOLATION, content)
+    if coeffs is None:
+        raise ValueError()
+    a, nx1, c, _ = int(coeffs[1]), int(coeffs[2]), int(coeffs[3]), int(coeffs[4])
+    grid = np.array(
+        [
+            [dtof(val) for val in line.strip().split()]
+            for line in content.split("Interpolation")[1].splitlines()[1:]
+        ]
+    )
+    parsed["grid"] = grid
+
+    if grid.shape != (a * nx1, c):
+        raise ValueError()
+
+    return Block(**parsed)

extras/heavy_em_n3lo/.gitignore

@@ -0,0 +1 @@
+*.npy

extras/heavy_em_n3lo/produce_grids.py

@@ -0,0 +1,18 @@
+import pathlib
+
+import numpy as np
+
+here = pathlib.Path(__file__).parent / "yadism_grids"
+here.mkdir(exist_ok=True)
+
+
+def produce_grids():
+    etalist = np.geomspace(1e-5, 9.9e5, 73, endpoint=True)
+    xilist = np.geomspace(1e-3, 1e4, 49, endpoint=True)
+
+    np.save(here / "eta.npy", etalist)
+    np.save(here / "xi.npy", xilist)
+
+
+if __name__ == "__main__":
+    produce_grids()

extras/heavy_em_n3lo/yad_grids.py

@@ -0,0 +1,111 @@
+import pathlib
+import sys
+import time
+from multiprocessing import Pool
+
+import adani
+import numpy as np
+
+here = pathlib.Path(__file__).parent / "yadism_grids"
+here.mkdir(exist_ok=True)
+
+nf = int(sys.argv[1])
+n_threads = int(sys.argv[2])
+channel = sys.argv[3]
+order = int(sys.argv[4])
+try:
+    variation = int(sys.argv[5])
+except IndexError:
+    variation = 0
+mufrac = 1.0
+verbose = True
+
+
+def x_eta(eta, m2Q2):
+    return 1.0 / (1.0 + 4.0 * m2Q2 * (eta + 1))
+
+
+def function_to_exe_in_parallel(pair):
+    eta, xi = pair
+    m2Q2 = 1 / xi
+    m2mu2 = 1 / xi
+    x = x_eta(eta, m2Q2)
+    if order == 3:
+        if channel == "2g":
+            return adani.C2_g3_approximation(
+                x, m2Q2, m2mu2, nf, v=variation, method_flag=1
+            )
+        elif channel == "2q":
+            return adani.C2_ps3_approximation(x, m2Q2, m2mu2, nf, v=variation)
+        elif channel == "Lg":
+            return adani.CL_g3_approximation(
+                x, m2Q2, m2mu2, nf, v=variation, method_flag=1
+            )
+        elif channel == "Lq":
+            return adani.CL_ps3_approximation(x, m2Q2, m2mu2, nf, v=variation)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+    # TODO: do we really want this?
+    ## NNLO approximated
+    elif order == 2:
+        if channel == "2g":
+            return adani.C2_g2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "2q":
+            return adani.C2_ps2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "Lg":
+            return adani.CL_g2_approximation(x, m2Q2, m2mu2, v=variation)
+        elif channel == "Lq":
+            return adani.CL_ps2_approximation(x, m2Q2, m2mu2, v=variation)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+    ## NNLO exact
+    elif order == 0:
+        if channel == "2g":
+            return adani.C2_g2(x, m2Q2, m2mu2)
+        elif channel == "2q":
+            return adani.C2_ps2(x, m2Q2, m2mu2)
+        elif channel == "Lg":
+            return adani.CL_g2(x, m2Q2, m2mu2)
+        elif channel == "Lq":
+            return adani.CL_ps2(x, m2Q2, m2mu2)
+        else:
+            raise ValueError("Set channel to one of these: 2g 2q Lg Lq")
+
+
+def run(n_threads, eta_grid, xi_grid):
+    grid = []
+    for xi in xi_grid:
+        for eta in eta_grid:
+            grid.append((eta, xi))
+    args = (function_to_exe_in_parallel, grid)
+    with Pool(n_threads) as pool:
+        result = pool.map(*args)
+    return result
+
+
+if __name__ == "__main__":
+    output_file = f"C{channel}_nf{nf}_var{variation}.npy"
+    etafname = here / "eta.npy"
+    eta_grid = np.load(etafname)
+    xifname = here / "xi.npy"
+    xi_grid = np.load(xifname)
+
+    if verbose:
+        print(
+            f"Computation of the grid for the coefficient function C{channel} for nf = {nf}, and µ/Q = {mufrac}, variation = {variation}"
+        )
+        print(f"Size of the grid (eta,xi) = ({len(eta_grid)},{len(xi_grid)})")
+        print(
+            "This may take a while (depending on the number of threads you choose). In order to spend this time, I would suggest you this interesting view:"
+        )
+        print("https://www.youtube.com/watch?v=53pG68KCUMI")
+
+    start = time.perf_counter()
+    res_vec = np.array(run(n_threads, eta_grid, xi_grid))
+    if verbose:
+        print("total running time: ", time.perf_counter() - start)
+
+    res_mat = res_vec.reshape(len(xi_grid), len(eta_grid))
+    if verbose:
+        print("Saving grid in ", output_file)
+    np.save(here / output_file, res_mat)