Source code for validphys.mc2hessian


This module containts the functionality to compute reduced set using the `mc2hessian` algorithm
(See section 2.1 of of `1602.00005 <>`_).
import logging
import numbers
import pathlib
import shutil

import numpy as np

from reportengine.checks import check, check_positive, make_argcheck
from validphys import lhaindex
from validphys.checks import check_pdf_is_montecarlo
from validphys.lhio import hessian_from_lincomb
from validphys.pdfgrids import xplotting_grid

log = logging.getLogger(__name__)

[docs]def gridname(pdf, Neig, mc2hname: (str, type(None)) = None): """If no custom `mc2hname' is specified, the name of the Hessian PDF is automatically generated.""" if mc2hname is None: grid_name = f"{}_hessian_{Neig}" else: grid_name = mc2hname return grid_name
@make_argcheck def _check_xminmax(xmin, xminlin, xmax): check(0 < xmin < xminlin < xmax <= 1, "Expecting 0 < xmin < xminlin < xmax <= 1")
[docs]@_check_xminmax @check_positive("nplog") @check_positive("nplin") def mc2hessian_xgrid( xmin: float = 1e-5, xminlin: float = 1e-1, xmax: numbers.Real = 1, nplog: int = 50, nplin: int = 50, ): """Provides the points in x to sample the PDF. `logspace` and `linspace` will be called with the respsctive parameters. Generates a grid with ``nplog`` logarithmically spaced points between ``xmin`` and ``xminlin`` followed by ``nplin`` linearly spaced points between ``xminlin`` and ``xmax`` """ return np.append( np.geomspace(xmin, xminlin, num=nplog, endpoint=False), np.linspace(xminlin, xmax, num=nplin, endpoint=False), )
[docs]@check_pdf_is_montecarlo def mc2hessian( pdf, Q, Neig: int, mc2hessian_xgrid, output_path, gridname, installgrid: bool = False ): """Produces a Hessian PDF by transfroming a Monte Carlo PDF set. Parameters ----------- pdf : validphys.core.PDF An existng validphys PDF object which will be converted into a Hessian PDF set Q : float Energy scale at which the Monte Carlo PDF is sampled Neig : int Number of basis eigenvectors in the Hessian PDF set mc2hessian_xgrid : numpy.ndarray The points in x at which to sample the Monte Carlo PDF set output path : pathlib.PosixPath The validphys output path where the PDF will be written gridname : str Name of the Hessian PDF set installgrid : bool, optional, default=``False`` Whether to copyt the Hessian grid to the LHAPDF path """ result_path = _create_mc2hessian( pdf, Q=Q, xgrid=mc2hessian_xgrid, Neig=Neig, output_path=output_path, name=gridname ) if installgrid: lhafolder = pathlib.Path(lhaindex.get_lha_datapath()) dest = lhafolder / gridname if lhaindex.isinstalled(gridname): log.warning( "Target directory for new PDF, %s, already exists. " "Removing contents.", dest, ) if dest.is_dir(): shutil.rmtree(str(dest)) else: dest.unlink() shutil.copytree(result_path, dest)"Hessian PDF set installed at %s", dest)
def _create_mc2hessian(pdf, Q, xgrid, Neig, output_path, name=None): X = _get_X(pdf, Q, xgrid, reshape=True) vec = _compress_X(X, Neig) norm = _pdf_normalization(pdf) return hessian_from_lincomb(pdf, vec / norm, folder=output_path, set_name=name) def _get_X(pdf, Q, xgrid, reshape=False): pdf_grid = xplotting_grid(pdf, Q, xgrid=xgrid) pdf_grid_values = pdf_grid.grid_values replicas = pdf_grid_values.error_members() mean = pdf_grid_values.central_value() Xt = replicas - mean if reshape: Xt = Xt.reshape(Xt.shape[0], Xt.shape[1] * Xt.shape[2]) return Xt.T def _compress_X(X, neig): _U, _S, V = np.linalg.svd(X, full_matrices=False) vec = V[:neig, :].T return vec def _pdf_normalization(pdf): """Extract the quantity by which we have to divide the eigenvectors to get the correct errors, depending on the `error_type` of `pdf`.""" nrep = len(pdf) - 1 if pdf.error_type == "replicas": norm = np.sqrt(nrep - 1) elif pdf.error_type in ("hessian", "symmhessian"): norm = 1 else: raise NotImplementedError( "This PDF error type is not supported." "PDF error: %s" % pdf.error_type ) return norm