Source code for validphys.asy_exponents

# -*- coding: utf-8 -*-
"""
Tools for computing and plotting asymptotic exponents.
"""
import logging
import numbers
import warnings

import numpy as np
import pandas as pd

from reportengine import collect
from reportengine.checks import check_positive
from reportengine.figure import figuregen
from reportengine.floatformatting import format_number
from reportengine.table import table
from validphys.checks import check_pdf_normalize_to, check_xlimits, make_argcheck
from validphys.core import PDF
from validphys.pdfbases import Basis, check_basis
import validphys.pdfgrids as pdfgrids
from validphys.pdfplots import BandPDFPlotter, PDFPlotter

log = logging.getLogger(__name__)


[docs]@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]) alphaGrid_values = -np.log(abs(pdfGrid_values)) alphaGrid_values = np.gradient(alphaGrid_values, dx, axis=2, edge_order=2) 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
[docs]@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] betaGrid_values = np.log(abs(pdfGrid_values)) betaGrid_values = (xGrid - 1.0) * np.gradient(betaGrid_values, dx, axis=2, edge_order=2) 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
[docs]class AsyExponentBandPlotter(BandPDFPlotter): """Class inheriting from BandPDFPlotter, changing title and ylabel to reflect the asymptotic exponent being plotted. """ def __init__(self, exponent, *args, **kwargs): self.exponent = exponent super().__init__(*args, **kwargs)
[docs] def get_title(self, parton_name): return fr"$\{self.exponent}_a$ for ${parton_name}$ at {format_number(self.Q, 3)} GeV"
[docs] 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}$"
alpha_asy_pdfs = collect('alpha_asy', ('pdfs',))
[docs]@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',))
[docs]@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 )
[docs]@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