From 31be4f533e5909c4d7de9d4bbd819c7de892a425 Mon Sep 17 00:00:00 2001
From: enocera <emanueleroberto.nocera@unito.it>
Date: Thu, 14 Jul 2022 16:58:21 +0100
Subject: [PATCH] Implemented function to compute asymptotic exponents
---
.../asy_exponents/asy_exponents_plots.yaml | 42 +++
.../asy_exponents/asy_exponents_table.yaml | 28 ++
validphys2/src/validphys/app.py | 1 +
validphys2/src/validphys/asy_exponents.py | 266 ++++++++++++++++++
validphys2/src/validphys/pdfbases.py | 22 ++
5 files changed, 359 insertions(+)
create mode 100644 validphys2/examples/asy_exponents/asy_exponents_plots.yaml
create mode 100644 validphys2/examples/asy_exponents/asy_exponents_table.yaml
create mode 100644 validphys2/src/validphys/asy_exponents.py
diff --git a/validphys2/examples/asy_exponents/asy_exponents_plots.yaml b/validphys2/examples/asy_exponents/asy_exponents_plots.yaml
new file mode 100644
index 0000000000..1a94fd35f3
--- /dev/null
+++ b/validphys2/examples/asy_exponents/asy_exponents_plots.yaml
@@ -0,0 +1,42 @@
+meta:
+ author: Emanuele R. Nocera
+ title: Asymptotic Exponents
+ keywords: [asy_exp]
+
+pdfs:
+ - {id: "NNPDF40_nnlo_as_01180", label: "NNPDF4.0"}
+
+Q: 5000 # GeV
+
+alpha:
+ xmin: 1e-5
+ xmax: 1e-3
+ ybottom: 0.0
+ ytop: 1.0
+ basis: flavour
+ npoints: 70
+
+beta:
+ xmin: 0.6
+ xmax: 0.9
+ ybottom: 0
+ ytop: 10
+ basis: flavour
+ npoints: 70
+
+show_mc_errors: False
+
+template_text: |
+ ---
+ title: Asymptotic exponents
+ keywords: [asy_exp]
+ ---
+ Asymptotic exponents plots
+ -----------
+ {@alpha plot_alpha_asy@}
+ {@beta plot_beta_asy@}
+ -----------
+
+actions_:
+ - report(main=true)
+
\ No newline at end of file
diff --git a/validphys2/examples/asy_exponents/asy_exponents_table.yaml b/validphys2/examples/asy_exponents/asy_exponents_table.yaml
new file mode 100644
index 0000000000..91baf90219
--- /dev/null
+++ b/validphys2/examples/asy_exponents/asy_exponents_table.yaml
@@ -0,0 +1,28 @@
+meta:
+ author: Emanuele R. Nocera
+ title: Asymptotic Exponents
+ keywords: [asy_exp]
+
+pdf: NNPDF40_nnlo_as_01180
+
+Q: 5000 # GeV
+
+x_alpha: 1e-5
+x_beta: 0.90
+npoints: 70
+
+basis: flavour
+
+template_text: |
+ ---
+ title: Asymptotic exponents
+ keywords: [asy_exp]
+ ---
+ Asymptotic exponents
+ -----------
+ {@asymptotic_exponents_table@}
+ -----------
+
+actions_:
+ - report(main=true)
+
\ No newline at end of file
diff --git a/validphys2/src/validphys/app.py b/validphys2/src/validphys/app.py
index 614022a9c7..4d40699df6 100644
--- a/validphys2/src/validphys/app.py
+++ b/validphys2/src/validphys/app.py
@@ -37,6 +37,7 @@ providers = [
"validphys.correlations",
"validphys.chi2grids",
"validphys.eff_exponents",
+ "validphys.asy_exponents",
"validphys.paramfits.dataops",
"validphys.paramfits.plots",
"validphys.theorycovariance.construction",
diff --git a/validphys2/src/validphys/asy_exponents.py b/validphys2/src/validphys/asy_exponents.py
new file mode 100644
index 0000000000..4add5e82b2
--- /dev/null
+++ b/validphys2/src/validphys/asy_exponents.py
@@ -0,0 +1,266 @@
+# -*- coding: utf-8 -*-
+"""
+Tools for computing and plotting effective exponents.
+"""
+from __future__ import generator_stop
+
+import logging
+import numbers
+import random
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from reportengine import collect
+from reportengine.compat import yaml
+from reportengine.figure import figuregen
+from reportengine.floatformatting import format_number, significant_digits, format_error_value_columns
+from reportengine.table import table
+
+from validphys.checks import check_positive, check_pdf_normalize_to, make_argcheck, check_xlimits
+from validphys.core import PDF, FitSpec
+from validphys.pdfbases import check_basis, Basis
+from validphys.pdfplots import BandPDFPlotter, PDFPlotter
+
+import validphys.pdfgrids as pdfgrids
+
+from findiff import FinDiff
+
+log = logging.getLogger(__name__)
+
+INTERNAL_LINESTYLE = ['-.', ':']
+INTERNAL_COLOR = plt.rcParams['axes.prop_cycle'].by_key()["color"]
+
+@check_positive('Q')
+@make_argcheck(check_basis)
+@check_xlimits
+def alpha_asy(pdf: PDF, *,
+ xmin: numbers.Real = 1e-6,
+ xmax: numbers.Real = 1e-3,
+ npoints: int = 100,
+ Q: numbers.Real = 1.65,
+ basis: (str, Basis),
+ flavours: (list, tuple, type(None)) = None):
+ """Returns a list of xplotting_grids containing the value of the asymptotic
+ exponent alpha, as defined by the first relationship in Eq. (4) of
+ [arXiv:1604.00024], at the specified value of Q (in GeV), in the interval [xmin, xmax].
+
+ basis: Is one of the bases defined in pdfbases.py. This includes 'flavour'
+ and 'evolution'.
+
+ flavours: A set of elements from the basis.
+ If None, the defaults for that basis will be selected.
+
+ npoints: the number of sub-intervals in the range [xmin, xmax] on which the
+ derivative is computed.
+ """
+ #Loading the filter map of the fit/PDF
+ checked = check_basis(basis, flavours)
+ basis = checked['basis']
+ flavours = checked['flavours']
+
+ if npoints == 2:
+ xGrid = np.array([xmin, xmax])
+ else:
+ xGrid = pdfgrids.xgrid(xmin, xmax, 'log', npoints)
+
+ pdfGrid = pdfgrids.xplotting_grid(
+ pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours)
+ pdfGrid_values = pdfGrid.grid_values.data
+ # NOTE: without this I get "setting an array element with a sequence"
+ xGrid = pdfGrid.xgrid
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', RuntimeWarning)
+ dx = np.log(xGrid[1]) - np.log(xGrid[0])
+ d_dx = FinDiff(2,dx,acc=4)
+ alphaGrid_values = -np.log(abs(pdfGrid_values))
+ alphaGrid_values = d_dx(alphaGrid_values)
+ alphaGrid_values[alphaGrid_values == -np.inf] = np.nan # when PDF_i =0
+ alphaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(alphaGrid_values))
+
+ return alphaGrid
+
+@check_positive('Q')
+@make_argcheck(check_basis)
+@check_xlimits
+def beta_asy(pdf, *,
+ xmin: numbers.Real = 0.6,
+ xmax: numbers.Real = 0.9,
+ npoints: int = 100,
+ Q: numbers.Real = 1.65,
+ basis: (str, Basis),
+ flavours: (list, tuple, type(None)) = None):
+ """Returns a list of xplotting_grids containing the value of the asymptotic
+ exponent beta, as defined by the second relationship in Eq. (4) of
+ [arXiv:1604.00024], at the specified value of Q (in GeV), in the interval [xmin, xmax].
+
+ basis: Is one of the bases defined in pdfbases.py. This includes 'flavour'
+ and 'evolution'.
+
+ flavours: A set of elements from the basis.
+ If None, the defaults for that basis will be selected.
+
+ npoints: the number of sub-intervals in the range [xmin, xmax] on which the
+ derivative is computed.
+ """
+ checked = check_basis(basis, flavours)
+ basis = checked['basis']
+ flavours = checked['flavours']
+
+ if npoints == 2:
+ xGrid = np.array([xmin, xmax])
+ else:
+ xGrid = pdfgrids.xgrid(xmin, xmax, 'linear', npoints)
+
+
+ pdfGrid = pdfgrids.xplotting_grid(
+ pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours)
+ pdfGrid_values = pdfGrid.grid_values.data
+ # NOTE: without this I get "setting an array element with a sequence"
+ xGrid = pdfGrid.xgrid
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', RuntimeWarning)
+ dx = xGrid[1] - xGrid[0]
+ d_dx = FinDiff(2,dx,acc=4)
+ betaGrid_values = np.log(abs(pdfGrid_values))
+ betaGrid_values = (xGrid - 1.) * d_dx(betaGrid_values)
+ betaGrid_values[betaGrid_values == -np.inf] = np.nan # when PDF_i =0
+ betaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(betaGrid_values))
+
+ return betaGrid
+
+class AsyExponentPlotter(PDFPlotter):
+ """ Class inherenting from BandPDFPlotter, changing title and ylabel to reflect the effective
+ exponent being plotted.
+ """
+
+ def __init__(self, exponent, *args, **kwargs):
+ self.exponent = exponent
+ super().__init__(*args, **kwargs)
+
+ def get_title(self, parton_name):
+ return fr"$\{self.exponent}_a$ for ${parton_name}$ at {format_number(self.Q, 3)} GeV"
+
+ def get_ylabel(self, parton_name):
+ if self.normalize_to is not None:
+ return "Ratio to {}".format(self.normalize_pdf.label)
+ else:
+ return fr"$\{self.exponent}_a$ for ${parton_name}$"
+
+class AsyExponentBandPlotter(BandPDFPlotter, AsyExponentPlotter):
+ def __init__(self, exponent, *args, **kwargs):
+ super().__init__(exponent, *args, **kwargs)
+
+alpha_asy_pdfs = collect('alpha_asy', ('pdfs',))
+
+@figuregen
+@check_pdf_normalize_to
+def plot_alpha_asy(
+ pdfs, alpha_asy_pdfs, pdfs_alpha_lines,
+ normalize_to: (int, str, type(None)) = None,
+ ybottom=None, ytop=None):
+ """ Plots the alpha asymptotic exponent """
+ yield from AsyExponentBandPlotter(
+ 'alpha',
+ pdfs,
+ alpha_asy_pdfs,
+ 'log',
+ normalize_to,
+ ybottom,
+ ytop,)
+
+beta_asy_pdfs = collect('beta_asy', ('pdfs',))
+
+@figuregen
+@check_pdf_normalize_to
+def plot_beta_asy(
+ pdfs, beta_asy_pdfs, pdfs_beta_lines,
+ normalize_to: (int, str, type(None)) = None,
+ ybottom=None, ytop=None):
+ """ Plots the beta asymptotic exponent """
+ yield from AsyExponentBandPlotter(
+ 'beta',
+ pdfs,
+ beta_asy_pdfs,
+ 'linear',
+ normalize_to,
+ ybottom,
+ ytop,)
+
+@table
+@make_argcheck(check_basis)
+def asymptotic_exponents_table(
+ pdf: PDF,
+ *,
+ x_alpha: numbers.Real = 1e-6,
+ x_beta: numbers.Real = 0.90,
+ Q: numbers.Real = 1.65,
+ basis:(str, Basis),
+ flavours: (list, tuple, type(None)) = None,
+ npoints=100,):
+ """ Returns a table with the values of the asymptotic exponents alpha and beta, as defined
+ in Eq. (4) of [arXiv:1604.00024], at the specified value of x and Q.
+
+ basis: Is one of the bases defined in pdfbases.py. This includes 'flavour'
+ and 'evolution'.
+
+ flavours: A set of elements from the basis.
+ If None, the defaults for that basis will be selected.
+
+ npoints: the number of sub-intervals in the range [xmin, xmax] on which the
+ derivative is computed.
+ """
+
+ alpha_a = alpha_asy(
+ pdf,
+ xmin=x_alpha,
+ xmax=1e-3,
+ npoints=npoints,
+ Q=Q,
+ basis=basis,
+ flavours=flavours)
+
+ beta_a = beta_asy(
+ pdf,
+ xmin=0.60,
+ xmax=x_beta,
+ npoints=npoints,
+ Q=Q,
+ basis=basis,
+ flavours=flavours)
+
+ alphastats = alpha_a.grid_values
+ betastats = beta_a.grid_values
+
+ with warnings.catch_warnings():
+ warnings.simplefilter('ignore', RuntimeWarning)
+
+ alpha_cv = alphastats.central_value()
+ beta_cv = betastats.central_value()
+
+ alpha_er = alphastats.std_error()
+ beta_er = betastats.std_error()
+
+ alpha_68 = alphastats.errorbar68()
+ beta_68 = betastats.errorbar68()
+
+ flavours_label = []
+ asy_exp_mean = []
+ asy_exp_err = []
+ asy_exp_min = []
+ asy_exp_max = []
+
+ for (j, fl) in enumerate(flavours):
+ asy_exp_mean.extend((alpha_cv[j,0],beta_cv[j,-1]))
+ asy_exp_err.extend((alpha_er[j,0],beta_er[j,-1]))
+ asy_exp_min.extend((alpha_68[0][j,0],beta_68[0][j,-1]))
+ asy_exp_max.extend((alpha_68[1][j,0],beta_68[1][j,-1]))
+ flavours_label.append(f"${basis.elementlabel(fl)}$")
+
+ asy_exp_data = {"mean": asy_exp_mean, "std": asy_exp_err, "min(68% CL)": asy_exp_min, "max(68% CL)": asy_exp_max}
+ ind = pd.MultiIndex.from_product([flavours_label, [r"$\alpha$", r"$\beta$"]])
+
+ df = pd.DataFrame(asy_exp_data, index=ind, columns=["mean","std","min(68% CL)","max(68% CL)"])
+ return df
diff --git a/validphys2/src/validphys/pdfbases.py b/validphys2/src/validphys/pdfbases.py
index 89321f098c..6154b59325 100644
--- a/validphys2/src/validphys/pdfbases.py
+++ b/validphys2/src/validphys/pdfbases.py
@@ -601,6 +601,28 @@ r'\bar{d}': {'dbar':1},
'c': {'c':1},
})
+@scalar_function_transformation(label="u_V")
+def u_valence(func, xmat, qmat):
+ gv = func([2, -2], xmat, qmat)
+ u = gv[:, [0], ...]
+ ubar = gv[:, [1], ...]
+ return u - ubar
+
+@scalar_function_transformation(label="d_V")
+def d_valence(func, xmat, qmat):
+ gv = func([1, -1], xmat, qmat)
+ d = gv[:, [0], ...]
+ dbar = gv[:, [1], ...]
+ return d - dbar
+
+@scalar_function_transformation(label="S")
+def total_sea(func, xmat, qmat):
+ gv = func([3, -2, -1, -3], xmat, qmat)
+ s = gv[:, [0], ...]
+ ubar = gv[:, [1], ...]
+ dbar = gv[:, [2], ...]
+ sbar = gv[:, [3], ...]
+ return 2.*(ubar + dbar) + s + sbar
@scalar_function_transformation(label="u/d")
def ud_ratio(func, xmat, qmat):
--
GitLab