"""
arclength.py
Module for the computation and presentation of arclengths.
"""
from collections import namedtuple
from collections.abc import Sequence
import numbers
import numpy as np
import pandas as pd
import scipy.integrate as integrate
from reportengine import collect
from reportengine.checks import check_positive, make_argcheck
from reportengine.figure import figure
from reportengine.table import table
from validphys import plotutils
from validphys.checks import check_pdf_normalize_to
from validphys.core import PDF
from validphys.pdfbases import Basis, check_basis
from validphys.pdfgrids import xgrid, xplotting_grid
ArcLengthGrid = namedtuple("ArcLengthGrid", ("pdf", "basis", "flavours", "stats"))
[docs]
@check_positive("Q")
@make_argcheck(check_basis)
def arc_lengths(
pdf: PDF,
Q: numbers.Real,
basis: (str, Basis) = "flavour",
flavours: (list, tuple, type(None)) = None,
):
"""
Compute arc lengths at scale Q
set up a grid with three segments and
compute the arclength for each segment.
Note: the variation of the PDF over the grid
is computed by computing the forward differences
between adjacent grid points.
Parameters
----------
pdf : validphys.core.PDF object
Q : float
scale at which to evaluate PDF
basis : default = "flavour"
flavours : default = None
Returns
-------
validphys.arclength.ArcLengthGrid object
object that contains the PDF, basis, flavours, and computed
arc length statistics.
"""
checked = check_basis(basis, flavours)
basis, flavours = checked["basis"], checked["flavours"]
# x-grid points and limits in three segments
npoints = 199 # 200 intervals
seg_min = [1e-6, 1e-4, 1e-2]
seg_max = [1e-4, 1e-2, 1.0]
res = np.zeros((pdf.get_members(), len(flavours)))
# Integrate the separate segments
for a, b in zip(seg_min, seg_max):
# Finite diff. step-size, x-grid
eps = (b - a) / npoints
ixgrid = xgrid(a, b, "linear", npoints)
# PDFs evaluated on grid, use the entire thing, the Stats class will chose later
xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data * ixgrid[1]
fdiff = np.diff(xfgrid) / eps # Compute forward differences
res += integrate.simpson(np.sqrt(1 + np.square(fdiff)), x=ixgrid[1][1:])
stats = pdf.stats_class(res)
return ArcLengthGrid(pdf, basis, flavours, stats)
# Collect arc_lengths over PDF list
pdfs_arc_lengths = collect(arc_lengths, ["pdfs"])
[docs]
@table
def arc_length_table(arc_lengths):
"""Return a table with the descriptive statistics of the arc lengths
over members of the PDF."""
arc_length_data = arc_lengths.stats.error_members()
arc_length_columns = [f"${arc_lengths.basis.elementlabel(fl)}$" for fl in arc_lengths.flavours]
return pd.DataFrame(arc_length_data, columns=arc_length_columns).describe().iloc[1:, :]
[docs]
@figure
@check_pdf_normalize_to
def plot_arc_lengths(
pdfs_arc_lengths: Sequence, Q: numbers.Real, normalize_to: (type(None), int) = None
):
"""Plot the arc lengths of provided pdfs"""
fig, ax = plotutils.subplots()
if normalize_to is not None:
ax.set_ylabel(f"Arc length $Q={Q}$ GeV (normalised)")
else:
ax.set_ylabel(f"Arc length $Q={Q}$ GeV")
for ipdf, arclengths in enumerate(pdfs_arc_lengths):
xvalues = np.array(range(len(arclengths.flavours)))
xlabels = ["$" + arclengths.basis.elementlabel(fl) + "$" for fl in arclengths.flavours]
ylower, yupper = arclengths.stats.errorbar68()
yvalues = (ylower + yupper) * 0.5
yerr = np.abs(yupper - ylower) * 0.5
if normalize_to is not None:
norm_cv = pdfs_arc_lengths[normalize_to].stats.central_value()
yvalues = np.divide(yvalues, norm_cv)
yerr = np.divide(yerr, norm_cv)
shift = (ipdf - (len(pdfs_arc_lengths) - 1) / 2.0) / 5.0
try:
ax.errorbar(
xvalues + shift,
yvalues,
yerr=yerr,
fmt='',
linestyle='',
label=arclengths.pdf.label,
)
except ValueError as e:
if any(yerr < 0):
raise ValueError(
f"The central value for pdf `{arclengths.pdf}` is outside of the error bands"
) from e
raise e
ax.set_xticks(xvalues)
ax.set_xticklabels(xlabels)
ax.legend()
return fig
# TODO: this should probably go somewhere else
[docs]
def integrability_number(
pdf: PDF,
Q: numbers.Real,
basis: (str, Basis) = "evolution",
flavours: (list, tuple, type(None)) = None,
):
r"""Return \sum_i |x_i*f(x_i)|, x_i = {1e-9, 1e-8, 1e-7}
for selected flavours
"""
checked = check_basis(basis, flavours)
basis, flavours = checked["basis"], checked["flavours"]
ixgrid = xgrid(1e-9, 1e-6, "log", 3)
xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data
res = np.sum(np.abs(xfgrid), axis=2)
return res.squeeze()