diff --git a/benchmarks/eko/benchmark_inverse_matching.py b/benchmarks/eko/benchmark_inverse_matching.py index 7ca3b1a6a..d9c16b0b2 100644 --- a/benchmarks/eko/benchmark_inverse_matching.py +++ b/benchmarks/eko/benchmark_inverse_matching.py @@ -35,6 +35,7 @@ xif=1.0, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), matching_order=[2, 0], + use_fhmruvv=False, ) # operator settings diff --git a/doc/source/shared/abbreviations.rst b/doc/source/shared/abbreviations.rst index 7c415269b..8a03601f5 100644 --- a/doc/source/shared/abbreviations.rst +++ b/doc/source/shared/abbreviations.rst @@ -77,6 +77,7 @@ .. |BFKL| replace:: :abbr:`BFKL (Balitsky-Fadin-Kuraev-Lipatov)` + .. external .. |yadism| replace:: @@ -91,7 +92,6 @@ .. |pineko| replace:: :pineko:`\ ` - .. |APFEL| raw:: html APFEL @@ -108,6 +108,9 @@ QCDNUM +.. |FHMRUVV| replace:: + :abbr:`FHMRUVV (Falcioni, Herzog, Moch, Ruijl, Ueda, Vermaseren and Vogt)` + .. |T| raw:: html ✓ diff --git a/extras/lh_bench_23/.gitignore b/extras/lh_bench_23/.gitignore index 11fa526db..2942ad98d 100644 --- a/extras/lh_bench_23/.gitignore +++ b/extras/lh_bench_23/.gitignore @@ -1,2 +1,3 @@ *.tar *.csv +*.dat diff --git a/extras/lh_bench_23/cfg.py b/extras/lh_bench_23/cfg.py index 51762de73..ef7c747ef 100644 --- a/extras/lh_bench_23/cfg.py +++ b/extras/lh_bench_23/cfg.py @@ -1,4 +1,5 @@ import copy +import pathlib from math import inf, nan import numpy as np @@ -8,8 +9,16 @@ from eko.io import runcards from eko.io.types import ReferenceRunning +here = pathlib.Path(__file__).parent +eko_dir = here / "ekos" +table_dir = here / "tables" + + _sqrt2 = float(np.sqrt(2)) +# setup x rotation +xgrid = np.array([1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9]) + # theory settings # --------------- _t_vfns = dict( @@ -30,8 +39,9 @@ matching_ratios=[1.0, 1.0, 1.0], ), xif=1.0, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), matching_order=[2, 0], + use_fhmruvv=False, ) @@ -59,6 +69,19 @@ def ffns_theory(xif=1.0): return runcards.TheoryCard.from_dict(tt) +def n3lo_theory(ad_variation, is_ffns, use_fhmruvv=False, xif=1.0): + """Generate an N3LO theory card.""" + base = _t_ffns if is_ffns else _t_vfns + tt = copy.deepcopy(base) + tt["xif"] = xif + tt["order"] = [4, 0] + # here we keep the NNLO matching + tt["matching_order"] = [2, 0] + tt["n3lo_ad_variation"] = ad_variation + tt["use_fhmruvv"] = use_fhmruvv + return runcards.TheoryCard.from_dict(tt) + + # operator settings # ----------------- _o_vfns = dict( diff --git a/extras/lh_bench_23/plot_bench.py b/extras/lh_bench_23/plot_bench.py new file mode 100644 index 000000000..6887e7ac3 --- /dev/null +++ b/extras/lh_bench_23/plot_bench.py @@ -0,0 +1,49 @@ +from cfg import here, table_dir, xgrid +from utils import ( + compute_n3lo_avg_err, + compute_n3lo_nnlo_diff, + load_n3lo_tables, + load_nnlo_table, + plot_diff_to_nnlo, + plot_pdfs, +) + +USE_LINX = True +REL_DIFF = True +SCHEME = "VFNS" +SV = "central" + +plot_dir = here / "plots" +n3lo_table_dir = table_dir # / SCHEME + + +# load tables +eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="EKO") +fhmv_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="FHMV") +nnlo_central = load_nnlo_table(table_dir, SCHEME, SV) + +# compute avg and std +eko_res = compute_n3lo_avg_err(eko_dfs) +fhmv_res = compute_n3lo_avg_err(fhmv_dfs) +# eko_4mom_res = = compute_n3lo_avg_err(eko_dfs_4mom) + +n3lo_dfs = [ + (eko_res, "aN3LO EKO"), + (fhmv_res, "aN3LO FHMV"), + # (eko_4mom_res, "aN3LO EKO 4 mom"), +] + +# PDFs plots +plot_pdfs(xgrid, n3lo_dfs, nnlo_central, SCHEME, USE_LINX, plot_dir) + +# relative diff plots +eko_diff = compute_n3lo_nnlo_diff(eko_res, nnlo_central, REL_DIFF) +fhmv_diff = compute_n3lo_nnlo_diff(fhmv_res, nnlo_central, REL_DIFF) +n3lo_dfs = [ + (eko_diff, "aN3LO EKO"), + (fhmv_diff, "aN3LO FHMV"), + # (eko_4mom_res, "aN3LO EKO 4 mom"), +] + +# relative, absolute diff plots +plot_diff_to_nnlo(xgrid, n3lo_dfs, SCHEME, USE_LINX, plot_dir, REL_DIFF) diff --git a/extras/lh_bench_23/plot_bench_evol.py b/extras/lh_bench_23/plot_bench_evol.py new file mode 100644 index 000000000..92972af14 --- /dev/null +++ b/extras/lh_bench_23/plot_bench_evol.py @@ -0,0 +1,45 @@ +from cfg import here, table_dir, xgrid +from utils import ( + compute_n3lo_avg_err, + compute_n3lo_nnlo_diff, + load_n3lo_tables, + load_nnlo_table, + plot_diff_to_nnlo, + plot_pdfs, +) + +USE_LINX = False +REL_DIFF = True +SCHEME = "VFNS" +SV = "central" + +plot_dir = here / "plots_evol" +n3lo_table_dir = table_dir # / SCHEME + + +# load tables +eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="EKO", rotate_to_evol=True) +fhmv_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="FHMV", rotate_to_evol=True) +nnlo_central = load_nnlo_table(table_dir, SCHEME, SV, rotate_to_evol=True) + +# compute avg and std +eko_res = compute_n3lo_avg_err(eko_dfs) +fhmv_res = compute_n3lo_avg_err(fhmv_dfs) + +n3lo_dfs = [ + (eko_res, "aN3LO EKO"), + (fhmv_res, "aN3LO FHMV"), +] + +# absolute plots +plot_pdfs(xgrid, n3lo_dfs, nnlo_central, SCHEME, USE_LINX, plot_dir) + +# relative, absolute diff plots +eko_diff = compute_n3lo_nnlo_diff(eko_res, nnlo_central, REL_DIFF) +fhmv_diff = compute_n3lo_nnlo_diff(fhmv_res, nnlo_central, REL_DIFF) +n3lo_dfs = [ + (eko_diff, "aN3LO EKO"), + (fhmv_diff, "aN3LO FHMV"), +] + +plot_diff_to_nnlo(xgrid, n3lo_dfs, SCHEME, USE_LINX, plot_dir, REL_DIFF) diff --git a/extras/lh_bench_23/plot_bench_msht.py b/extras/lh_bench_23/plot_bench_msht.py new file mode 100644 index 000000000..0314082ee --- /dev/null +++ b/extras/lh_bench_23/plot_bench_msht.py @@ -0,0 +1,62 @@ +from cfg import here, table_dir, xgrid +from utils import ( + compute_n3lo_avg_err, + compute_n3lo_nnlo_diff, + load_msht, + load_n3lo_tables, + load_nnlo_table, + plot_diff_to_nnlo, + plot_pdfs, +) + +USE_LINX = False +REL_DIFF = True +SCHEME = "VFNS" +SV = "central" + +plot_dir = here / "plots_msht" +n3lo_table_dir = table_dir # / SCHEME +msht_table_dir = table_dir + + +# load tables +eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="EKO") +fhmv_eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, approx="FHMV") +msht_dfs = load_msht(msht_table_dir, SCHEME, approx="MSHT") +fhmv_msht_dfs = load_msht(msht_table_dir, SCHEME, approx="FHMV") +nnlo_central = load_nnlo_table(table_dir, SCHEME, SV) + +# compute avg and std +eko_res = compute_n3lo_avg_err(eko_dfs) +fhmv_eko_res = compute_n3lo_avg_err(fhmv_eko_dfs) +msht_res = compute_n3lo_avg_err(msht_dfs) +fhmv_msht_res = compute_n3lo_avg_err(fhmv_msht_dfs) +# eko_4mom_res = = compute_n3lo_avg_err(eko_dfs_4mom) + +n3lo_dfs = [ + (eko_res, "EKO"), + (fhmv_eko_res, "FHMV EKO"), + (msht_res, "MSHT"), + (fhmv_msht_res, "FHMV MSHT") + # (eko_4mom_res, "aN3LO EKO 4 mom"), +] + +# PDFs plots +plot_pdfs(xgrid, n3lo_dfs, nnlo_central, SCHEME, USE_LINX, plot_dir) + +# relative diff plots +eko_diff = compute_n3lo_nnlo_diff(eko_res, nnlo_central, REL_DIFF) +fhmv_eko_diff = compute_n3lo_nnlo_diff(fhmv_eko_res, nnlo_central, REL_DIFF) +msht_diff = compute_n3lo_nnlo_diff(msht_res, nnlo_central, REL_DIFF) +fhmv_msht_diff = compute_n3lo_nnlo_diff(fhmv_msht_res, nnlo_central, REL_DIFF) + +n3lo_dfs = [ + (eko_diff, "EKO"), + (fhmv_eko_diff, "FHMV EKO"), + (msht_diff, "MSHT"), + (fhmv_msht_diff, "FHMV MSHT") + # (eko_4mom_res, "aN3LO EKO 4 mom"), +] + +# relative, absolute diff plots +plot_diff_to_nnlo(xgrid, n3lo_dfs, SCHEME, USE_LINX, plot_dir, REL_DIFF) diff --git a/extras/lh_bench_23/plotstyle.mplstyle b/extras/lh_bench_23/plotstyle.mplstyle new file mode 100644 index 000000000..d82ce87b4 --- /dev/null +++ b/extras/lh_bench_23/plotstyle.mplstyle @@ -0,0 +1,89 @@ +### FIGURE +figure.facecolor : white +figure.dpi : 150 +figure.titlesize : 16 +figure.titleweight : bold +figure.figsize : 6.4, 4.8 + + +### GRIDS +grid.color : 0.8 +grid.linestyle : -- +grid.alpha : 0.5 + + +### PATCHES +patch.facecolor : 4C72B0 +patch.linewidth : 0.7 + +#### IMAGE +image.cmap : Greys + +### LEGEND +legend.fontsize : 12 +legend.loc : best +legend.fancybox : True +legend.framealpha : 0.3 +legend.facecolor : white +legend.frameon : True +legend.handleheight : 1 +legend.numpoints : 1 +legend.scatterpoints : 1 + + +### XTICKS +xtick.color : ".6" +xtick.labelcolor : ".4" +xtick.labelsize : 10.0 +xtick.direction : in +xtick.major.size : 4.0 +xtick.major.width : 1.0 +xtick.major.top : true +xtick.major.bottom : true +xtick.minor.size : 2.0 +xtick.minor.width : 0.5 +xtick.minor.top : true +xtick.minor.bottom : true + +### XTICKS +ytick.color : ".6" +ytick.labelcolor : ".4" +ytick.labelsize : 10.0 +ytick.direction : in +ytick.major.size : 4.0 +ytick.major.width : 1.0 +ytick.major.left : true +ytick.major.right : true +ytick.minor.size : 2.0 +ytick.minor.width : 0.5 +ytick.minor.left : true +ytick.minor.right : true + +### FONTS +font.family : serif +font.sans-serif : Arial +font.size : 13 +font.weight : bold + +### TEXTS +text.usetex : True +text.latex.preamble : \usepackage{amsmath} +text.color : ".15" + +### AXES +axes.grid : True +axes.axisbelow : True +axes.facecolor : white +axes.edgecolor : 0.8 +axes.linewidth : 1.0 +axes.titlesize : 16 +axes.titleweight : bold +axes.labelsize : 16 +axes.labelcolor : 0.15 +axes.titlepad : 9.0 +axes.xmargin : 0 +axes.labelweight : bold +axes.spines.top : True +axes.spines.right : True +# axes.prop_cycle : cycler('color', ['001D66', 'da3b46', 'F6AE2D', '343434', '2ca02c', '9467bd', '8c564b']) +# axes.prop_cycle : cycler('color', ["001D66", "B01C64", "3D2E85", "FFBF33", "9467bd", "FF6000"]) diff --git a/extras/lh_bench_23/run-n3lo.py b/extras/lh_bench_23/run-n3lo.py new file mode 100644 index 000000000..4ee00d520 --- /dev/null +++ b/extras/lh_bench_23/run-n3lo.py @@ -0,0 +1,137 @@ +import argparse +import logging +import pathlib +import sys + +import pandas as pd +import yaml +from banana import toy +from cfg import ( + _sqrt2, + eko_dir, + ffns_labels, + ffns_operator, + ffns_rotate_to_LHA, + n3lo_theory, + table_dir, + vfns_labels, + vfns_operator, + vfns_rotate_to_LHA, + xgrid, +) + +import eko +from eko.runner.managed import solve +from ekobox import apply +from ekomark.benchmark.external.LHA_utils import here as there + +# reference values +with open(there / "LHA.yaml", encoding="utf-8") as o: + ref_data = yaml.safe_load(o) + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("scheme", help="FFNS or VFNS?") + parser.add_argument("sv", help="scale variation: up, central, or down") + parser.add_argument( + "ad_variation", + help="n3lo anomalous dimension variation: (gg, gq, qg, qq, nsp, nsm, nsv)", + nargs="*", + type=int, + ) + parser.add_argument( + "--use_fhmruvv", + help="Use the FHMRUVV N3LO ad approximation", + action="store_true", + ) + parser.add_argument("--rerun", help="Rerun eko", action="store_true") + parser.add_argument( + "-v", "--verbose", help="Print eko log to screen", action="store_true" + ) + args = parser.parse_args() + + # determine xif + if "central".startswith(args.sv): + xif = 1.0 + sv = "central" + part = 1 + elif "up".startswith(args.sv): + xif = _sqrt2 + sv = "up" + part = 2 + elif "down".startswith(args.sv): + xif = 1.0 / _sqrt2 + sv = "down" + part = 3 + else: + raise ValueError( + "sv has to be up, central, or down - or any abbreviation there of" + ) + # determine scheme + if args.scheme == "FFNS": + scheme = "FFNS" + t = n3lo_theory( + ad_variation=args.ad_variation, + is_ffns=True, + use_fhmruvv=args.use_fhmruvv, + xif=xif, + ) + o = ffns_operator + tab = 14 + lab = ffns_labels + rot = ffns_rotate_to_LHA + elif args.scheme == "VFNS": + scheme = "VFNS" + t = n3lo_theory( + ad_variation=args.ad_variation, + is_ffns=False, + use_fhmruvv=args.use_fhmruvv, + xif=xif, + ) + o = vfns_operator + tab = 15 + lab = vfns_labels + rot = vfns_rotate_to_LHA + else: + raise ValueError("scheme has to be FFNS or VFNS") + + # eko path + eko_dir.mkdir(exist_ok=True) + approx_name = "FHMRUVV" if args.use_fhmruvv else "EKO" + var_name = "-".join([str(a) for a in args.ad_variation]) + p = pathlib.Path(f"{eko_dir}/{scheme}-{sv}-{var_name}-{approx_name}.tar") + + # recompute? + if not p.exists() or args.rerun: + print("(Re)running eko ...") + p.unlink(True) + if args.verbose: + logStdout = logging.StreamHandler(sys.stdout) + logStdout.setLevel(logging.INFO) + logStdout.setFormatter(logging.Formatter("%(message)s")) + logging.getLogger("eko").handlers = [] + logging.getLogger("eko").addHandler(logStdout) + logging.getLogger("eko").setLevel(logging.INFO) + solve(t, o, p) + + # apply PDF + out = {} + with eko.EKO.read(p) as eko_: + pdf = apply.apply_pdf_flavor(eko_, toy.mkPDF("ToyLH", 0), xgrid, rot, lab) + for lab, f in list(pdf.values())[0]["pdfs"].items(): + out[lab] = xgrid * f + + # display result + pd.set_option("display.float_format", "{:.4e}".format) + me = pd.DataFrame(out) + print("EKO") + print(me) + # dump to file + table_dir.mkdir(exist_ok=True) + me.to_csv(f"{table_dir}/table{scheme}-{sv}-{var_name}-{approx_name}.csv") + + # load reference + ref = pd.DataFrame(ref_data[f"table{tab}"][f"part{part}"]) + print() + print("rel. distance to reference") + print((me - ref) / ref) diff --git a/extras/lh_bench_23/run-nnlo.py b/extras/lh_bench_23/run-nnlo.py index 1a006f7bc..85ab2f984 100644 --- a/extras/lh_bench_23/run-nnlo.py +++ b/extras/lh_bench_23/run-nnlo.py @@ -8,14 +8,18 @@ import yaml from banana import toy from cfg import ( + _sqrt2, + eko_dir, ffns_labels, ffns_operator, ffns_rotate_to_LHA, ffns_theory, + table_dir, vfns_labels, vfns_operator, vfns_rotate_to_LHA, vfns_theory, + xgrid, ) import eko @@ -23,13 +27,6 @@ from ekobox import apply from ekomark.benchmark.external.LHA_utils import here as there -_sqrt2 = float(np.sqrt(2)) - - -# setup x rotation -xgrid = np.array([1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 0.3, 0.5, 0.7, 0.9]) - - # reference values with open(there / "LHA.yaml", encoding="utf-8") as o: ref_data = yaml.safe_load(o) @@ -81,7 +78,8 @@ raise ValueError("scheme has to be FFNS or VFNS") # eko path - p = pathlib.Path(f"{scheme}-{sv}.tar") + eko_dir.mkdir(exist_ok=True) + p = pathlib.Path(f"{eko_dir}/{scheme}-{sv}.tar") # recompute? if not p.exists() or args.rerun: @@ -109,7 +107,8 @@ print("EKO") print(me) # dump to file - me.to_csv(f"table{tab}-part{part}.csv") + table_dir.mkdir(exist_ok=True) + me.to_csv(f"{table_dir}/table{tab}-part{part}.csv") # load reference ref = pd.DataFrame(ref_data[f"table{tab}"][f"part{part}"]) diff --git a/extras/lh_bench_23/run_eko.sh b/extras/lh_bench_23/run_eko.sh new file mode 100755 index 000000000..579ccbb50 --- /dev/null +++ b/extras/lh_bench_23/run_eko.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +SCHEME="FFNS" +SV="central" +AD_VAR=(0 0 0 0 0 0 0) +AD_MAX_VAR_LIST=(19 21 15 6) + +# run the central +python run-n3lo.py $SCHEME $SV "${AD_VAR[@]}" + +# loop on gammas +for I in {0..3} + do + # loop on variations + VARIATIONS=${AD_MAX_VAR_LIST[$I]} + VAR_LIST=($(seq 1 1 $VARIATIONS)) + for VAR in "${VAR_LIST[@]}" + do + AD_VAR[$I]=$VAR + python run-n3lo.py $SCHEME $SV "${AD_VAR[@]}" + AD_VAR[$I]=0 + done + done diff --git a/extras/lh_bench_23/run_fhmv.sh b/extras/lh_bench_23/run_fhmv.sh new file mode 100755 index 000000000..45807ea4d --- /dev/null +++ b/extras/lh_bench_23/run_fhmv.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +SCHEME="FFNS" +SV="central" +AD_VAR=(0 0 0 0 0 0 0) + + +# run the central +python run-n3lo.py $SCHEME $SV "${AD_VAR[@]}" "--use_fhmv" + +# loop on gammas +for I in {0..6} + do + # loop on variations + for VAR in {1..2} + do + AD_VAR[$I]=$VAR + python run-n3lo.py $SCHEME $SV "${AD_VAR[@]}" "--use_fhmv" + AD_VAR[$I]=0 + done + done diff --git a/extras/lh_bench_23/utils.py b/extras/lh_bench_23/utils.py new file mode 100644 index 000000000..dc83a2c7e --- /dev/null +++ b/extras/lh_bench_23/utils.py @@ -0,0 +1,317 @@ +"""Plotting utils.""" + +import pathlib + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +FMT_LIST = ["x", "o", "v", "*"] +LHA_LABELS_MAP = { + "L_m": r"L^- = \bar{d} - \bar{u}", + "L_p": r"L^+ = 2(\bar{u} + \bar{d})", + "b_p": "b^+", + "c_p": "c^+", + "s_p": "s^+", +} + +HERE = pathlib.Path(__file__).parent +plt.style.use(HERE / "plotstyle.mplstyle") + + +def lha_labels(scheme: str) -> list: + """PDFs labels in the LHA basis.""" + if scheme == "FFNS": + pdf_labels = [ + "u_v", + "d_v", + r"L^- = \bar{d} - \bar{u}", + r"L^+ = 2(\bar{u} + \bar{d})", + "s_v", + "s^+", + "c^+", + "g", + ] + elif scheme == "VFNS": + pdf_labels = [ + "u_v", + "d_v", + r"L^- = \bar{d} - \bar{u}", + r"L^+ = 2(\bar{u} + \bar{d})", + "s^+", + "c^+", + "b^+", + "g", + ] + return pdf_labels + + +def evol_labels(scheme: str) -> list: + """PDFs labels in the Evolution basis.""" + if scheme == "FFNS": + pdf_labels = [ + "V", + "V_3", + "V_8", + "T_3", + "T_8", + "T_{15}", + r"\Sigma", + "g", + ] + elif scheme == "VFNS": + pdf_labels = [ + "V", + "V_3", + "T_3", + "T_8", + "T_{15}", + "T_{24}", + r"\Sigma", + "g", + ] + return pdf_labels + + +def rotate_lha_to_evol(df: pd.DataFrame, scheme: str) -> pd.DataFrame: + """Rotation from LHA to Evolution basis.""" + if scheme == "FFNS": + matrix = [ + [1, 1, 0, 0, 0, 0, 0, 0], # V + [1, -1, 0, 0, 0, 0, 0, 0], # V3 + [1, 1, 0, 0, -2, 0, 0, 0], # V8 + [1, -1, -2, 0, 0, 0, 0, 0], # T3 + [1, 1, 0, 1, 0, -2, 0, 0], # T8 + [1, 1, 0, 1, 0, 1, -3, 0], # T15 + [1, 1, 0, 1, 0, 1, 1, 0], # S + [0, 0, 0, 0, 0, 0, 0, 1], # g + ] + elif scheme == "VFNS": + matrix = [ + [1, 1, 0, 0, 0, 0, 0, 0], # V, V8 + [1, -1, 0, 0, 0, 0, 0, 0], # V3 + [1, -1, -2, 0, 0, 0, 0, 0], # T3 + [1, 1, 0, 1, -2, 0, 0, 0], # T8 + [1, 1, 0, 1, 1, -3, 0, 0], # T15 + [1, 1, 0, 1, 1, 1, -4, 0], # T24 + [1, 1, 0, 1, 1, 1, 1, 0], # S + [0, 0, 0, 0, 0, 0, 0, 1], # g + ] + rotated = (matrix @ df.values.T).T + return pd.DataFrame(rotated, columns=evol_labels(scheme)) + + +def load_n3lo_tables( + n3lo_table_dir: pathlib.Path, scheme: str, approx: str, rotate_to_evol: bool = False +) -> list: + """Load the N3LO tables.""" + dfs = [] + for p in n3lo_table_dir.iterdir(): + if scheme not in p.stem: + continue + if approx in p.stem: + table = pd.read_csv(p, index_col=0) + table.rename(columns=LHA_LABELS_MAP, inplace=True) + if rotate_to_evol: + table = rotate_lha_to_evol(table, scheme) + dfs.append(table) + return dfs + + +def load_nnlo_table( + table_dir, scheme, sv="central", rotate_to_evol: bool = False +) -> pd.DataFrame: + """Load the NNLO tables.""" + if scheme == "FFNS": + tab = 14 + elif scheme == "VFNS": + tab = 15 + if sv == "central": + part = 1 + + table = pd.read_csv(f"{table_dir}/table{tab}-part{part}.csv", index_col=0) + table.rename(columns=LHA_LABELS_MAP, inplace=True) + if rotate_to_evol: + table = rotate_lha_to_evol(table, scheme) + return table + + +def load_msht( + table_dir: pathlib.Path, scheme: str, approx: str, rotate_to_evol: bool = False +) -> list: + """Load MSHT files.""" + + if scheme != "VFNS": + raise ValueError(f"{scheme} not provided by MSHT, comment it out") + APPROX_MAP = { + "FHMV": "Moch", + "MSHT": "Posterior", + } + fhmv_msht_table_dir = table_dir / f"{scheme}_{APPROX_MAP[approx]}_numbers" + + columns = lha_labels(scheme) + # columns.insert(0,'x') + # columns.insert(0,'Q') + dfs = [] + + for p in fhmv_msht_table_dir.iterdir(): + data = np.loadtxt(p) + data = pd.DataFrame(data[:, 2:], columns=columns) + if rotate_to_evol: + data = rotate_lha_to_evol(data, scheme) + dfs.append(data) + return dfs + + +def compute_n3lo_avg_err(dfs: list) -> tuple: + """Compute N3LO average and error.""" + df_central = np.mean(dfs, axis=0) + df_central = pd.DataFrame(df_central, columns=dfs[0].columns) + # TODO: improve errors. + df_std = np.std(dfs, axis=0) + df_std = pd.DataFrame(df_std, columns=dfs[0].columns) + return df_central, df_std + + +def compute_n3lo_nnlo_diff(n3lo: tuple, nnlo: pd.DataFrame, rel_diff: bool) -> tuple: + """Compute relative / absolute differece to NNLO.""" + norm = nnlo + n3lo_central, n3lo_std = n3lo + if not rel_diff: + norm = pd.DataFrame(np.full(nnlo.shape, 1), columns=nnlo.columns) + diff = (n3lo_central - nnlo) / norm + diff_std = np.abs(n3lo_std / norm) + return diff, diff_std + + +def plot_pdfs( + xgrid: np.array, + n3lo_dfs: tuple, + nnlo_df: pd.DataFrame, + scheme: str, + use_linx: bool, + plot_dir: pathlib.Path, +) -> None: + """Absolute PDFs plots.""" + + fig, axs = plt.subplots(2, 4, figsize=(15, 7)) + + xcut = 4 if use_linx else 0 + xgrid = xgrid[xcut:] + + xscale = "linx" if use_linx else "logx" + plot_name = f"lh_n3lo_bench_{scheme}_{xscale}" + plot_dir.mkdir(exist_ok=True) + + fig.suptitle(f"{scheme}" + " $Q: \\sqrt{2} \\to 100 \\ GeV$") + + # loop on PDFs + for i, ax in enumerate( + axs.reshape( + 8, + ) + ): + # loop on n3lo + for j, (tabs, approx_label) in enumerate(n3lo_dfs): + central, err = tabs + ax.errorbar( + xgrid, + central.values[xcut:, i], + yerr=err.values[xcut:, i], + fmt=FMT_LIST[j], + label=approx_label, + capsize=5, + ) + # ax.errorbar( + # xgrid, + # eko_4mom_central[:, i], + # yerr=eko_4mom_std[:, i], + # fmt="x", + # label="aN3LO EKO (4 moments)", + # capsize=5, + # ) + ax.errorbar(xgrid, nnlo_df.values[xcut:, i], fmt=".", label="NNLO") + ax.hlines( + 0, + xgrid.min() - xgrid.min() / 3, + 1, + linestyles="dotted", + color="black", + linewidth=0.5, + ) + if not use_linx: + ax.set_xscale("log") + ax.set_xlabel("$x$") + ax.set_ylabel(f"${nnlo_df.columns[i]}$") + ax.set_xlim(xgrid.min() - xgrid.min() / 3, 1) + + plt.legend() + plt.tight_layout() + plt.savefig(f"{plot_dir}/{plot_name}.pdf") + + +def plot_diff_to_nnlo( + xgrid: np.array, + n3lo_dfs: tuple, + scheme: str, + use_linx: bool, + plot_dir: pathlib.Path, + rel_dff: bool, +) -> None: + """Difference w.r.t NNLO PDFs plots.""" + + fig, axs = plt.subplots(2, 4, figsize=(15, 7)) + + xcut = 4 if use_linx else 0 + xgrid = xgrid[xcut:] + xscale = "linx" if use_linx else "logx" + + diff_type = "rel_diff" if rel_dff else "abs_diff" + plot_name = f"lh_n3lo_bench_{scheme}_{xscale}_{diff_type}" + + diff_type = "Relative" if rel_dff else "Absolute" + fig.suptitle( + f"{diff_type} difference to NNLO, {scheme}" + " $Q: \\sqrt{2} \\to 100 \\ GeV$" + ) + + for i, ax in enumerate( + axs.reshape( + 8, + ) + ): + # loop on n3lo + for j, (tabs, approx_label) in enumerate(n3lo_dfs): + central, err = tabs + ax.errorbar( + xgrid, + central.values[xcut:, i], + yerr=err.values[xcut:, i], + fmt=FMT_LIST[j], + label=approx_label, + capsize=5, + ) + # ax.errorbar( + # xgrid, + # eko_4mom_diff.values[:, i], + # yerr=eko_4mom_diff_std.values[:, i], + # fmt="x", + # label="aN3LO EKO (4 moments)", + # capsize=5, + # ) + ax.hlines( + 0, + xgrid.min() - xgrid.min() / 3, + 1, + linestyles="dotted", + color="black", + linewidth=0.5, + ) + if not use_linx: + ax.set_xscale("log") + ax.set_xlabel("$x$") + ax.set_ylabel(f"${central.columns[i]}$") + ax.set_xlim(xgrid.min() - xgrid.min() / 3, 1) + + plt.legend() + plt.tight_layout() + plt.savefig(f"{plot_dir}/{plot_name}.pdf") diff --git a/extras/n3lo_bench/plot_msht.py b/extras/n3lo_bench/plot_msht.py index c922c1c03..8ddc92d2a 100644 --- a/extras/n3lo_bench/plot_msht.py +++ b/extras/n3lo_bench/plot_msht.py @@ -74,7 +74,7 @@ def msht_splitting(entry, x, nf, posterior): ) qqns = np.array( [ - msht.pqqps3a(x, ans0), + msht.p3nsa(x, ans0, 1, nf), msht.p3nsa(x, ans1, 1, nf), msht.p3nsa(x, ans2, 1, nf), ] diff --git a/extras/n3lo_bench/splitting_function_utils.py b/extras/n3lo_bench/splitting_function_utils.py index 6917e6bd8..6b32c86cc 100644 --- a/extras/n3lo_bench/splitting_function_utils.py +++ b/extras/n3lo_bench/splitting_function_utils.py @@ -12,12 +12,12 @@ map_non_singlet_modes = {"+": 10101, "-": 10201, "v": 10200} -def compute_ad(nf, n_grid, ns_mode=None, n3lo_variation="best"): +def compute_ad(nf, n_grid, ns_mode=None, n3lo_variation=(0, 0, 0, 0, 0, 0, 0)): ns_mode = map_non_singlet_modes.get(ns_mode, None) gs_list = [] for n in n_grid: if ns_mode is not None: - gs_list.append(gamma_ns((4, 0), ns_mode, n, nf).real) + gs_list.append(gamma_ns((4, 0), ns_mode, n, nf, n3lo_variation).real) else: gs_list.append(gamma_singlet((4, 0), n, nf, n3lo_variation).real) return np.array(gs_list) @@ -53,7 +53,7 @@ def integrand(u, x, order, entry, nf, ns_mode, n3lo_variation, L): idx1, idx2 = map_singlet_entries[entry] gamma = gamma[order, idx1, idx2] else: - gamma = gamma_ns((order + 1, 0), ns_mode, path.n, nf) + gamma = gamma_ns((order + 1, 0), ns_mode, path.n, nf, n3lo_variation) if L != 0: gamma = sv.exponentiated.gamma_variation(gamma, (order + 1, 0), nf, L) gamma = gamma[order] diff --git a/poetry.lock b/poetry.lock index 98f2235ce..f1898610a 100644 --- a/poetry.lock +++ b/poetry.lock @@ -130,13 +130,13 @@ files = [ [[package]] name = "banana-hep" -version = "0.6.9" +version = "0.6.12" description = "Benchmark QCD physics" optional = true python-versions = ">=3.8.0,<3.12" files = [ - {file = "banana_hep-0.6.9-py3-none-any.whl", hash = "sha256:e857b777eb0f64b0bf5e0a2b7cc6afcb5680a60dc3bdd7ce9c34c4bf3b94828b"}, - {file = "banana_hep-0.6.9.tar.gz", hash = "sha256:0b9f5cd3e4e44d5c8117bab93c832784009feb6d87d5a874da0d47fc91ab0514"}, + {file = "banana_hep-0.6.12-py3-none-any.whl", hash = "sha256:3cb065926a774689f6dee606dde0453e4e856dc1372c99e12f6b54aca26c41dd"}, + {file = "banana_hep-0.6.12.tar.gz", hash = "sha256:0cf4a78f5affa5e1afe138ae6c9df5414deb9785901b585fcbdf6bbf531a1b9b"}, ] [package.dependencies] @@ -430,7 +430,6 @@ files = [ {file = "contourpy-1.1.0-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:18a64814ae7bce73925131381603fff0116e2df25230dfc80d6d690aa6e20b37"}, {file = "contourpy-1.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:90c81f22b4f572f8a2110b0b741bb64e5a6427e0a198b2cdc1fbaf85f352a3aa"}, {file = "contourpy-1.1.0-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:53cc3a40635abedbec7f1bde60f8c189c49e84ac180c665f2cd7c162cc454baa"}, - {file = "contourpy-1.1.0-cp310-cp310-win32.whl", hash = "sha256:9b2dd2ca3ac561aceef4c7c13ba654aaa404cf885b187427760d7f7d4c57cff8"}, {file = "contourpy-1.1.0-cp310-cp310-win_amd64.whl", hash = "sha256:1f795597073b09d631782e7245016a4323cf1cf0b4e06eef7ea6627e06a37ff2"}, {file = "contourpy-1.1.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:0b7b04ed0961647691cfe5d82115dd072af7ce8846d31a5fac6c142dcce8b882"}, {file = "contourpy-1.1.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:27bc79200c742f9746d7dd51a734ee326a292d77e7d94c8af6e08d1e6c15d545"}, @@ -439,7 +438,6 @@ files = [ {file = "contourpy-1.1.0-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:e5cec36c5090e75a9ac9dbd0ff4a8cf7cecd60f1b6dc23a374c7d980a1cd710e"}, {file = "contourpy-1.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:1f0cbd657e9bde94cd0e33aa7df94fb73c1ab7799378d3b3f902eb8eb2e04a3a"}, {file = "contourpy-1.1.0-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:181cbace49874f4358e2929aaf7ba84006acb76694102e88dd15af861996c16e"}, - {file = "contourpy-1.1.0-cp311-cp311-win32.whl", hash = "sha256:edb989d31065b1acef3828a3688f88b2abb799a7db891c9e282df5ec7e46221b"}, {file = "contourpy-1.1.0-cp311-cp311-win_amd64.whl", hash = "sha256:fb3b7d9e6243bfa1efb93ccfe64ec610d85cfe5aec2c25f97fbbd2e58b531256"}, {file = "contourpy-1.1.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:bcb41692aa09aeb19c7c213411854402f29f6613845ad2453d30bf421fe68fed"}, {file = "contourpy-1.1.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:5d123a5bc63cd34c27ff9c7ac1cd978909e9c71da12e05be0231c608048bb2ae"}, @@ -448,7 +446,6 @@ files = [ {file = "contourpy-1.1.0-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:317267d915490d1e84577924bd61ba71bf8681a30e0d6c545f577363157e5e94"}, {file = "contourpy-1.1.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d551f3a442655f3dcc1285723f9acd646ca5858834efeab4598d706206b09c9f"}, {file = "contourpy-1.1.0-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:e7a117ce7df5a938fe035cad481b0189049e8d92433b4b33aa7fc609344aafa1"}, - {file = "contourpy-1.1.0-cp38-cp38-win32.whl", hash = "sha256:108dfb5b3e731046a96c60bdc46a1a0ebee0760418951abecbe0fc07b5b93b27"}, {file = "contourpy-1.1.0-cp38-cp38-win_amd64.whl", hash = "sha256:d4f26b25b4f86087e7d75e63212756c38546e70f2a92d2be44f80114826e1cd4"}, {file = "contourpy-1.1.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:bc00bb4225d57bff7ebb634646c0ee2a1298402ec10a5fe7af79df9a51c1bfd9"}, {file = "contourpy-1.1.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:189ceb1525eb0655ab8487a9a9c41f42a73ba52d6789754788d1883fb06b2d8a"}, @@ -457,7 +454,6 @@ files = [ {file = "contourpy-1.1.0-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:143dde50520a9f90e4a2703f367cf8ec96a73042b72e68fcd184e1279962eb6f"}, {file = "contourpy-1.1.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e94bef2580e25b5fdb183bf98a2faa2adc5b638736b2c0a4da98691da641316a"}, {file = "contourpy-1.1.0-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:ed614aea8462735e7d70141374bd7650afd1c3f3cb0c2dbbcbe44e14331bf002"}, - {file = "contourpy-1.1.0-cp39-cp39-win32.whl", hash = "sha256:71551f9520f008b2950bef5f16b0e3587506ef4f23c734b71ffb7b89f8721999"}, {file = "contourpy-1.1.0-cp39-cp39-win_amd64.whl", hash = "sha256:438ba416d02f82b692e371858143970ed2eb6337d9cdbbede0d8ad9f3d7dd17d"}, {file = "contourpy-1.1.0-pp38-pypy38_pp73-macosx_10_9_x86_64.whl", hash = "sha256:a698c6a7a432789e587168573a864a7ea374c6be8d4f31f9d87c001d5a843493"}, {file = "contourpy-1.1.0-pp38-pypy38_pp73-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:397b0ac8a12880412da3551a8cb5a187d3298a72802b45a3bd1805e204ad8439"}, @@ -808,7 +804,6 @@ files = [ {file = "greenlet-2.0.2-cp27-cp27m-win32.whl", hash = "sha256:6c3acb79b0bfd4fe733dff8bc62695283b57949ebcca05ae5c129eb606ff2d74"}, {file = "greenlet-2.0.2-cp27-cp27m-win_amd64.whl", hash = "sha256:283737e0da3f08bd637b5ad058507e578dd462db259f7f6e4c5c365ba4ee9343"}, {file = "greenlet-2.0.2-cp27-cp27mu-manylinux2010_x86_64.whl", hash = "sha256:d27ec7509b9c18b6d73f2f5ede2622441de812e7b1a80bbd446cb0633bd3d5ae"}, - {file = "greenlet-2.0.2-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:d967650d3f56af314b72df7089d96cda1083a7fc2da05b375d2bc48c82ab3f3c"}, {file = "greenlet-2.0.2-cp310-cp310-macosx_11_0_x86_64.whl", hash = "sha256:30bcf80dda7f15ac77ba5af2b961bdd9dbc77fd4ac6105cee85b0d0a5fcf74df"}, {file = "greenlet-2.0.2-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:26fbfce90728d82bc9e6c38ea4d038cba20b7faf8a0ca53a9c07b67318d46088"}, {file = "greenlet-2.0.2-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:9190f09060ea4debddd24665d6804b995a9c122ef5917ab26e1566dcc712ceeb"}, @@ -817,7 +812,6 @@ files = [ {file = "greenlet-2.0.2-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:76ae285c8104046b3a7f06b42f29c7b73f77683df18c49ab5af7983994c2dd91"}, {file = "greenlet-2.0.2-cp310-cp310-win_amd64.whl", hash = "sha256:2d4686f195e32d36b4d7cf2d166857dbd0ee9f3d20ae349b6bf8afc8485b3645"}, {file = "greenlet-2.0.2-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:c4302695ad8027363e96311df24ee28978162cdcdd2006476c43970b384a244c"}, - {file = "greenlet-2.0.2-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:d4606a527e30548153be1a9f155f4e283d109ffba663a15856089fb55f933e47"}, {file = "greenlet-2.0.2-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:c48f54ef8e05f04d6eff74b8233f6063cb1ed960243eacc474ee73a2ea8573ca"}, {file = "greenlet-2.0.2-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:a1846f1b999e78e13837c93c778dcfc3365902cfb8d1bdb7dd73ead37059f0d0"}, {file = "greenlet-2.0.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:3a06ad5312349fec0ab944664b01d26f8d1f05009566339ac6f63f56589bc1a2"}, @@ -847,7 +841,6 @@ files = [ {file = "greenlet-2.0.2-cp37-cp37m-win32.whl", hash = "sha256:3f6ea9bd35eb450837a3d80e77b517ea5bc56b4647f5502cd28de13675ee12f7"}, {file = "greenlet-2.0.2-cp37-cp37m-win_amd64.whl", hash = "sha256:7492e2b7bd7c9b9916388d9df23fa49d9b88ac0640db0a5b4ecc2b653bf451e3"}, {file = "greenlet-2.0.2-cp38-cp38-macosx_10_15_x86_64.whl", hash = "sha256:b864ba53912b6c3ab6bcb2beb19f19edd01a6bfcbdfe1f37ddd1778abfe75a30"}, - {file = "greenlet-2.0.2-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:1087300cf9700bbf455b1b97e24db18f2f77b55302a68272c56209d5587c12d1"}, {file = "greenlet-2.0.2-cp38-cp38-manylinux2010_x86_64.whl", hash = "sha256:ba2956617f1c42598a308a84c6cf021a90ff3862eddafd20c3333d50f0edb45b"}, {file = "greenlet-2.0.2-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fc3a569657468b6f3fb60587e48356fe512c1754ca05a564f11366ac9e306526"}, {file = "greenlet-2.0.2-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:8eab883b3b2a38cc1e050819ef06a7e6344d4a990d24d45bc6f2cf959045a45b"}, @@ -856,7 +849,6 @@ files = [ {file = "greenlet-2.0.2-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:b0ef99cdbe2b682b9ccbb964743a6aca37905fda5e0452e5ee239b1654d37f2a"}, {file = "greenlet-2.0.2-cp38-cp38-win32.whl", hash = "sha256:b80f600eddddce72320dbbc8e3784d16bd3fb7b517e82476d8da921f27d4b249"}, {file = "greenlet-2.0.2-cp38-cp38-win_amd64.whl", hash = "sha256:4d2e11331fc0c02b6e84b0d28ece3a36e0548ee1a1ce9ddde03752d9b79bba40"}, - {file = "greenlet-2.0.2-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:8512a0c38cfd4e66a858ddd1b17705587900dd760c6003998e9472b77b56d417"}, {file = "greenlet-2.0.2-cp39-cp39-macosx_11_0_x86_64.whl", hash = "sha256:88d9ab96491d38a5ab7c56dd7a3cc37d83336ecc564e4e8816dbed12e5aaefc8"}, {file = "greenlet-2.0.2-cp39-cp39-manylinux2010_x86_64.whl", hash = "sha256:561091a7be172ab497a3527602d467e2b3fbe75f9e783d8b8ce403fa414f71a6"}, {file = "greenlet-2.0.2-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:971ce5e14dc5e73715755d0ca2975ac88cfdaefcaab078a284fea6cfabf866df"}, @@ -1405,16 +1397,6 @@ files = [ {file = "MarkupSafe-2.1.3-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:5bbe06f8eeafd38e5d0a4894ffec89378b6c6a625ff57e3028921f8ff59318ac"}, {file = "MarkupSafe-2.1.3-cp311-cp311-win32.whl", hash = "sha256:dd15ff04ffd7e05ffcb7fe79f1b98041b8ea30ae9234aed2a9168b5797c3effb"}, {file = "MarkupSafe-2.1.3-cp311-cp311-win_amd64.whl", hash = "sha256:134da1eca9ec0ae528110ccc9e48041e0828d79f24121a1a146161103c76e686"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:f698de3fd0c4e6972b92290a45bd9b1536bffe8c6759c62471efaa8acb4c37bc"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:aa57bd9cf8ae831a362185ee444e15a93ecb2e344c8e52e4d721ea3ab6ef1823"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ffcc3f7c66b5f5b7931a5aa68fc9cecc51e685ef90282f4a82f0f5e9b704ad11"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:47d4f1c5f80fc62fdd7777d0d40a2e9dda0a05883ab11374334f6c4de38adffd"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-manylinux_2_5_i686.manylinux1_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:1f67c7038d560d92149c060157d623c542173016c4babc0c1913cca0564b9939"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-musllinux_1_1_aarch64.whl", hash = "sha256:9aad3c1755095ce347e26488214ef77e0485a3c34a50c5a5e2471dff60b9dd9c"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:14ff806850827afd6b07a5f32bd917fb7f45b046ba40c57abdb636674a8b559c"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:8f9293864fe09b8149f0cc42ce56e3f0e54de883a9de90cd427f191c346eb2e1"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-win32.whl", hash = "sha256:715d3562f79d540f251b99ebd6d8baa547118974341db04f5ad06d5ea3eb8007"}, - {file = "MarkupSafe-2.1.3-cp312-cp312-win_amd64.whl", hash = "sha256:1b8dd8c3fd14349433c79fa8abeb573a55fc0fdd769133baac1f5e07abf54aeb"}, {file = "MarkupSafe-2.1.3-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:8e254ae696c88d98da6555f5ace2279cf7cd5b3f52be2b5cf97feafe883b58d2"}, {file = "MarkupSafe-2.1.3-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:cb0932dc158471523c9637e807d9bfb93e06a95cbf010f1a38b98623b929ef2b"}, {file = "MarkupSafe-2.1.3-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9402b03f1a1b4dc4c19845e5c749e3ab82d5078d16a2a4c2cd2df62d57bb0707"}, @@ -2314,7 +2296,6 @@ files = [ {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:69b023b2b4daa7548bcfbd4aa3da05b3a74b772db9e23b982788168117739938"}, {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:81e0b275a9ecc9c0c0c07b4b90ba548307583c125f54d5b6946cfee6360c733d"}, {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ba336e390cd8e4d1739f42dfe9bb83a3cc2e80f567d8805e11b46f4a943f5515"}, - {file = "PyYAML-6.0.1-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:326c013efe8048858a6d312ddd31d56e468118ad4cdeda36c719bf5bb6192290"}, {file = "PyYAML-6.0.1-cp310-cp310-win32.whl", hash = "sha256:bd4af7373a854424dabd882decdc5579653d7868b8fb26dc7d0e99f823aa5924"}, {file = "PyYAML-6.0.1-cp310-cp310-win_amd64.whl", hash = "sha256:fd1592b3fdf65fff2ad0004b5e363300ef59ced41c2e6b3a99d4089fa8c5435d"}, {file = "PyYAML-6.0.1-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:6965a7bc3cf88e5a1c3bd2e0b5c22f8d677dc88a455344035f03399034eb3007"}, @@ -2322,15 +2303,8 @@ files = [ {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:42f8152b8dbc4fe7d96729ec2b99c7097d656dc1213a3229ca5383f973a5ed6d"}, {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:062582fca9fabdd2c8b54a3ef1c978d786e0f6b3a1510e0ac93ef59e0ddae2bc"}, {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d2b04aac4d386b172d5b9692e2d2da8de7bfb6c387fa4f801fbf6fb2e6ba4673"}, - {file = "PyYAML-6.0.1-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:e7d73685e87afe9f3b36c799222440d6cf362062f78be1013661b00c5c6f678b"}, {file = "PyYAML-6.0.1-cp311-cp311-win32.whl", hash = "sha256:1635fd110e8d85d55237ab316b5b011de701ea0f29d07611174a1b42f1444741"}, {file = "PyYAML-6.0.1-cp311-cp311-win_amd64.whl", hash = "sha256:bf07ee2fef7014951eeb99f56f39c9bb4af143d8aa3c21b1677805985307da34"}, - {file = "PyYAML-6.0.1-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:855fb52b0dc35af121542a76b9a84f8d1cd886ea97c84703eaa6d88e37a2ad28"}, - {file = "PyYAML-6.0.1-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:40df9b996c2b73138957fe23a16a4f0ba614f4c0efce1e9406a184b6d07fa3a9"}, - {file = "PyYAML-6.0.1-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:6c22bec3fbe2524cde73d7ada88f6566758a8f7227bfbf93a408a9d86bcc12a0"}, - {file = "PyYAML-6.0.1-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:8d4e9c88387b0f5c7d5f281e55304de64cf7f9c0021a3525bd3b1c542da3b0e4"}, - {file = "PyYAML-6.0.1-cp312-cp312-win32.whl", hash = "sha256:d483d2cdf104e7c9fa60c544d92981f12ad66a457afae824d146093b8c294c54"}, - {file = "PyYAML-6.0.1-cp312-cp312-win_amd64.whl", hash = "sha256:0d3304d8c0adc42be59c5f8a4d9e3d7379e6955ad754aa9d6ab7a398b59dd1df"}, {file = "PyYAML-6.0.1-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:50550eb667afee136e9a77d6dc71ae76a44df8b3e51e41b77f6de2932bfe0f47"}, {file = "PyYAML-6.0.1-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:1fe35611261b29bd1de0070f0b2f47cb6ff71fa6595c077e42bd0c419fa27b98"}, {file = "PyYAML-6.0.1-cp36-cp36m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:704219a11b772aea0d8ecd7058d0082713c3562b4e271b849ad7dc4a5c90c13c"}, @@ -2347,7 +2321,6 @@ files = [ {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a0cd17c15d3bb3fa06978b4e8958dcdc6e0174ccea823003a106c7d4d7899ac5"}, {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:28c119d996beec18c05208a8bd78cbe4007878c6dd15091efb73a30e90539696"}, {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7e07cbde391ba96ab58e532ff4803f79c4129397514e1413a7dc761ccd755735"}, - {file = "PyYAML-6.0.1-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:49a183be227561de579b4a36efbb21b3eab9651dd81b1858589f796549873dd6"}, {file = "PyYAML-6.0.1-cp38-cp38-win32.whl", hash = "sha256:184c5108a2aca3c5b3d3bf9395d50893a7ab82a38004c8f61c258d4428e80206"}, {file = "PyYAML-6.0.1-cp38-cp38-win_amd64.whl", hash = "sha256:1e2722cc9fbb45d9b87631ac70924c11d3a401b2d7f410cc0e3bbf249f2dca62"}, {file = "PyYAML-6.0.1-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:9eb6caa9a297fc2c2fb8862bc5370d0303ddba53ba97e71f08023b6cd73d16a8"}, @@ -2355,7 +2328,6 @@ files = [ {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5773183b6446b2c99bb77e77595dd486303b4faab2b086e7b17bc6bef28865f6"}, {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:b786eecbdf8499b9ca1d697215862083bd6d2a99965554781d0d8d1ad31e13a0"}, {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bc1bf2925a1ecd43da378f4db9e4f799775d6367bdb94671027b73b393a7c42c"}, - {file = "PyYAML-6.0.1-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:04ac92ad1925b2cff1db0cfebffb6ffc43457495c9b3c39d3fcae417d7125dc5"}, {file = "PyYAML-6.0.1-cp39-cp39-win32.whl", hash = "sha256:faca3bdcf85b2fc05d06ff3fbc1f83e1391b3e724afa3feba7d13eeab355484c"}, {file = "PyYAML-6.0.1-cp39-cp39-win_amd64.whl", hash = "sha256:510c9deebc5c0225e8c96813043e62b680ba2f9c50a08d3724c7f28a747d1486"}, {file = "PyYAML-6.0.1.tar.gz", hash = "sha256:bfdf460b1736c775f2ba9f6a92bca30bc2095067b8a9d77876d1fad6cc3b4a43"}, @@ -2884,7 +2856,6 @@ files = [ {file = "SQLAlchemy-1.4.49-cp27-cp27mu-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:03db81b89fe7ef3857b4a00b63dedd632d6183d4ea5a31c5d8a92e000a41fc71"}, {file = "SQLAlchemy-1.4.49-cp310-cp310-macosx_11_0_x86_64.whl", hash = "sha256:95b9df9afd680b7a3b13b38adf6e3a38995da5e162cc7524ef08e3be4e5ed3e1"}, {file = "SQLAlchemy-1.4.49-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a63e43bf3f668c11bb0444ce6e809c1227b8f067ca1068898f3008a273f52b09"}, - {file = "SQLAlchemy-1.4.49-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ca46de16650d143a928d10842939dab208e8d8c3a9a8757600cae9b7c579c5cd"}, {file = "SQLAlchemy-1.4.49-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:f835c050ebaa4e48b18403bed2c0fda986525896efd76c245bdd4db995e51a4c"}, {file = "SQLAlchemy-1.4.49-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9c21b172dfb22e0db303ff6419451f0cac891d2e911bb9fbf8003d717f1bcf91"}, {file = "SQLAlchemy-1.4.49-cp310-cp310-win32.whl", hash = "sha256:5fb1ebdfc8373b5a291485757bd6431de8d7ed42c27439f543c81f6c8febd729"}, @@ -2894,35 +2865,26 @@ files = [ {file = "SQLAlchemy-1.4.49-cp311-cp311-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:5debe7d49b8acf1f3035317e63d9ec8d5e4d904c6e75a2a9246a119f5f2fdf3d"}, {file = "SQLAlchemy-1.4.49-cp311-cp311-win32.whl", hash = "sha256:82b08e82da3756765c2e75f327b9bf6b0f043c9c3925fb95fb51e1567fa4ee87"}, {file = "SQLAlchemy-1.4.49-cp311-cp311-win_amd64.whl", hash = "sha256:171e04eeb5d1c0d96a544caf982621a1711d078dbc5c96f11d6469169bd003f1"}, - {file = "SQLAlchemy-1.4.49-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:f23755c384c2969ca2f7667a83f7c5648fcf8b62a3f2bbd883d805454964a800"}, - {file = "SQLAlchemy-1.4.49-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:8396e896e08e37032e87e7fbf4a15f431aa878c286dc7f79e616c2feacdb366c"}, - {file = "SQLAlchemy-1.4.49-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:66da9627cfcc43bbdebd47bfe0145bb662041472393c03b7802253993b6b7c90"}, - {file = "SQLAlchemy-1.4.49-cp312-cp312-win32.whl", hash = "sha256:9a06e046ffeb8a484279e54bda0a5abfd9675f594a2e38ef3133d7e4d75b6214"}, - {file = "SQLAlchemy-1.4.49-cp312-cp312-win_amd64.whl", hash = "sha256:7cf8b90ad84ad3a45098b1c9f56f2b161601e4670827d6b892ea0e884569bd1d"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-macosx_10_14_x86_64.whl", hash = "sha256:36e58f8c4fe43984384e3fbe6341ac99b6b4e083de2fe838f0fdb91cebe9e9cb"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b31e67ff419013f99ad6f8fc73ee19ea31585e1e9fe773744c0f3ce58c039c30"}, - {file = "SQLAlchemy-1.4.49-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ebc22807a7e161c0d8f3da34018ab7c97ef6223578fcdd99b1d3e7ed1100a5db"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:c14b29d9e1529f99efd550cd04dbb6db6ba5d690abb96d52de2bff4ed518bc95"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:c40f3470e084d31247aea228aa1c39bbc0904c2b9ccbf5d3cfa2ea2dac06f26d"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-win32.whl", hash = "sha256:706bfa02157b97c136547c406f263e4c6274a7b061b3eb9742915dd774bbc264"}, {file = "SQLAlchemy-1.4.49-cp36-cp36m-win_amd64.whl", hash = "sha256:a7f7b5c07ae5c0cfd24c2db86071fb2a3d947da7bd487e359cc91e67ac1c6d2e"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-macosx_11_0_x86_64.whl", hash = "sha256:4afbbf5ef41ac18e02c8dc1f86c04b22b7a2125f2a030e25bbb4aff31abb224b"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:24e300c0c2147484a002b175f4e1361f102e82c345bf263242f0449672a4bccf"}, - {file = "SQLAlchemy-1.4.49-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:393cd06c3b00b57f5421e2133e088df9cabcececcea180327e43b937b5a7caa5"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:201de072b818f8ad55c80d18d1a788729cccf9be6d9dc3b9d8613b053cd4836d"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7653ed6817c710d0c95558232aba799307d14ae084cc9b1f4c389157ec50df5c"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-win32.whl", hash = "sha256:647e0b309cb4512b1f1b78471fdaf72921b6fa6e750b9f891e09c6e2f0e5326f"}, {file = "SQLAlchemy-1.4.49-cp37-cp37m-win_amd64.whl", hash = "sha256:ab73ed1a05ff539afc4a7f8cf371764cdf79768ecb7d2ec691e3ff89abbc541e"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-macosx_11_0_x86_64.whl", hash = "sha256:37ce517c011560d68f1ffb28af65d7e06f873f191eb3a73af5671e9c3fada08a"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a1878ce508edea4a879015ab5215546c444233881301e97ca16fe251e89f1c55"}, - {file = "SQLAlchemy-1.4.49-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:95ab792ca493891d7a45a077e35b418f68435efb3e1706cb8155e20e86a9013c"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:0e8e608983e6f85d0852ca61f97e521b62e67969e6e640fe6c6b575d4db68557"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ccf956da45290df6e809ea12c54c02ace7f8ff4d765d6d3dfb3655ee876ce58d"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-win32.whl", hash = "sha256:f167c8175ab908ce48bd6550679cc6ea20ae169379e73c7720a28f89e53aa532"}, {file = "SQLAlchemy-1.4.49-cp38-cp38-win_amd64.whl", hash = "sha256:45806315aae81a0c202752558f0df52b42d11dd7ba0097bf71e253b4215f34f4"}, {file = "SQLAlchemy-1.4.49-cp39-cp39-macosx_11_0_x86_64.whl", hash = "sha256:b6d0c4b15d65087738a6e22e0ff461b407533ff65a73b818089efc8eb2b3e1de"}, {file = "SQLAlchemy-1.4.49-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a843e34abfd4c797018fd8d00ffffa99fd5184c421f190b6ca99def4087689bd"}, - {file = "SQLAlchemy-1.4.49-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:738d7321212941ab19ba2acf02a68b8ee64987b248ffa2101630e8fccb549e0d"}, {file = "SQLAlchemy-1.4.49-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:1c890421651b45a681181301b3497e4d57c0d01dc001e10438a40e9a9c25ee77"}, {file = "SQLAlchemy-1.4.49-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d26f280b8f0a8f497bc10573849ad6dc62e671d2468826e5c748d04ed9e670d5"}, {file = "SQLAlchemy-1.4.49-cp39-cp39-win32.whl", hash = "sha256:ec2268de67f73b43320383947e74700e95c6770d0c68c4e615e9897e46296294"}, @@ -3234,4 +3196,4 @@ mark = ["banana-hep", "matplotlib", "pandas", "sqlalchemy"] [metadata] lock-version = "2.0" python-versions = "^3.8,<3.12" -content-hash = "1b953c2ab9f1d4629f923a290a9db0b092ce2ff711b2e6009372cf1d602d5da3" +content-hash = "c47f4aa70036a1a754fc2528a335539342d0919a66dfe1051fe63cbce01338a2" diff --git a/pyproject.toml b/pyproject.toml index 31be6cb0c..d02ee9f80 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,7 +44,7 @@ PyYAML = "^6.0" lz4 = "^4.0.2" numba = "^0.57.0" # ekomark -banana-hep = { version = "^0.6.9", optional = true } +banana-hep = { version = "^0.6.12", optional = true } sqlalchemy = { version = "^1.4.21", optional = true } pandas = { version = "^1.3.0", optional = true } matplotlib = { version = "^3.5.1", optional = true } diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 29e67c19c..242fee66e 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -209,6 +209,7 @@ def quad_ker( n3lo_ad_variation, is_polarized, is_time_like, + use_fhmruvv, ): """Raw evolution kernel inside quad. @@ -255,11 +256,13 @@ def quad_ker( is_threshold : boolean is this an intermediate threshold operator? n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` is_polarized : boolean is polarized evolution ? is_time_like : boolean is time-like evolution ? + use_fhmruvv : bool + if True use the |FHMRUVV| |N3LO| anomalous dimension Returns ------- @@ -288,6 +291,7 @@ def quad_ker( is_polarized, is_time_like, n3lo_ad_variation, + use_fhmruvv, ) else: ker = quad_ker_qed( @@ -308,6 +312,7 @@ def quad_ker( sv_mode, is_threshold, n3lo_ad_variation, + use_fhmruvv, ) # recombine everything @@ -332,6 +337,7 @@ def quad_ker_qcd( is_polarized, is_time_like, n3lo_ad_variation, + use_fhmruvv, ): """Raw evolution kernel inside quad. @@ -364,7 +370,9 @@ def quad_ker_qcd( is_threshold : boolean is this an itermediate threshold operator? n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv : bool + if True use the |FHMRUVV| |N3LO| anomalous dimensions Returns ------- @@ -383,7 +391,7 @@ def quad_ker_qcd( gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf) else: gamma_singlet = ad_us.gamma_singlet( - order, ker_base.n, nf, n3lo_ad_variation + order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: @@ -416,7 +424,9 @@ def quad_ker_qcd( if is_time_like: gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf) else: - gamma_ns = ad_us.gamma_ns(order, mode0, ker_base.n, nf) + gamma_ns = ad_us.gamma_ns( + order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) if sv_mode == sv.Modes.exponentiated: gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) ker = ns.dispatcher( @@ -452,6 +462,7 @@ def quad_ker_qed( sv_mode, is_threshold, n3lo_ad_variation, + use_fhmruvv, ): """Raw evolution kernel inside quad. @@ -492,7 +503,9 @@ def quad_ker_qed( is_threshold : boolean is this an itermediate threshold operator? n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv : bool + if True use the |FHMRUVV| |N3LO| anomalous dimensions Returns ------- @@ -501,7 +514,9 @@ def quad_ker_qed( """ # compute the actual evolution kernel for QEDxQCD if ker_base.is_QEDsinglet: - gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation) + gamma_s = ad_us.gamma_singlet_qed( + order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_s = sv.exponentiated.gamma_variation_qed( @@ -529,7 +544,9 @@ def quad_ker_qed( ) @ np.ascontiguousarray(ker) ker = select_QEDsinglet_element(ker, mode0, mode1) elif ker_base.is_QEDvalence: - gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) + gamma_v = ad_us.gamma_valence_qed( + order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_v = sv.exponentiated.gamma_variation_qed( @@ -554,7 +571,9 @@ def quad_ker_qed( ) @ np.ascontiguousarray(ker) ker = select_QEDvalence_element(ker, mode0, mode1) else: - gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) + gamma_ns = ad_us.gamma_ns_qed( + order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_ns = sv.exponentiated.gamma_variation_qed( @@ -825,6 +844,7 @@ def quad_ker(self, label, logx, areas): n3lo_ad_variation=self.config["n3lo_ad_variation"], is_polarized=self.config["polarized"], is_time_like=self.config["time_like"], + use_fhmruvv=self.config["use_fhmruvv"], ) def initialize_op_members(self): @@ -946,11 +966,12 @@ def compute(self): self.a_em[1], ) logger.info( - "%s: order: (%d, %d), solution strategy: %s", + "%s: order: (%d, %d), solution strategy: %s, use fhmruvv: %s", self.log_label, self.order[0], self.order[1], self.config["method"], + self.config["use_fhmruvv"], ) self.integrate() diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index 3e066a6de..d1ddb23bc 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -60,6 +60,7 @@ def __init__( atlas: Atlas, couplings: Couplings, interpol_dispatcher: InterpolatorDispatcher, + use_fhmruvv: bool, ): # check config = {} @@ -69,6 +70,7 @@ def __init__( config["HQ"] = mass_scheme config["ModSV"] = configs.scvar_method config["n3lo_ad_variation"] = n3lo_ad_variation + config["use_fhmruvv"] = use_fhmruvv for i, q in enumerate("cbt"): config[f"m{q}"] = masses[i] diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 63aa5c215..755fcaf72 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -46,7 +46,9 @@ class TheoryCard(DictLike): xif: float """Ratio between factorization scale and process scale.""" n3lo_ad_variation: N3LOAdVariation - """|N3LO| anomalous dimension variation: ``(gg, gq, qg,qq)``.""" + """|N3LO| anomalous dimension variation: ``(gg, gq, qg, qq, nsp, nsm, nsv)``.""" + use_fhmruvv: Optional[bool] = False + """If True use the |FHMRUVV| |N3LO| anomalous dimensions""" matching_order: Optional[Order] = None """Matching conditions perturbative order tuple, ``(QCD, QED)``.""" @@ -217,10 +219,11 @@ def new_theory(self): raise ValueError(f"Unknown mass scheme '{old['HQ']}'") new["xif"] = old["XIF"] - new["n3lo_ad_variation"] = old.get("n3lo_ad_variation", (0, 0, 0, 0)) + + new["n3lo_ad_variation"] = old.get("n3lo_ad_variation", (0, 0, 0, 0, 0, 0, 0)) # here PTO: 0 means truly LO, no QED matching is available so far. new["matching_order"] = old.get("PTO_matching", [old["PTO"], 0]) - + new["use_fhmruvv"] = old.get("use_fhmruvv", False) return TheoryCard.from_dict(new) @property diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py index 991ec2bd0..39659f121 100644 --- a/src/eko/runner/legacy.py +++ b/src/eko/runner/legacy.py @@ -77,6 +77,7 @@ def __init__( couplings=cs, interpol_dispatcher=bfd, n3lo_ad_variation=new_theory.n3lo_ad_variation, + use_fhmruvv=new_theory.use_fhmruvv, matching_order=new_theory.matching_order, ) diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py index f3d93a558..f003d955d 100644 --- a/src/eko/runner/parts.py +++ b/src/eko/runner/parts.py @@ -84,6 +84,7 @@ def evolve_configs(eko: EKO) -> dict: n_integration_cores=ocard.configs.n_integration_cores, ModSV=ocard.configs.scvar_method, n3lo_ad_variation=tcard.n3lo_ad_variation, + use_fhmruvv=tcard.use_fhmruvv, # Here order is shifted by one, no QED matching is available so far. matching_order=tcard.matching_order, ) diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 559790d30..c7536bc50 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -26,8 +26,9 @@ matching_ratios=[1.0, 1.0, 1.0], ), xif=1.0, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), matching_order=[0, 0], + use_fhmruvv=False, ) _operator = dict( diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 6358a5af8..32687e489 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -25,7 +25,7 @@ @nb.njit(cache=True) -def gamma_ns(order, mode, n, nf): +def gamma_ns(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False): r"""Compute the tower of the non-singlet anomalous dimensions. Parameters @@ -38,6 +38,10 @@ def gamma_ns(order, mode, n, nf): Mellin variable nf : int Number of active flavors + n3lo_ad_variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -70,18 +74,32 @@ def gamma_ns(order, mode, n, nf): gamma_ns[2] = gamma_ns_2 # N3LO if order[0] >= 4: - if mode == 10101: - gamma_ns_3 = as4.gamma_nsp(n, nf, cache) - elif mode == 10201: - gamma_ns_3 = as4.gamma_nsm(n, nf, cache) - elif mode == 10200: - gamma_ns_3 = as4.gamma_nsv(n, nf, cache) + if use_fhmruvv: + if mode == 10101: + gamma_ns_3 = as4.fhmruvv.gamma_nsp( + n, nf, cache, variation=n3lo_ad_variation[4] + ) + elif mode == 10201: + gamma_ns_3 = as4.fhmruvv.gamma_nsm( + n, nf, cache, variation=n3lo_ad_variation[5] + ) + elif mode == 10200: + gamma_ns_3 = as4.fhmruvv.gamma_nsv( + n, nf, cache, variation=n3lo_ad_variation[6] + ) + else: + if mode == 10101: + gamma_ns_3 = as4.gamma_nsp(n, nf, cache) + elif mode == 10201: + gamma_ns_3 = as4.gamma_nsm(n, nf, cache) + elif mode == 10200: + gamma_ns_3 = as4.gamma_nsv(n, nf, cache) gamma_ns[3] = gamma_ns_3 return gamma_ns @nb.njit(cache=True) -def gamma_singlet(order, n, nf, n3lo_ad_variation): +def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): r"""Compute the tower of the singlet anomalous dimensions matrices. Parameters @@ -93,7 +111,9 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation): nf : int Number of active flavors n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -109,25 +129,32 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation): if order[0] >= 3: gamma_s[2] = as3.gamma_singlet(n, nf, cache) if order[0] >= 4: - gamma_s[3] = as4.gamma_singlet(n, nf, cache, n3lo_ad_variation) + if use_fhmruvv: + gamma_s[3] = as4.fhmruvv.gamma_singlet(n, nf, cache, n3lo_ad_variation) + else: + gamma_s[3] = as4.gamma_singlet(n, nf, cache, n3lo_ad_variation) return gamma_s @nb.njit(cache=True) -def gamma_ns_qed(order, mode, n, nf): +def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED non-singlet anomalous dimensions. Parameters ---------- - order : tuple(int,int) - perturbative orders - mode : 10102 | 10103 | 10202 | 10203 - sector identifier - n : complex - Mellin variable - nf : int - Number of active flavors + order : tuple(int,int) + perturbative orders + mode : 10102 | 10103 | 10202 | 10203 + sector identifier + n : complex + Mellin variable + nf : int + Number of active flavors + n3lo_ad_variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -156,10 +183,20 @@ def gamma_ns_qed(order, mode, n, nf): elif mode in [10202, 10203]: gamma_ns[3, 0] = as3.gamma_nsm(n, nf, cache) if order[0] >= 4: - if mode in [10102, 10103]: - gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache) - elif mode in [10202, 10203]: - gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache) + if use_fhmruvv: + if mode in [10102, 10103]: + gamma_ns[4, 0] = as4.fhmruvv.gamma_nsp( + n, nf, cache, n3lo_ad_variation[4] + ) + elif mode in [10202, 10203]: + gamma_ns[4, 0] = as4.fhmruvv.gamma_nsm( + n, nf, cache, n3lo_ad_variation[5] + ) + else: + if mode in [10102, 10103]: + gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache) + elif mode in [10202, 10203]: + gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache) return gamma_ns @@ -255,7 +292,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache): @nb.njit(cache=True) -def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): +def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED singlet anomalous dimensions matrices. @@ -268,7 +305,9 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): nf : int Number of active flavors n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -288,12 +327,17 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): if order[0] >= 3: gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, cache) if order[0] >= 4: - gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation) + if use_fhmruvv: + gamma_s[4, 0] = as4.fhmruvv.gamma_singlet_qed( + n, nf, cache, n3lo_ad_variation + ) + else: + gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation) return gamma_s @nb.njit(cache=True) -def gamma_valence_qed(order, n, nf): +def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED valence anomalous dimensions matrices. @@ -305,6 +349,11 @@ def gamma_valence_qed(order, n, nf): Mellin variable nf : int Number of active flavors + n3lo_ad_variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions + Returns ------- @@ -324,5 +373,10 @@ def gamma_valence_qed(order, n, nf): if order[0] >= 3: gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, cache) if order[0] >= 4: - gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache) + if use_fhmruvv: + gamma_v[4, 0] = as4.fhmruvv.gamma_valence_qed( + n, nf, cache, n3lo_ad_variation + ) + else: + gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache) return gamma_v diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py index 56c2c16ff..9b3b37517 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py @@ -6,6 +6,7 @@ import numba as nb import numpy as np +from . import fhmruvv from .ggg import gamma_gg from .ggq import gamma_gq from .gnsm import gamma_nsm @@ -34,7 +35,7 @@ def gamma_singlet(N, nf, cache, variation): cache: numpy.ndarray Harmonic sum cache variation : tuple - |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` Returns ------- diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py new file mode 100644 index 000000000..e92a55ffa --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py @@ -0,0 +1,143 @@ +"""The |FHMRUVV| |N3LO| Altarelli-Parisi splitting kernels approximations. + +Authors follow Pegasus convention and so there is an additional global minus sign with respect to our conventions. +""" +import numba as nb +import numpy as np + +from .ggg import gamma_gg +from .ggq import gamma_gq +from .gnsm import gamma_nsm +from .gnsp import gamma_nsp +from .gnsv import gamma_nsv +from .gps import gamma_ps +from .gqg import gamma_qg + + +@nb.njit(cache=True) +def gamma_singlet(N, nf, cache, variation): + r"""Compute the |N3LO| singlet anomalous dimension matrix. + + .. math:: + \gamma_S^{(3)} = \left(\begin{array}{cc} + \gamma_{qq}^{(3)} & \gamma_{qg}^{(3)}\\ + \gamma_{gq}^{(3)} & \gamma_{gg}^{(3)} + \end{array}\right) + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` + + Returns + ------- + numpy.ndarray + |N3LO| singlet anomalous dimension matrix + :math:`\gamma_{S}^{(3)}(N)` + + """ + gamma_qq = gamma_nsp(N, nf, cache, variation[3]) + gamma_ps( + N, nf, cache, variation[3] + ) + gamma_S_0 = np.array( + [ + [gamma_qq, gamma_qg(N, nf, cache, variation[2])], + [ + gamma_gq(N, nf, cache, variation[1]), + gamma_gg(N, nf, cache, variation[0]), + ], + ], + np.complex_, + ) + return gamma_S_0 + + +@nb.njit(cache=True) +def gamma_singlet_qed(N, nf, cache, variation): + r"""Compute the leading-order singlet anomalous dimension matrix for the unified evolution basis. + + .. math:: + \\gamma_S^{(3,0)} = \\left(\begin{array}{cccc} + \\gamma_{gg}^{(3,0)} & 0 & \\gamma_{gq}^{(3,0)} & 0\\ + 0 & 0 & 0 & 0 \\ + \\gamma_{qg}^{(3,0)} & 0 & \\gamma_{qq}^{(3,0)} & 0 \\ + 0 & 0 & 0 & \\gamma_{qq}^{(3,0)} \\ + \\end{array}\right) + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` + + Returns + ------- + numpy.ndarray + Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)` + + """ + gamma_np_p = gamma_nsp(N, nf, cache, variation[3]) + gamma_qq = gamma_np_p + gamma_ps(N, nf, cache, variation[3]) + gamma_S = np.array( + [ + [ + gamma_gg(N, nf, cache, variation[0]), + 0.0 + 0.0j, + gamma_gq(N, nf, cache, variation[1]), + 0.0 + 0.0j, + ], + [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j], + [gamma_qg(N, nf, cache, variation[2]), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j], + [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_np_p], + ], + np.complex_, + ) + return gamma_S + + +@nb.njit(cache=True) +def gamma_valence_qed(N, nf, cache, variation): + r"""Compute the leading-order valence anomalous dimension matrix for the unified evolution basis. + + .. math:: + \\gamma_V^{(3,0)} = \\left(\begin{array}{cc} + \\gamma_{nsV}^{(3,0)} & 0\\ + 0 & \\gamma_{ns-}^{(3,0)} + \\end{array}\right) + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : tuple + |N3LO| anomalous dimension variation ``(nsm, nsv)`` + + Returns + ------- + numpy.ndarray + Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)` + + """ + gamma_V = np.array( + [ + [gamma_nsv(N, nf, cache, variation[-1]), 0.0], + [0.0, gamma_nsm(N, nf, cache, variation[-2])], + ], + np.complex_, + ) + return gamma_V diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py new file mode 100644 index 000000000..1d2a2555b --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py @@ -0,0 +1,130 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{gg}^{(3)}`.""" +import numba as nb + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm11, lm12m1, lm13m1 + + +@nb.njit(cache=True) +def gamma_gg(n, nf, cache, variation): + r"""Compute the |N3LO| gluon-gluon singlet anomalous dimension. + + The routine is taken from :cite:`Moch:2023tdj`. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| gluon-gluon singlet anomalous dimension + :math:`\gamma_{gg}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + + nf2 = nf * nf + nf3 = nf * nf2 + + # The known large-x coefficients [except delta(1-x)] + A4gluon = 40880.330 - 11714.246 * nf + 440.04876 * nf2 + 7.3627750 * nf3 + B4gluon = 68587.64 - 18143.983 * nf + 423.81135 * nf2 + 9.0672154 * 0.1 * nf3 + + Ccoeff = 8.5814120 * 10**4 - 1.3880515 * 10**4 * nf + 1.3511111 * 10**2 * nf2 + Dcoeff = 5.4482808 * 10**4 - 4.3411337 * 10**3 * nf - 2.1333333 * 10 * nf2 + + # The known coefficients of 1/x*ln^a x terms, a = 3,2 + bfkl0 = -8.308617314 * 10**3 + bfkl1 = -1.069119905 * 10**5 - 9.963830436 * 10**2 * nf + + # The resulting part of the function + P3gg01 = ( + +bfkl0 * (-(6 / (-1 + n) ** 4)) + + bfkl1 * 2 / (-1 + n) ** 3 + + A4gluon * (-S1) + + B4gluon + + Ccoeff * lm11(n, S1) + + Dcoeff * 1 / n + ) + + # The selected approximations for nf = 3, 4, 5 + if nf == 3: + P3ggApp1 = ( + P3gg01 + + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 345063.0 * 1 / ((-1 + n) * n) + + 86650.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + + 158160.0 * (-(1 / n**2)) + - 15741.0 * lm12m1(n, S1, S2) + - 9417.0 * lm13m1(n, S1, S2, S3) + ) + P3ggApp2 = ( + P3gg01 + + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 1265632.0 * 1 / ((-1 + n) * n) + - 656644.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + - 1352233.0 * (-(1 / n**2)) + + 203298.0 * lm12m1(n, S1, S2) + + 39112.0 * lm13m1(n, S1, S2, S3) + ) + elif nf == 4: + P3ggApp1 = ( + P3gg01 + + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 342625.0 * 1 / ((-1 + n) * n) + + 100372.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + + 189167.0 * (-(1 / n**2)) + - 29762.0 * lm12m1(n, S1, S2) + - 12102.0 * lm13m1(n, S1, S2, S3) + ) + P3ggApp2 = ( + P3gg01 + + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 1271540.0 * 1 / ((-1 + n) * n) + - 649661.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + - 1334919.0 * (-(1 / n**2)) + + 191263.0 * lm12m1(n, S1, S2) + + 36867.0 * lm13m1(n, S1, S2, S3) + ) + elif nf == 5: + P3ggApp1 = ( + P3gg01 + + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 337540.0 * 1 / ((-1 + n) * n) + + 119366.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + + 223769.0 * (-(1 / n**2)) + - 45129.0 * lm12m1(n, S1, S2) + - 15046.0 * lm13m1(n, S1, S2, S3) + ) + P3ggApp2 = ( + P3gg01 + + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2) + - 1274800.0 * 1 / ((-1 + n) * n) + - 637406.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n)) + - 1314010.0 * (-(1 / n**2)) + + 177882.0 * lm12m1(n, S1, S2) + + 34362.0 * lm13m1(n, S1, S2, S3) + ) + else: + raise NotImplementedError("nf=6 is not available at N3LO") + + # We return (for now) one of the two error-band representatives + # or the present best estimate, their average + if variation == 1: + P3GGA = P3ggApp1 + elif variation == 2: + P3GGA = P3ggApp2 + else: + P3GGA = 0.5 * (P3ggApp1 + P3ggApp2) + + return -P3GGA diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py new file mode 100644 index 000000000..97904f459 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py @@ -0,0 +1,124 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{gq}^{(3)}`.""" +import numba as nb + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm12, lm13, lm14, lm15 + + +@nb.njit(cache=True) +def gamma_gq(n, nf, cache, variation): + r"""Compute the |N3LO| gluon-quark singlet anomalous dimension. + + The routine is taken from :cite:`Moch:2023tdj`. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| gluon-quark singlet anomalous dimension + :math:`\gamma_{gq}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + S4 = c.get(c.S4, cache, n) + S5 = c.get(c.S5, cache, n) + nf2 = nf * nf + + # Known large-x coefficients + x1L5cff = 1.3443073 * 10 - 5.4869684 * 0.1 * nf + x1L4cff = 3.7539831 * 10**2 - 3.4494742 * 10 * nf + 8.7791495 * 0.1 * nf2 + + # Small-x, Casimir scaled from P_gg (approx. for bfkl1) + bfkl0 = -8.3086173 * 10**3 / 2.25 + bfkl1 = (-1.0691199 * 10**5 - nf * 9.9638304 * 10**2) / 2.25 + + # The resulting part of the function + P3GQ01 = ( + +bfkl0 * (-(6 / (-1 + n) ** 4)) + + bfkl1 * 2 / (-1 + n) ** 3 + + x1L4cff * lm14(n, S1, S2, S3, S4) + + x1L5cff * lm15(n, S1, S2, S3, S4, S5) + ) + + # The selected approximations for nf = 3, 4, 5 + if nf == 3: + P3gqApp1 = ( + P3GQ01 + + 3.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 161562.0 * 1 / ((-1 + n) * n) + + 36469.0 * 1 / n + + 72317.0 * (-(1 / n**2)) + - 3977.3 * lm12(n, S1, S2) + + 484.4 * lm13(n, S1, S2, S3) + ) + P3gqApp2 = ( + P3GQ01 + + 5.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 546482.0 * 1 / ((-1 + n) * n) + - 39464.0 * 1 / n + - 401000.0 * (-(1 / n**2)) + + 13270.0 * lm12(n, S1, S2) + + 3289.0 * lm13(n, S1, S2, S3) + ) + elif nf == 4: + P3gqApp1 = ( + P3GQ01 + + 3.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 158805.0 * 1 / ((-1 + n) * n) + + 35098.0 * 1 / n + + 87258.0 * (-(1 / n**2)) + - 4834.1 * lm12(n, S1, S2) + + 176.6 * lm13(n, S1, S2, S3) + ) + P3gqApp2 = ( + P3GQ01 + + 5.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 547215.0 * 1 / ((-1 + n) * n) + - 41523.0 * 1 / n + - 390350.0 * (-(1 / n**2)) + + 12571.0 * lm12(n, S1, S2) + + 3007.0 * lm13(n, S1, S2, S3) + ) + elif nf == 5: + P3gqApp1 = ( + P3GQ01 + + 3.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 154336.0 * 1 / ((-1 + n) * n) + + 33889.0 * 1 / n + + 103440.0 * (-(1 / n**2)) + - 5745.8 * lm12(n, S1, S2) + - 128.6 * lm13(n, S1, S2, S3) + ) + P3gqApp2 = ( + P3GQ01 + + 5.4 * bfkl1 * (-(1 / (-1 + n) ** 2)) + - 546236.0 * 1 / ((-1 + n) * n) + - 43421.0 * 1 / n + - 378460.0 * (-(1 / n**2)) + + 11816.0 * lm12(n, S1, S2) + + 2727.3 * lm13(n, S1, S2, S3) + ) + else: + raise NotImplementedError("nf=6 is not available at N3LO") + + # We return (for now) one of the two error-band representatives + # or the present best estimate, their average + if variation == 1: + P3GQA = P3gqApp1 + elif variation == 2: + P3GQA = P3gqApp2 + else: + P3GQA = 0.5 * (P3gqApp1 + P3gqApp2) + return -P3GQA diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsm.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsm.py new file mode 100644 index 000000000..5d80c9c75 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsm.py @@ -0,0 +1,221 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{ns,-}^{(3)}`.""" +import numba as nb + +from eko.constants import CF, zeta3 + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm11, lm11m1, lm12m1, lm13m1 + + +@nb.njit(cache=True) +def gamma_nsm(n, nf, cache, variation): + r"""Compute the |N3LO| valence-like non-singlet anomalous dimension. + + The routine is taken from :cite:`Moch:2017uml`. + + The :math:`nf^{0,1}` leading large-nc contributions and the :math:`nf^2` part are + high-accuracy (0.1% or better) parametrizations of the exact + results. The :math:`nf^3` expression is exact up to numerical truncations. + + The remaining :math:`nf^{0,1}` terms are approximations based on the first + eight odd moments together with small-x and large-x constraints. + The two sets spanning the error estimate are called via IMOD = 1 + and IMOD = 2. Any other value of IMOD invokes their average. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| valence-like non-singlet anomalous dimension + :math:`\gamma_{ns,-}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + S4 = c.get(c.S4, cache, n) + + # Leading large-n_c, nf^0 and nf^1, parametrized + P3NSA0 = ( + 360.0 / n**7 + - 1920.0 / n**6 + + 7147.812 / n**5 + - 17179.356 / n**4 + + 34241.9 / n**3 + - 51671.329999999994 / n**2 + + 19069.8 * lm11(n, S1) + - (491664.8019540468 / n) + - 4533.0 / (1 + n) ** 3 + - 11825.0 / (1 + n) ** 2 + + 129203.0 / (1 + n) + - 254965.0 / (2 + n) + + 83377.5 / (3 + n) + - 45750.0 / (4 + n) + + (49150.0 * (6.803662258392675 + n) * S1) / (n**2 * (1.0 + n)) + + (334400.0 * S2) / n + ) + P3NSA1 = ( + 160.0 / n**6 + - 864.0 / n**5 + + 2583.1848 / n**4 + - 5834.624 / n**3 + + 9239.374 / n**2 + - 3079.76 * lm11(n, S1) + - (114047.0 / n) + - 465.0 / (1 + n) ** 4 + - 1230.0 / (1 + n) ** 3 + + 7522.5 / (1 + n) ** 2 + + 55669.3 / (1 + n) + - 43057.8 / (2 + n) + + 13803.8 / (3 + n) + - 7896.0 / (4 + n) + - (120.0 * (-525.063 + n) * S1) / (n**2 * (1.0 + n)) + + (63007.5 * S2) / n + ) + + # Nonleading large-n_c, nf^0 and nf^1: two approximations + P3NMA01 = ( + 0.4964335 * (720 / n**7 - 720.0 / n**6) + - 13.5288 / n**4 + + 1618.07 / n**2 + - 2118.8669999999997 * lm11(n, S1) + + 31897.8 * lm11m1(n, S1) + + 4653.76 * lm12m1(n, S1, S2) + + 3902.3590000000004 / n + + 5992.88 / (1 + n) + + 19335.7 / (2 + n) + - 31321.4 / (3 + n) + ) + P3NMA02 = ( + +0.4964335 * (720 / n**7 - 2160.0 / n**6) + - 189.6138 / n**4 + + 3065.92 / n**3 + - 2118.8669999999997 * lm11(n, S1) + - 3997.39 * lm11m1(n, S1) + + 511.567 * lm13m1(n, S1, S2, S3) + - (2099.268 / n) + + 4043.59 / (1 + n) + - 19430.190000000002 / (2 + n) + + 15386.6 / (3 + n) + ) + + P3NMA11 = ( + +64.7083 / n**5 + - 254.024 / n**3 + + 337.931 * lm11(n, S1) + + 1856.63 * lm11m1(n, S1) + + 440.17 * lm12m1(n, S1, S2) + + 419.53485 / n + + 114.457 / (1 + n) + + 2341.816 / (2 + n) + - 2570.73 / (3 + n) + ) + + P3NMA12 = ( + -17.0616 / n**6 + - 19.53254 / n**3 + + 337.931 * lm11(n, S1) + - 1360.04 * lm11m1(n, S1) + + 38.7337 * lm13m1(n, S1, S2, S3) + - (367.64646999999997 / n) + + 335.995 / (1 + n) + - 1269.915 / (2 + n) + + 1605.91 / (3 + n) + ) + + # nf^2 (parametrized) and nf^3 (exact) + P3NSMA2 = -( + -193.84583328013258 + - 23.7037032 / n**5 + + 117.5967 / n**4 + - 256.5896 / n**3 + + 437.881 / n**2 + + 720.385709813466 / n + - 48.720000000000006 / (1 + n) ** 4 + + 189.51000000000002 / (1 + n) ** 3 + + 391.02500000000003 / (1 + n) ** 2 + + 367.4750000000001 / (1 + n) + + 404.47249999999997 / (2 + n) + - 2063.325 / ((1 + n) ** 2 * (2 + n)) + - (1375.55 * n) / ((1 + n) ** 2 * (2 + n)) + + 687.775 / ((1 + n) * (2 + n)) + - 81.71999999999998 / (3 + n) + + 114.9225 / (4 + n) + + 195.5772 * S1 + - (817.725 * S1) / n**2 + + (714.46361 * S1) / n + - (687.775 * S1) / (1 + n) + - (817.725 * S2) / n + ) + eta = 1 / n * 1 / (n + 1) + P3NSA3 = -CF * ( + -32 / 27 * zeta3 * eta + - 16 / 9 * zeta3 + - 16 / 27 * eta**4 + - 16 / 81 * eta**3 + + 80 / 27 * eta**2 + - 320 / 81 * eta + + 32 / 27 * 1 / (n + 1) ** 4 + + 128 / 27 * 1 / (n + 1) ** 2 + + 64 / 27 * S1 * zeta3 + - 32 / 81 * S1 + - 32 / 81 * S2 + - 160 / 81 * S3 + + 32 / 27 * S4 + + 131 / 81 + ) + + # Assembly regular piece. + P3NSMAI = P3NSA0 + nf * P3NSA1 + nf**3 * P3NSA3 + nf**2 * P3NSMA2 + if variation == 1: + P3NSMA = P3NSMAI + P3NMA01 + nf * P3NMA11 + elif variation == 2: + P3NSMA = P3NSMAI + P3NMA02 + nf * P3NMA12 + else: + P3NSMA = P3NSMAI + 0.5 * ((P3NMA01 + P3NMA02) + nf * (P3NMA11 + P3NMA12)) + + # The singular piece. + A4qI = ( + 2.120902 * 10**4 + - 5.179372 * 10**3 * nf + # + 1.955772 * 10**2 * nf**2 + # + 3.272344 * nf**3 + ) + A4ap1 = -511.228 + 7.08645 * nf + A4ap2 = -502.481 + 7.82077 * nf + D1 = 1 / n - S1 + if variation == 1: + P3NSMB = (A4qI + A4ap1) * D1 + elif variation == 2: + P3NSMB = (A4qI + A4ap2) * D1 + else: + P3NSMB = (A4qI + 0.5 * (A4ap1 + A4ap2)) * D1 + + # The local piece. + B4qI = ( + 2.579609 * 10**4 + + 0.08 + - (5.818637 * 10**3 + 0.97) * nf + # + (1.938554 * 10**2 + 0.0037) * nf**2 + # + 3.014982 * nf**3 + ) + B4ap1 = -2426.05 + 266.674 * nf - 0.05 * nf + B4ap2 = -2380.255 + 270.518 * nf - 0.05 * nf + if variation == 1: + P3NSMC = B4qI + B4ap1 + elif variation == 2: + P3NSMC = B4qI + B4ap2 + else: + P3NSMC = +B4qI + 0.5 * (B4ap1 + B4ap2) + + return -(P3NSMA + P3NSMB + P3NSMC) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsp.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsp.py new file mode 100644 index 000000000..de69d5313 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsp.py @@ -0,0 +1,228 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{ns,+}^{(3)}`.""" +import numba as nb + +from eko.constants import CF, zeta3 + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm11, lm11m1, lm12m1, lm13m1 + + +@nb.njit(cache=True) +def gamma_nsp(n, nf, cache, variation): + r"""Compute the |N3LO| singlet-like non-singlet anomalous dimension. + + The routine is taken from :cite:`Moch:2017uml`. + + The :math:`nf^{0,1}` leading large-nc contributions and the :math:`nf^2` part + are high-accuracy (0.1% or better) parametrizations of the exact + results. The :math:`nf^3` expression is exact up to numerical truncations. + + The remaining :math:`nf^{0,1}` terms are approximations based on the first + eight even moments together with small-x and large-x constraints. + The two sets spanning the error estimate are called via IMOD = 1 + and IMOD = 2. Any other value of IMOD invokes their average. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| singlet-like non-singlet anomalous dimension + :math:`\gamma_{ns,+}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + S4 = c.get(c.S4, cache, n) + + # Leading large-n_c, nf^0 and nf^1, parametrized + P3NSA0 = ( + 360.0 / n**7 + - 1920.0 / n**6 + + 7147.812 / n**5 + - 17179.356 / n**4 + + 34241.9 / n**3 + - 51671.329999999994 / n**2 + + 19069.8 * lm11(n, S1) + - (491664.8019540468 / n) + - 4533.0 / (1 + n) ** 3 + - 11825.0 / (1 + n) ** 2 + + 129203.0 / (1 + n) + - 254965.0 / (2 + n) + + 83377.5 / (3 + n) + - 45750.0 / (4 + n) + + (49150.0 * (6.803662258392675 + n) * S1) / (n**2 * (1.0 + n)) + + (334400.0 * S2) / n + ) + P3NSA1 = ( + 160.0 / n**6 + - 864.0 / n**5 + + 2583.1848 / n**4 + - 5834.624 / n**3 + + 9239.374 / n**2 + - 3079.76 * lm11(n, S1) + - (114047.0 / n) + - 465.0 / (1 + n) ** 4 + - 1230.0 / (1 + n) ** 3 + + 7522.5 / (1 + n) ** 2 + + 55669.3 / (1 + n) + - 43057.8 / (2 + n) + + 13803.8 / (3 + n) + - 7896.0 / (4 + n) + - (120.0 * (-525.063 + n) * S1) / (n**2 * (1.0 + n)) + + (63007.5 * S2) / n + ) + + # Nonleading large-n_c, nf^0 and nf^1: two approximations + P3NPA01 = ( + -(107.16 / n**7) + + 339.753 / n**6 + - 1341.01 / n**5 + + 2412.94 / n**4 + - 3678.88 / n**3 + - 2118.87 * lm11(n, S1) + - 1777.27 * lm12m1(n, S1, S2) + - 204.183 * lm13m1(n, S1, S2, S3) + + 1853.56 / n + - 8877.38 / (1 + n) + + 7393.83 / (2 + n) + - 2464.61 / (3 + n) + ) + P3NPA02 = ( + -(107.16 / n**7) + + 339.753 / n**6 + - 1341.01 / n**5 + + 379.152 / n**3 + - 1389.73 / n**2 + - 2118.87 * lm11(n, S1) + - 173.936 * lm12m1(n, S1, S2) + + 223.078 * lm13m1(n, S1, S2, S3) + - (2096.54 / n) + + 8698.39 / (1 + n) + - 19188.9 / (2 + n) + + 10490.5 / (3 + n) + ) + + P3NPA11 = ( + -(33.5802 / n**6) + + 111.802 / n**5 + + 50.772 / n**4 + - 118.608 / n**3 + + 337.931 * lm11(n, S1) + - 143.813 * lm11m1(n, S1) + - 18.8803 * lm13m1(n, S1, S2, S3) + + 304.82503 / n + - 1116.34 / (1 + n) + + 2187.58 / (2 + n) + - 1071.24 / (3 + n) + ) + P3NPA12 = ( + -(33.5802 / n**6) + + 111.802 / n**5 + - 204.341 / n**4 + + 267.404 / n**3 + + 337.931 * lm11(n, S1) + - 745.573 * lm11m1(n, S1) + + 8.61438 * lm13m1(n, S1, S2, S3) + - (385.52331999999996 / n) + + 690.151 / (1 + n) + - 656.386 / (2 + n) + + 656.386 / (3 + n) + ) + + # nf^2 (parametrized) and nf^3 (exact) + P3NSPA2 = -( + -193.85906555742952 + - 18.962964 / n**5 + + 99.1605 / n**4 + - 225.141 / n**3 + + 393.0056000000001 / n**2 + - 403.50217685814835 / n + - 34.425000000000004 / (1 + n) ** 4 + + 108.42 / (1 + n) ** 3 + - 93.8225 / (1 + n) ** 2 + + 534.725 / (1 + n) + + 246.50250000000003 / (2 + n) + - 25.455 / ((1 + n) ** 2 * (2 + n)) + - (16.97 * n) / ((1 + n) ** 2 * (2 + n)) + + 8.485 / ((1 + n) * (2 + n)) + - 110.015 / (3 + n) + + 78.9875 / (4 + n) + + 195.5772 * S1 + - (101.0775 * S1) / n**2 + + (35.17361 * S1) / n + - (8.485 * S1) / (1 + n) + - (101.0775 * S2) / n + ) + eta = 1 / n * 1 / (n + 1) + P3NSA3 = -CF * ( + -32 / 27 * zeta3 * eta + - 16 / 9 * zeta3 + - 16 / 27 * eta**4 + - 16 / 81 * eta**3 + + 80 / 27 * eta**2 + - 320 / 81 * eta + + 32 / 27 * 1 / (n + 1) ** 4 + + 128 / 27 * 1 / (n + 1) ** 2 + + 64 / 27 * S1 * zeta3 + - 32 / 81 * S1 + - 32 / 81 * S2 + - 160 / 81 * S3 + + 32 / 27 * S4 + + 131 / 81 + ) + + # Assembly regular piece. + P3NSPAI = P3NSA0 + nf * P3NSA1 + nf**2 * P3NSPA2 + nf**3 * P3NSA3 + if variation == 1: + P3NSPA = P3NSPAI + P3NPA01 + nf * P3NPA11 + elif variation == 2: + P3NSPA = P3NSPAI + P3NPA02 + nf * P3NPA12 + else: + P3NSPA = P3NSPAI + 0.5 * ((P3NPA01 + P3NPA02) + nf * (P3NPA11 + P3NPA12)) + + # The singular piece. + A4qI = ( + 2.120902 * 10**4 + - 5.179372 * 10**3 * nf + # + 1.955772 * 10**2 * nf**2 + # + 3.272344 * nf**3 + ) + A4ap1 = -507.152 + 7.33927 * nf + A4ap2 = -505.209 + 7.53662 * nf + D1 = 1 / n - S1 + if variation == 1: + P3NSPB = (A4qI + A4ap1) * D1 + elif variation == 2: + P3NSPB = (A4qI + A4ap2) * D1 + else: + P3NSPB = (A4qI + 0.5 * (A4ap1 + A4ap2)) * D1 + + # ..The local piece. + B4qI = ( + 2.579609 * 10**4 + + 0.08 + - (5.818637 * 10**3 + 0.97) * nf + # + (1.938554 * 10**2 + 0.0037) * nf**2 + # + 3.014982 * nf**3 + ) + B4ap1 = -2405.03 + 267.965 * nf + B4ap2 = -2394.47 + 269.028 * nf + if variation == 1: + P3NSPC = B4qI + B4ap1 + elif variation == 2: + P3NSPC = B4qI + B4ap2 + else: + P3NSPC = +B4qI + 0.5 * (B4ap1 + B4ap2) + + return -(P3NSPA + P3NSPB + P3NSPC) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsv.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsv.py new file mode 100644 index 000000000..9b7dcfecd --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gnsv.py @@ -0,0 +1,123 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{ns,v}^{(3)}`.""" +import numba as nb + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm11m1, lm12m1, lm13m1 +from .gnsm import gamma_nsm + + +@nb.njit(cache=True) +def gamma_nss(n, nf, cache, variation): + r"""Compute the |N3LO| sea non-single anomalous dimension. + + The routine is taken from :cite:`Moch:2017uml`. + + The :math:`nf^2` part is a high-accuracy (0.1% or better) parametrization + of the exact expression obtained in :cite:`Davies:2016jie`, see xpns3m.f + + The :math:`nf^1` part is an approximation based on the first 9 odd moments. + The two sets spanning the error estimate are called via IMOD = 1 + and IMOD = 2. Any other value of IMOD invokes their average. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache : numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| sea non-singlet anomalous dimension :math:`\gamma_{ns,s}^{(3)}` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + + # nf^1: two approximations + P3NSA11 = ( + 2880 / n**7 + - 11672.4 / n**6 + + 12802.560000000001 / n**5 + - 7626.66 / n**4 + + 6593.2 / n**3 + - 3687.6 / n**2 + + 4989.2 / (1 + n) + - 6596.93 / (2 + n) + + 1607.73 / (3 + n) + + 60.4 * lm12m1(n, S1, S2) + + 4.685 * lm13m1(n, S1, S2, S3) + ) + P3NSA12 = ( + -2880 / n**7 + + 4066.32 / n**6 + - 5682.24 / n**5 + + 5540.88 / n**4 + + 546.1 / n**3 + - 2987.83 / n**2 + + 2533.54 / n + - 1502.75 / (1 + n) + - 2297.56 / (2 + n) + + 1266.77 / (3 + n) + - 254.63 * lm11m1(n, S1) + - 0.28953 * lm13m1(n, S1, S2, S3) + ) + + # nf^2 (parametrized) + P3NSSA2 = ( + 47.4074 / n**6 + - 142.222 / n**5 + + 32.1201 / n**4 + - 132.824 / n**3 + + 647.397 / n**2 + + 19.7 * lm11m1(n, S1) + - 3.43547 * lm12m1(n, S1, S2) + - 1262.0951538579698 / n + - 187.17000000000002 / (1 + n) ** 4 + + 453.885 / (1 + n) ** 3 + + 147.01749999999998 / (1 + n) ** 2 + + 1614.1000000000001 / (1 + n) + - 380.12500000000006 / (2 + n) + - 42.575 / (3 + n) + + (42.977500000000006 * S2) / n + + (0.0900000000000047 * (477.52777777775293 + n) * S1) / (n**2 * (1 + n)) + ) + + if variation == 1: + P3NSSA = nf * P3NSA11 + nf**2 * P3NSSA2 + elif variation == 2: + P3NSSA = nf * P3NSA12 + nf**2 * P3NSSA2 + else: + P3NSSA = 0.5 * nf * (P3NSA11 + P3NSA12) + nf**2 * P3NSSA2 + return -P3NSSA + + +@nb.njit(cache=True) +def gamma_nsv(n, nf, cache, variation): + r"""Compute the |N3LO| valence non-singlet anomalous dimension. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| valence non-singlet anomalous dimension + :math:`\gamma_{ns,v}^{(3)}(N)` + + """ + return gamma_nsm(n, nf, cache, variation) + gamma_nss(n, nf, cache, variation) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gps.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gps.py new file mode 100644 index 000000000..87a8ba15d --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gps.py @@ -0,0 +1,165 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{ps}^{(3)}`.""" +import numba as nb + +from ......harmonics import cache as c +from ......harmonics.log_functions import ( + lm11m1, + lm12m1, + lm12m2, + lm13m1, + lm13m2, + lm14m1, + lm14m2, +) + + +@nb.njit(cache=True) +def gamma_ps(n, nf, cache, variation): + r"""Compute the |N3LO| pure singlet quark-quark anomalous dimension. + + The routine is taken from :cite:`Falcioni:2023luc`. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| pure singlet quark-quark anomalous dimension + :math:`\gamma_{ps}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + S4 = c.get(c.S4, cache, n) + nf2 = nf * nf + nf3 = nf * nf2 + xm1lm1 = -(1 / (-1 + n) ** 2) + 1 / n**2 + + # Known large-x coefficients + x1L4cff = -5.6460905 * 10 * nf + 3.6213992 * nf2 + x1L3cff = -2.4755054 * 10**2 * nf + 4.0559671 * 10 * nf2 - 1.5802469 * nf3 + y1L4cff = -1.3168724 * 10 * nf + y1L3cff = -1.9911111 * 10**2 * nf + 1.3695473 * 10 * nf2 + + # Known small-x coefficients + bfkl1 = 1.7492273 * 10**3 * nf + x0L6cff = -7.5061728 * nf + 7.9012346 * 0.1 * nf2 + x0L5cff = 2.8549794 * 10 * nf + 3.7925926 * nf2 + x0L4cff = -8.5480010 * 10**2 * nf + 7.7366255 * 10 * nf2 - 1.9753086 * 0.1 * nf3 + + # The resulting part of the function + P3ps01 = ( + +bfkl1 * 2 / (-1 + n) ** 3 + + x0L6cff * 720 / n**7 + + x0L5cff * -120 / n**6 + + x0L4cff * 24 / n**5 + + x1L3cff * lm13m1(n, S1, S2, S3) + + x1L4cff * lm14m1(n, S1, S2, S3, S4) + + y1L3cff * lm13m2(n, S1, S2, S3) + + y1L4cff * lm14m2(n, S1, S2, S3, S4) + ) + + # The selected approximations for nf = 3, 4, 5 + if nf == 3: + P3psApp1 = ( + P3ps01 + + 67731.0 * xm1lm1 + + 274100.0 * 1 / ((-1 + n) * n) + - 104493.0 * (1 / n - n / (2 + 3 * n + n**2)) + + 34403.0 * 1 / (6 + 5 * n + n**2) + + 353656.0 * (-(1 / n**2) + 1 / (1 + n) ** 2) + + 10620.0 * 2 / n**3 + + 40006.0 * -6 / n**4 + - 7412.1 * lm11m1(n, S1) + - 2365.1 * lm12m1(n, S1, S2) + + 1533.0 * lm12m2(n, S1, S2) + ) + P3psApp2 = ( + P3ps01 + + 54593.0 * xm1lm1 + + 179748.0 * 1 / ((-1 + n) * n) + - 195263.0 * 1 / (n + n**2) + + 12789.0 * 2 / (3 + 4 * n + n**2) + + 4700.0 * (-(1 / n**2) + 1 / (1 + n) ** 2) + - 103604.0 * 2 / n**3 + - 2758.3 * -6 / n**4 + - 2801.2 * lm11m1(n, S1) + - 1986.9 * lm12m1(n, S1, S2) + - 6005.9 * lm12m2(n, S1, S2) + ) + elif nf == 4: + P3psApp1 = ( + P3ps01 + + 90154.0 * xm1lm1 + + 359084.0 * 1 / ((-1 + n) * n) + - 136319.0 * (1 / n - n / (2 + 3 * n + n**2)) + + 45379.0 * 1 / (6 + 5 * n + n**2) + + 461167.0 * (-(1 / n**2) + 1 / (1 + n) ** 2) + + 13869.0 * 2 / n**3 + + 52525.0 * -6 / n**4 + - 7498.2 * lm11m1(n, S1) + - 2491.5 * lm12m1(n, S1, S2) + + 1727.2 * lm12m2(n, S1, S2) + ) + P3psApp2 = ( + P3ps01 + + 72987.0 * xm1lm1 + + 235802.0 * 1 / ((-1 + n) * n) + - 254921.0 * 1 / (n + n**2) + + 17138.0 * 2 / (3 + 4 * n + n**2) + + 5212.9 * (-(1 / n**2) + 1 / (1 + n) ** 2) + - 135378.0 * 2 / n**3 + - 3350.9 * -6 / n**4 + - 1472.7 * lm11m1(n, S1) + - 1997.2 * lm12m1(n, S1, S2) + - 8123.3 * lm12m2(n, S1, S2) + ) + elif nf == 5: + P3psApp1 = ( + P3ps01 + + 112481.0 * xm1lm1 + + 440555.0 * 1 / ((-1 + n) * n) + - 166581.0 * (1 / n - n / (2 + 3 * n + n**2)) + + 56087.0 * 1 / (6 + 5 * n + n**2) + + 562992.0 * (-(1 / n**2) + 1 / (1 + n) ** 2) + + 16882.0 * 2 / n**3 + + 64577.0 * -6 / n**4 + - 6570.1 * lm11m1(n, S1) + - 2365.7 * lm12m1(n, S1, S2) + + 1761.7 * lm12m2(n, S1, S2) + ) + P3psApp2 = ( + P3ps01 + + 91468.0 * xm1lm1 + + 289658.0 * 1 / ((-1 + n) * n) + - 311749.0 * 1 / (n + n**2) + + 21521.0 * 2 / (3 + 4 * n + n**2) + + 4908.9 * (-(1 / n**2) + 1 / (1 + n) ** 2) + - 165795.0 * 2 / n**3 + - 3814.9 * -6 / n**4 + + 804.5 * lm11m1(n, S1) + - 1760.8 * lm12m1(n, S1, S2) + - 10295.0 * lm12m2(n, S1, S2) + ) + else: + raise NotImplementedError("nf=6 is not available at N3LO") + + # We return (for now) one of the two error-band boundaries + # or the present best estimate, their average + if variation == 1: + P3psA = P3psApp1 + elif variation == 2: + P3psA = P3psApp2 + else: + P3psA = 0.5 * (P3psApp1 + P3psApp2) + return -P3psA diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gqg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gqg.py new file mode 100644 index 000000000..dd3244fd7 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/gqg.py @@ -0,0 +1,164 @@ +r"""The unpolarized, space-like anomalous dimension :math:`\gamma_{qg}^{(3)}`.""" +import numba as nb +import numpy as np + +from ......harmonics import cache as c +from ......harmonics.log_functions import lm11, lm12, lm13, lm14, lm14m1, lm15, lm15m1 + + +@nb.njit(cache=True) +def gamma_qg(n, nf, cache, variation): + r"""Compute the |N3LO| quark-gluon singlet anomalous dimension. + + The routine is taken from :cite:`Falcioni:2023vqq`. + + Parameters + ---------- + n : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : int + |N3LO| anomalous dimension variation + + Returns + ------- + complex + |N3LO| quark-gluon singlet anomalous dimension + :math:`\gamma_{qg}^{(3)}(N)` + + """ + S1 = c.get(c.S1, cache, n) + S2 = c.get(c.S2, cache, n) + S3 = c.get(c.S3, cache, n) + S4 = c.get(c.S4, cache, n) + S5 = c.get(c.S5, cache, n) + nf2 = nf * nf + nf3 = nf * nf2 + + # Known large-x coefficients + x1L5cff = 1.8518519 * nf - 4.1152263 * 0.1 * nf2 + x1L4cff = 3.5687794 * 10 * nf - 3.5116598 * nf2 - 8.2304527 * 0.01 * nf3 + y1L5cff = 2.8806584 * nf + 8.2304527 * 0.1 * nf2 + y1L4cff = -4.0511391 * 10 * nf + 5.5418381 * nf2 + 1.6460905 * 0.1 * nf3 + + # Known small-x coefficients + bfkl1 = 3.9357613 * 10**3 * nf + x0L6cff = -1.9588477 * 10 * nf + 2.7654321 * nf2 + x0L5cff = 2.1573663 * 10 * nf + 1.7244444 * 10 * nf2 + x0L4cff = -2.8667643 * 10**3 * nf + 3.0122403 * 10**2 * nf2 + 4.1316872 * nf3 + + # The resulting part of the function + P3QG01 = ( + +bfkl1 * 2 / (-1 + n) ** 3 + + x0L6cff * 720 / n**7 + + x0L5cff * -120 / n**6 + + x0L4cff * 24 / n**5 + + x1L4cff * lm14(n, S1, S2, S3, S4) + + x1L5cff * lm15(n, S1, S2, S3, S4, S5) + + y1L4cff * lm14m1(n, S1, S2, S3, S4) + + y1L5cff * lm15m1(n, S1, S2, S3, S4, S5) + ) + + # The selected approximations for nf = 3, 4, 5 + if nf == 3: + P3qgApp1 = ( + P3QG01 + + 187500.0 * -(1 / (-1 + n) ** 2) + + 826060.0 * 1 / ((-1 + n) * n) + - 150474.0 * 1 / n + + 226254.0 * (3 + n) / (2 + 3 * n + n**2) + + 577733.0 * -1 / n**2 + - 180747.0 * 2 / n**3 + + 95411.0 * -6 / n**4 + + 119.8 * lm13(n, S1, S2, S3) + + 7156.3 * lm12(n, S1, S2) + + 45790.0 * lm11(n, S1) + - 95682.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + P3qgApp2 = ( + P3QG01 + + 135000.0 * -(1 / (-1 + n) ** 2) + + 484742.0 * 1 / ((-1 + n) * n) + - 11627.0 * 1 / n + - 187478.0 * (3 + n) / (2 + 3 * n + n**2) + + 413512.0 * -1 / n**2 + - 82500.0 * 2 / n**3 + + 29987.0 * -6 / n**4 + - 850.1 * lm13(n, S1, S2, S3) + - 11425.0 * lm12(n, S1, S2) + - 75323.0 * lm11(n, S1) + + 282836.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + elif nf == 4: + P3qgApp1 = ( + P3QG01 + + 250000.0 * -(1 / (-1 + n) ** 2) + + 1089180.0 * 1 / ((-1 + n) * n) + - 241088.0 * 1 / n + + 342902.0 * (3 + n) / (2 + 3 * n + n**2) + + 720081.0 * -1 / n**2 + - 247071.0 * 2 / n**3 + + 126405.0 * -6 / n**4 + + 272.4 * lm13(n, S1, S2, S3) + + 10911.0 * lm12(n, S1, S2) + + 60563.0 * lm11(n, S1) + - 161448.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + P3qgApp2 = ( + P3QG01 + + 180000.0 * -(1 / (-1 + n) ** 2) + + 634090.0 * 1 / ((-1 + n) * n) + - 55958.0 * 1 / n + - 208744.0 * (3 + n) / (2 + 3 * n + n**2) + + 501120.0 * -1 / n**2 + - 116073.0 * 2 / n**3 + + 39173.0 * -6 / n**4 + - 1020.8 * lm13(n, S1, S2, S3) + - 13864.0 * lm12(n, S1, S2) + - 100922.0 * lm11(n, S1) + + 343243.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + elif nf == 5: + P3qgApp1 = ( + P3QG01 + + 312500.0 * -(1 / (-1 + n) ** 2) + + 1345700.0 * 1 / ((-1 + n) * n) + - 350466.0 * 1 / n + + 480028.0 * (3 + n) / (2 + 3 * n + n**2) + + 837903.0 * -1 / n**2 + - 315928.0 * 2 / n**3 + + 157086.0 * -6 / n**4 + + 472.7 * lm13(n, S1, S2, S3) + + 15415.0 * lm12(n, S1, S2) + + 75644.0 * lm11(n, S1) + - 244869.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + P3qgApp2 = ( + P3QG01 + + 225000.0 * -(1 / (-1 + n) ** 2) + + 776837.0 * 1 / ((-1 + n) * n) + - 119054.0 * 1 / n + - 209530.0 * (3 + n) / (2 + 3 * n + n**2) + + 564202.0 * -1 / n**2 + - 152181.0 * 2 / n**3 + + 48046.0 * -6 / n**4 + - 1143.8 * lm13(n, S1, S2, S3) + - 15553.0 * lm12(n, S1, S2) + - 126212.0 * lm11(n, S1) + + 385995.0 * (S1 - n * (np.pi**2 / 6 - S2)) / n**2 + ) + else: + raise NotImplementedError("nf=6 is not available at N3LO") + + # We return one of the two error-band representatives + # or the present best estimate, their average + if variation == 1: + P3QGA = P3qgApp1 + elif variation == 2: + P3QGA = P3qgApp2 + else: + P3QGA = 0.5 * (P3qgApp1 + P3qgApp2) + return -P3QGA diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index 47d5e7001..fa8eb77a0 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -42,7 +42,8 @@ def test_quad_ker_errors(): is_threshold=False, is_polarized=True, is_time_like=True, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + use_fhmruvv=False, ) @@ -94,7 +95,8 @@ def test_quad_ker(monkeypatch): is_threshold=False, is_polarized=p, is_time_like=t, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + use_fhmruvv=False, ) np.testing.assert_allclose(res_ns, res) for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: @@ -122,7 +124,8 @@ def test_quad_ker(monkeypatch): is_threshold=False, is_polarized=polarized, is_time_like=False, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + use_fhmruvv=False, ) np.testing.assert_allclose(res_sv, 1.0) for label in [ @@ -155,9 +158,10 @@ def test_quad_ker(monkeypatch): ev_op_max_order=(1, 0), sv_mode=sv, is_threshold=False, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), is_polarized=False, is_time_like=False, + use_fhmruvv=False, ) np.testing.assert_allclose(res_sv, 1.0) @@ -182,9 +186,10 @@ def test_quad_ker(monkeypatch): ev_op_max_order=(0, 0), sv_mode=1, is_threshold=False, - n3lo_ad_variation=(0, 0, 0, 0), + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), is_polarized=False, is_time_like=False, + use_fhmruvv=False, ) np.testing.assert_allclose(res_ns, 0.0) @@ -403,7 +408,9 @@ def quad_ker_pegasus( phi = 3 / 4 * np.pi c = 1.9 n = complex(c + u * np.exp(1j * phi)) - gamma_ns = ad.gamma_ns(order, mode0, n, nf) + gamma_ns = ad.gamma_ns( + order, mode0, n, nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0) + ) ker = ns.dispatcher( order, method, @@ -460,6 +467,7 @@ def quad_ker_pegasus( (0, 0, 0, 0), False, False, + False, ), epsabs=1e-12, epsrel=1e-5, diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index e6179ba99..bb18c3e17 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -77,7 +77,7 @@ def test_zero_true_gamma(): for qed in range(1, 2 + 1): order = (qcd, qed) n = np.random.rand() - gamma_ns = ad.gamma_ns_qed(order, mode, n, nf) + gamma_ns = ad.gamma_ns_qed(order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0)) for method in methods: for running in running_alpha: np.testing.assert_allclose( @@ -248,7 +248,9 @@ def test_ode_true_gamma(): ) a1 = sc.a_s(mu2_to) n = 3 + np.random.rand() - gamma_ns = ad.gamma_ns_qed(order, mode, n, nf) + gamma_ns = ad.gamma_ns_qed( + order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0) + ) gammatot = 0.0 for i in range(0, order[0] + 1): for j in range(0, order[1] + 1): diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py index a27bc6ef6..29f0f3d7a 100644 --- a/tests/eko/kernels/test_kernels_QEDvalence.py +++ b/tests/eko/kernels/test_kernels_QEDvalence.py @@ -67,7 +67,7 @@ def test_zero_true_gamma(monkeypatch): order = (qcd, qed) n = np.random.rand() gamma_v = anomalous_dimensions.unpolarized.space_like.gamma_valence_qed( - order, n, nf + order, n, nf, (0, 0, 0, 0, 0, 0, 0) ) for method in methods: np.testing.assert_allclose( diff --git a/tests/eko/kernels/test_ns.py b/tests/eko/kernels/test_ns.py index 013fc5433..86abef4bf 100644 --- a/tests/eko/kernels/test_ns.py +++ b/tests/eko/kernels/test_ns.py @@ -207,7 +207,7 @@ def test_error(monkeypatch): 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) with pytest.raises(NotImplementedError): - ad.gamma_ns((2, 0), 10202, 1, 3) + ad.gamma_ns((2, 0), 10202, 1, (0, 0, 0, 0, 0, 0, 0), 3) def test_gamma_usage(): diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 7ca91333d..d85b8926a 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -159,7 +159,13 @@ def scheme_diff(g, k, pto, is_singlet): # for order in [(2, 0), (3, 0), (4, 0)]: for order in [(2, 0)]: # Non singlet kernels - gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf) + gns = gamma_ns( + order, + br.non_singlet_pids_map["ns+"], + n, + nf, + n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), + ) ker = non_singlet.dispatcher( order, method, gns, a1, a0, nf, ev_op_iterations=1 ) @@ -177,7 +183,7 @@ def scheme_diff(g, k, pto, is_singlet): ) # Singlet kernels - gs = gamma_singlet(order, n, nf, n3lo_ad_variation=(0, 0, 0, 0)) + gs = gamma_singlet(order, n, nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0)) ker = singlet.dispatcher( order, method, gs, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1 ) diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py new file mode 100644 index 000000000..8bef494c7 --- /dev/null +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py @@ -0,0 +1,272 @@ +# Test N3LO anomalous dimensions +import numpy as np +import pytest + +from ekore import harmonics as h +from ekore.anomalous_dimensions.unpolarized.space_like.as4.fhmruvv import ( + gamma_singlet, + ggg, + ggq, + gnsm, + gnsp, + gnsv, + gps, + gqg, +) + +NF = 5 + + +def test_quark_number_conservation(): + N = 1 + sx_cache = h.cache.reset() + for imod in [0, 1, 2]: + # (ns,s) + refs = [-0.01100459, -0.007799, -0.01421] + np.testing.assert_allclose( + gnsv.gamma_nss(N, NF, sx_cache, imod), refs[imod], rtol=6e-05 + ) + + # (ns,m) + refs = [0.06776363, 0.064837, 0.07069] + np.testing.assert_allclose( + gnsm.gamma_nsm(N, NF, sx_cache, imod), refs[imod], rtol=6e-06 + ) + + +def test_momentum_conservation(): + N = 2 + sx_cache = h.cache.reset() + + g_singlet = np.zeros((3, 2, 2), dtype=complex) + for imod in [(0, 0, 0, 0), (1, 1, 1, 1), (2, 2, 2, 2)]: + g_singlet[imod, :, :] = gamma_singlet(N, NF, sx_cache, imod) + + # total + np.testing.assert_allclose( + g_singlet[:, 0, 0] + g_singlet[:, 1, 0], + [0.08617, 0.220242, -0.047901], + atol=2e-5, + ) + np.testing.assert_allclose( + g_singlet[:, 0, 1] + g_singlet[:, 1, 1], + [-0.134766, 0.465174, -0.734706], + atol=2e-5, + ) + + +def test_vogt_parametriztions(): + def qg3_moment(N, nf): + mom_list = [ + -654.4627782205557 * nf + + 245.61061978871788 * nf**2 + - 0.9249909688301847 * nf**3, + 290.31106867034487 * nf + - 76.51672403736478 * nf**2 + - 4.911625629947491 * nf**3, + 335.80080466045274 * nf + - 124.57102255718002 * nf**2 + - 4.193871425027802 * nf**3, + 294.58768309440677 * nf + - 135.3767647714609 * nf**2 + - 3.609775642729055 * nf**3, + 241.6153399044715 * nf + - 135.18742470907011 * nf**2 + - 3.189394834180898 * nf**3, + 191.97124640777176 * nf + - 131.16316638326697 * nf**2 + - 2.8771044305171913 * nf**3, + 148.5682948286098 * nf + - 125.82310814280595 * nf**2 + - 2.635918561148907 * nf**3, + 111.34042526856348 * nf + - 120.16819876888667 * nf**2 + - 2.4433790398202664 * nf**3, + 79.51561588665083 * nf + - 114.61713540075442 * nf**2 + - 2.28548686108789 * nf**3, + 52.24329555231736 * nf + - 109.34248910828198 * nf**2 + - 2.1531537251387527 * nf**3, + ] + return mom_list[int((N - 2) / 2)] + + def qq3ps_moment(N, nf): + mom_list = [ + -691.5937093082381 * nf + + 84.77398149891167 * nf**2 + + 4.4669568492355864 * nf**3, + -109.33023358432462 * nf + + 8.77688525974872 * nf**2 + + 0.3060771365698822 * nf**3, + -46.030613749542226 * nf + + 4.744075766957513 * nf**2 + + 0.042548957282380874 * nf**3, + -24.01455020567638 * nf + + 3.235193483272451 * nf**2 + - 0.007889256298951614 * nf**3, + -13.730393879922417 * nf + + 2.3750187592472374 * nf**2 + - 0.02102924056123573 * nf**3, + -8.152592251923657 * nf + + 1.8199581788320662 * nf**2 + - 0.024330231290833188 * nf**3, + -4.8404471801109565 * nf + + 1.4383273806219803 * nf**2 + - 0.024479943136069916 * nf**3, + -2.7511363301137024 * nf + + 1.164299642517469 * nf**2 + - 0.023546009234463816 * nf**3, + -1.375969240387974 * nf + + 0.9608733183576097 * nf**2 + - 0.022264393374041958 * nf**3, + -0.4426815682220422 * nf + + 0.8057453328332964 * nf**2 + - 0.02091826436475512 * nf**3, + ] + return mom_list[int((N - 2) / 2)] + + def gg3_moment(N, nf): + mom_list = [ + 654.4627782205557 * nf + - 245.6106197887179 * nf**2 + + 0.9249909688301847 * nf**3, + 39876.123276008046 + - 10103.4511350227 * nf + + 437.0988475397789 * nf**2 + + 12.955565459350593 * nf**3, + 53563.84353419538 + - 14339.131035160317 * nf + + 652.7773306808972 * nf**2 + + 16.654103652963503 * nf**3, + 62279.7437813437 + - 17150.696783851945 * nf + + 785.8806126875509 * nf**2 + + 18.933103109772713 * nf**3, + ] + return mom_list[int((N - 2) / 2)] + + def gq3_moment(N, nf): + mom_list = [ + -16663.225488 + + 4439.143749608238 * nf + - 202.55547919891168 * nf**2 + - 6.375390720235586 * nf**3, + -6565.7531450230645 + + 1291.0067460871576 * nf + - 16.146190170051486 * nf**2 + - 0.8397634037808341 * nf**3, + -3937.479370556893 + + 679.7185057363981 * nf + - 1.3720775271604673 * nf**2 + - 0.13979432728276966 * nf**3, + -2803.644107251366 + + 436.39305738710254 * nf + + 1.8149462465491055 * nf**2 + + 0.07358858022119033 * nf**3, + ] + return mom_list[int((N - 2) / 2)] + + for N in [2, 4, 6, 8, 10, 12, 14, 16, 18, 20]: + cache = h.cache.reset() + for variation in [0, 1, 2]: + for nf in [3, 4, 5]: + np.testing.assert_allclose( + gps.gamma_ps(N, nf, cache, variation), + qq3ps_moment(N, nf), + rtol=4e-4, + ) + np.testing.assert_allclose( + gqg.gamma_qg(N, nf, cache, variation), + qg3_moment(N, nf), + rtol=8e-4, + ) + if N <= 8: + np.testing.assert_allclose( + ggg.gamma_gg(N, nf, cache, variation), + gg3_moment(N, nf), + rtol=4e-4, + ) + np.testing.assert_allclose( + ggq.gamma_gq(N, nf, cache, variation), + gq3_moment(N, nf), + rtol=2e-4, + ) + + with pytest.raises(NotImplementedError): + gps.gamma_ps(N, 6, cache, variation) + gqg.gamma_qg(N, 6, cache, variation) + ggg.gamma_gg(N, 6, cache, variation) + ggq.gamma_gq(N, 6, cache, variation) + + +def test_non_singlet_reference_moments(): + NF = 4 + nsm_nf4_refs = [ + 4322.890485339998, + 5491.581109692005, + 6221.256799360004, + 6774.606221595994, + 7229.056043916002, + 7618.358743427995, + 7960.658678124, + ] + nss_nf4_refs = [ + 50.10532524, + 39.001939964, + 21.141505811200002, + 12.4834195012, + 8.0006134908, + 5.4610639744, + 3.9114290952, + 2.90857799, + ] + for N in [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0]: + sx_cache = h.cache.reset() + idx = int((N - 3) / 2) + if N != 17: + np.testing.assert_allclose( + gnsm.gamma_nsm(N, NF, sx_cache, variation=0), + nsm_nf4_refs[idx], + rtol=8e-5, + ) + np.testing.assert_allclose( + gnsv.gamma_nsv(N, NF, sx_cache, variation=0), + nss_nf4_refs[idx] + nsm_nf4_refs[idx], + rtol=8e-5, + ) + np.testing.assert_allclose( + gnsv.gamma_nss(N, NF, sx_cache, variation=0), nss_nf4_refs[idx], rtol=5e-5 + ) + + +def test_singlet_reference_moments(): + NF = 4 + nsp_nf4_refs = [ + 3679.6690577439995, + 5066.339235808004, + 5908.005605364002, + 6522.700744595994, + 7016.383458928004, + 7433.340927783997, + 7796.397038483998, + 8119.044600816003, + ] + for N in [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0]: + sx_cache = h.cache.reset() + np.testing.assert_allclose( + gnsp.gamma_nsp(N, NF, sx_cache, variation=0), + nsp_nf4_refs[int((N - 2) / 2)], + rtol=4e-5, + ) + + +def test_gamma_ps_extrapolation(): + # Test the prediction of N=22 wrt to :cite:`Falcioni:2023luc` + n22_ref = [6.2478570, 10.5202730, 15.6913948] + N = 22 + sx_cache = h.cache.reset() + my_res = [] + for nf in [3, 4, 5]: + my_res.append(gps.gamma_ps(N, nf, sx_cache, 0)) + np.testing.assert_allclose(n22_ref, n22_ref) diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py index c7c955960..5fdc1ee06 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py @@ -70,7 +70,7 @@ def test_eigensystem_gamma_singlet_projectors_EV(): # ignore Runtime Warnings warnings.simplefilter("ignore", RuntimeWarning) for gamma_S in ad_us.gamma_singlet( - o, N, nf, n3lo_ad_variation=(0, 0, 0, 0) + o, N, nf, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0) ): _exp, l_p, l_m, e_p, e_m = ad.exp_matrix_2D(gamma_S) # projectors behave as P_a . P_b = delta_ab P_a @@ -84,113 +84,153 @@ def test_eigensystem_gamma_singlet_projectors_EV(): def test_gamma_ns(): nf = 3 + n3lo_ad_variation = (0, 0, 0, 0, 0, 0, 0) # ad_us.as1 assert_almost_equal( - ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0 + ad_us.gamma_ns( + (3, 0), br.non_singlet_pids_map["ns+"], 1, nf, n3lo_ad_variation + )[0], + 0.0, ) # ad_us.as2 assert_allclose( - ad_us.gamma_ns((2, 0), br.non_singlet_pids_map["ns-"], 1, nf), + ad_us.gamma_ns( + (2, 0), br.non_singlet_pids_map["ns-"], 1, nf, n3lo_ad_variation + ), np.zeros(2), atol=2e-6, ) # ad_us.as3 assert_allclose( - ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["ns-"], 1, nf), + ad_us.gamma_ns( + (3, 0), br.non_singlet_pids_map["ns-"], 1, nf, n3lo_ad_variation + ), np.zeros(3), atol=2e-4, ) assert_allclose( - ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["nsV"], 1, nf), + ad_us.gamma_ns( + (3, 0), br.non_singlet_pids_map["nsV"], 1, nf, n3lo_ad_variation + ), np.zeros(3), atol=8e-4, ) # as4 assert_allclose( - ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["ns-"], 1, nf), + ad_us.gamma_ns( + (4, 0), br.non_singlet_pids_map["ns-"], 1, nf, n3lo_ad_variation + ), np.zeros(4), atol=2e-4, ) # N3LO valence has a spurious pole, need to add a small shift assert_allclose( - ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf), + ad_us.gamma_ns( + (4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf, n3lo_ad_variation + ), np.zeros(4), atol=5e-4, ) assert_raises( AssertionError, assert_allclose, - ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["ns+"], 1, nf), + ad_us.gamma_ns( + (4, 0), br.non_singlet_pids_map["ns+"], 1, nf, n3lo_ad_variation + ), np.zeros(4), ) with pytest.raises(NotImplementedError): - ad_us.gamma_ns((2, 0), 10106, 2.0, nf) + ad_us.gamma_ns((2, 0), 10106, 2.0, nf, n3lo_ad_variation) def test_gamma_ns_qed(): nf = 3 # aem1 assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns+u"], 1, nf)[0, 1], + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns+u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + )[0, 1], 0, decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns+d"], 1, nf)[0, 1], + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns+d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + )[0, 1], 0, decimal=5, ) # as1aem1 assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) # aem2 assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) # ad_us.as2 assert_almost_equal( - ad_us.gamma_ns_qed((2, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (2, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((3, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((2, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (2, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((3, 2)), decimal=5, ) # ad_us.as3 assert_almost_equal( - ad_us.gamma_ns_qed((3, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (3, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((4, 2)), decimal=3, ) assert_almost_equal( - ad_us.gamma_ns_qed((3, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (3, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((4, 2)), decimal=3, ) @@ -222,7 +262,7 @@ def test_dim_valence(): nf = 3 N = 2 cache = h.cache.reset() - gamma_valence = ad_us.gamma_valence_qed((3, 2), N, nf) + gamma_valence = ad_us.gamma_valence_qed((3, 2), N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_valence.shape == (4, 3, 2, 2) gamma_valence_as1 = ad_us.as1.gamma_valence_qed(N, cache) assert gamma_valence_as1.shape == (2, 2) @@ -235,9 +275,9 @@ def test_dim_valence(): def test_dim_nsp(): nf = 3 N = 2 - gamma_nsup = ad_us.gamma_ns_qed((3, 2), 10102, N, nf) + gamma_nsup = ad_us.gamma_ns_qed((3, 2), 10102, N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_nsup.shape == (4, 3) - gamma_nsdp = ad_us.gamma_ns_qed((3, 2), 10103, N, nf) + gamma_nsdp = ad_us.gamma_ns_qed((3, 2), 10103, N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_nsdp.shape == (4, 3) with pytest.raises(NotImplementedError): - ad_us.gamma_ns_qed((2, 0), 10106, N, nf) + ad_us.gamma_ns_qed((2, 0), 10106, N, nf, (0, 0, 0, 0, 0, 0, 0))