N3LO benchmark study

extras/lh_bench_23/.gitignore

@@ -1,3 +1,4 @@
 *.tar
 *.csv
 *.dat
+*.tex

extras/lh_bench_23/cfg.py

@@ -62,10 +62,12 @@ def vfns_theory(xif=1.0):
 ]
 
 
-def ffns_theory(xif=1.0):
+def ffns_theory(xif=1.0, pto=2):
     """Generate a VFNS theory card."""
     tt = copy.deepcopy(_t_ffns)
     tt["xif"] = xif
+    tt["order"] = (pto + 1, 0)
+    tt["matching_order"] = (pto, 0)
     return runcards.TheoryCard.from_dict(tt)
 
 
@@ -107,9 +109,15 @@ def n3lo_theory(ad_variation, is_ffns, use_fhmruvv=False, xif=1.0):
 )
 vfns_operator = runcards.OperatorCard.from_dict(_o_vfns)
 
-_o_ffns = copy.deepcopy(_o_vfns)
-_o_ffns["mugrid"] = [(100.0, 4)]
-ffns_operator = runcards.OperatorCard.from_dict(_o_ffns)
+
+def ffns_operator(ev_method="iterate-exact"):
+    """Generate a FFNS theory card."""
+    op = copy.deepcopy(_o_vfns)
+    op["mugrid"] = [(100.0, 4)]
+    op["configs"]["evolution_method"] = ev_method
+    if ev_method == "truncated":
+        op["configs"]["ev_op_iterations"] = 1
+    return runcards.OperatorCard.from_dict(op)
 
 
 # flavor rotations

extras/lh_bench_23/parse_to_latex.py

@@ -0,0 +1,179 @@
+from cfg import here, table_dir, xgrid
+from utils import compute_n3lo_avg_err, load_n3lo_tables
+
+n3lo_table_dir = table_dir
+
+latex_tab = here / "latex_tab"
+latex_tab.mkdir(exist_ok=True)
+
+SVS = ["central", "up", "down"]
+
+MIDRULE1 = r"""
+\hline \hline
+\multicolumn{9}{||c||}{} \\[-3mm]
+\multicolumn{9}{||c||}{"""
+
+MIDRULE2 = r"""} \\
+\multicolumn{9}{||c||}{} \\[-0.3cm]
+\hline \hline
+ & & & & & & & \\[-0.3cm]
+"""
+
+BOTTOMRULE = r"""
+\hline \hline
+\end{tabular}
+\end{center}
+\end{table}
+"""
+
+VFNS_LABELS = r"""
+    \multicolumn{1}{c|} {$xu_v$} &
+    \multicolumn{1}{c|} {$xd_v$} &
+    \multicolumn{1}{c|} {$xL_-$} &
+    \multicolumn{1}{c|} {$xL_+$} &
+    \multicolumn{1}{c|} {$xs_+$} &
+    \multicolumn{1}{c|} {$xc_+$} &
+    \multicolumn{1}{c|} {$xb_+$} &
+    \multicolumn{1}{c||}{$xg$} \\[0.5mm]
+    """
+
+FFNS_LABELS = r"""
+    \multicolumn{1}{c|} {$xu_v$} &
+    \multicolumn{1}{c|} {$xd_v$} &
+    \multicolumn{1}{c|} {$xL_-$} &
+    \multicolumn{1}{c|} {$xL_+$} &
+    \multicolumn{1}{c|} {$xs_v$} &
+    \multicolumn{1}{c|} {$xs_+$} &
+    \multicolumn{1}{c|} {$xc_+$} &
+    \multicolumn{1}{c||}{$xg$}
+    """
+
+
+def insert_headrule(scheme, approx, caption):
+    """Insert the middle rule."""
+    label = r"\label{tab:" + f"n3lo_{scheme.lower()}_{approx.lower()}" + "}"
+    scheme_label = (
+        r", $\, N_{\rm f} = 3\ldots 5\,$,"
+        if scheme == "VFNS"
+        else r"$\, N_{\rm f} = 4$,"
+    )
+    HEADRULE = (
+        r"""
+    \begin{table}[htp]
+    \caption{"""
+        + caption
+        + r"""}
+    """
+        + label
+        + r"""
+    \begin{center}
+    \vspace{5mm}
+    \begin{tabular}{||c||r|r|r|r|r|r|r|r||}
+    \hline \hline
+    \multicolumn{9}{||c||}{} \\[-3mm]
+    \multicolumn{9}{||c||}{"""
+        # + r"""aN$^3$LO, """
+        + approx
+        + scheme_label
+        + r"""$\,\mu_{\rm f}^2 = 10^4 \mbox{ GeV}^2$} \\
+    \multicolumn{9}{||c||}{} \\[-0.3cm]
+    \hline \hline
+    \multicolumn{9}{||c||}{} \\[-3mm]
+    \multicolumn{1}{||c||}{$x$} &
+    """
+    )
+    HEADRULE += VFNS_LABELS if scheme == "VFNS" else FFNS_LABELS
+    HEADRULE += r"""\\[0.5mm]"""
+    return HEADRULE
+
+
+def insert_midrule(sv):
+    """Insert the middle rule."""
+    # TODO: is this mapping correct or the other way round ??
+    # xif2 = 2 -> up
+    # xif2 = 1/2 -> down
+    label = {
+        "central": r"$\mu_{\rm r}^2 = \ \mu_{\rm f}^2$",
+        "down": r"$\mu_{\rm r}^2 = 0.5 \ \mu_{\rm f}^2$",
+        "up": r"$\mu_{\rm r}^2 = 2 \ \mu_{\rm f}^2$",
+    }
+    return MIDRULE1 + label[sv] + MIDRULE2
+
+
+def format_float(values):
+    """Clean float format."""
+    values = values.replace("0000", "0")
+    values = values.replace("e-0", "$^{-")
+    values = values.replace("e-10", "$^{-10")
+    values = values.replace("e+0", "$^{+")
+    values = values.replace("&", "}$ &")
+    values = values.replace(r"\\", r"}$ \\")
+    return values
+
+
+def dump_table(scheme: str, approx: str, caption: str):
+    """Write a nice latex table."""
+    final_tab = insert_headrule(scheme, approx.replace("EKO", "NNPDF"), caption)
+    # loop on scales
+    for sv in SVS:
+        # load tables
+        dfs = load_n3lo_tables(n3lo_table_dir, scheme, sv=sv, approx=approx)
+
+        central, err = compute_n3lo_avg_err(dfs)
+
+        central.insert(0, "x", xgrid)
+        values = central.to_latex(float_format="{:.4e}".format, index=False)
+        values = "".join(e for e in values.split("\n")[4:-3])
+        final_tab += insert_midrule(sv) + format_float(values)
+
+    final_tab += BOTTOMRULE
+
+    # write
+    with open(
+        latex_tab / f"table-{scheme}-{approx.replace('EKO', 'NNPDF')}.tex",
+        "w",
+        encoding="utf-8",
+    ) as f:
+        f.writelines(final_tab)
+
+
+if __name__ == "__main__":
+    approx = "FHMRUVV"
+    scheme = "FFNS"
+    caption = r"""
+        Results for the FFNS aN$^3$LO evolution
+        for the initial conditions and the input parton distributions
+        given in Sec.~\ref{sec:toy_pdf},
+        with the FHMRUVV splitting functions approximation and the NNPDF code.
+    """
+    dump_table(scheme, approx, caption)
+
+    approx = "FHMRUVV"
+    scheme = "VFNS"
+    caption = r"""
+        Results for the VFNS aN$^3$LO evolution
+        for the initial conditions and the input parton distributions
+        given in Sec.~\ref{sec:toy_pdf},
+        with the FHMRUVV splitting functions approximation and the NNPDF code.
+    """
+    dump_table(scheme, approx, caption)
+
+    approx = "EKO"
+    scheme = "FFNS"
+    caption = r"""
+        Results for the FFNS aN$^3$LO evolution
+        for the initial conditions and the input parton distributions
+        given in Sec.~\ref{sec:toy_pdf},
+        with the NNPDF splitting functions approximation.
+    """
+    dump_table(scheme, approx, caption)
+
+    approx = "EKO"
+    scheme = "VFNS"
+    caption = r"""
+        Results for the VFNS aN$^3$LO evolution
+        for the initial conditions and the input parton distributions
+        given in Sec.~\ref{sec:toy_pdf},
+        with the NNPDF splitting functions approximation.
+    """
+    dump_table(scheme, approx, caption)

extras/lh_bench_23/plot_bench_evol.py

@@ -18,28 +18,32 @@
 
 
 # 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)
+eko_dfs = load_n3lo_tables(
+    n3lo_table_dir, SCHEME, sv="central", approx="EKO", rotate_to_evol=True
+)
+fhmruvv_dfs = load_n3lo_tables(
+    n3lo_table_dir, SCHEME, sv="central", approx="FHMRUVV", 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)
+fhmruvv_res = compute_n3lo_avg_err(fhmruvv_dfs)
 
 n3lo_dfs = [
     (eko_res, "aN3LO EKO"),
-    (fhmv_res, "aN3LO FHMV"),
+    (fhmruvv_res, "aN3LO FHMRUVV"),
 ]
 
 # 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)
+fhmruvv_diff = compute_n3lo_nnlo_diff(fhmruvv_res, nnlo_central, REL_DIFF)
 n3lo_dfs = [
     (eko_diff, "aN3LO EKO"),
-    (fhmv_diff, "aN3LO FHMV"),
+    (fhmruvv_diff, "aN3LO FHMRUVV"),
 ]
 
 plot_diff_to_nnlo(xgrid, n3lo_dfs, SCHEME, USE_LINX, plot_dir, REL_DIFF)

extras/lh_bench_23/plot_bench_msht.py

@@ -15,48 +15,51 @@
 SV = "central"
 
 plot_dir = here / "plots_msht"
-n3lo_table_dir = table_dir  # / SCHEME
+n3lo_table_dir = table_dir
 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")
+eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, SV, approx="EKO")
+
+fhmruvv_eko_dfs = load_n3lo_tables(n3lo_table_dir, SCHEME, SV, approx="FHMRUVV")
+fhmruvv_msht_dfs = load_msht(msht_table_dir, SCHEME, approx="FHMRUVV")
+
+msht_post_dfs = load_msht(msht_table_dir, SCHEME, approx="MSHTposterior")
+msht_prior_dfs = load_msht(msht_table_dir, SCHEME, approx="MSHTprior")
 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)
+fhmruvv_eko_res = compute_n3lo_avg_err(fhmruvv_eko_dfs)
+fhmruvv_msht_res = compute_n3lo_avg_err(fhmruvv_msht_dfs)
+msht_post_res = compute_n3lo_avg_err(msht_post_dfs)
+msht_prior_res = compute_n3lo_avg_err(msht_prior_dfs)
 
-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"),
-]
+# compute average of FHMRUVV
+fhmruvv_res = []
+for a, b in zip(fhmruvv_msht_res, fhmruvv_eko_res):
+    fhmruvv_res.append((a + b) / 2)
 
 # PDFs plots
+n3lo_dfs = [
+    (fhmruvv_res, "FHMRUVV"),
+    (msht_prior_res, "MSHT (prior)"),
+    (msht_post_res, "MSHT (posterior)"),
+    (eko_res, "NNPDF"),
+]
 plot_pdfs(xgrid, n3lo_dfs, nnlo_central, SCHEME, USE_LINX, plot_dir)
 
-# relative diff plots
+# relative, absolute 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)
+fhmruvv_diff = compute_n3lo_nnlo_diff(fhmruvv_res, nnlo_central, REL_DIFF)
+msht_prior_diff = compute_n3lo_nnlo_diff(msht_prior_res, nnlo_central, REL_DIFF)
+msht_post_diff = compute_n3lo_nnlo_diff(msht_post_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"),
+    (fhmruvv_diff, "FHMRUVV"),
+    (msht_prior_diff, "MSHT (prior)"),
+    (msht_post_diff, "MSHT (posterior)"),
+    (eko_diff, "NNPDF"),
 ]
-
-# relative, absolute diff plots
 plot_diff_to_nnlo(xgrid, n3lo_dfs, SCHEME, USE_LINX, plot_dir, REL_DIFF)