Swap libNNPDF LHAPDFSet for a python PDF object

doc/sphinx/source/vp/examples.rst

@@ -69,6 +69,6 @@ NNLO pdf                             NNPDF40_nnlo_as_01180               NNPDF4.
 NNLO pdf hessian                     NNPDF40_nnlo_as_01180_hessian       NNPDF4.0 NNLO hessian PDF set generated from replicas
 NLO fit                              NNPDF40_nlo_as_01180                NNPDF4.0 NLO fit with 100 replicas (+ central replica)
 NNLO fit                             NNPDF40_nnlo_lowprecision           NNPDF4.0 NNLO low precision fit (theory 162) with 50 replicas (+ central replica)
-NNLO fit (iterated)                  NNPDF40_nnlo_lowprecision_iterated  Iteration of NNPDF40_nnlo_lowprecision
+NNLO fit (iterated)                  NNPDF40_nnlo_low_precision_iterated  Iteration of NNPDF40_nnlo_lowprecision
 fit                                  NNPDF40_example_closure_test        ``n3fit`` closure test fit with 30 replicas before and after postfit
 ===================================  =================================== ==================================================================

doc/sphinx/source/vp/pydataobjs.rst

@@ -255,3 +255,20 @@ from the API::
         "theoryid": 162
     }
     total_cov = API.dataset_inputs_covmat_from_systematics(**inp)
+
+Loading LHAPDF PDFs
+-------------------
+
+A wrapper class for LHAPDF PDFs is implemented in the :py:mod:`validphys.lhapdfset` module.
+An instance of this module will provide with a handful of useful wrappers to the underlying
+LHAPDF python interface. This is also the output of the ``pdf.load()`` method.
+
+For example, the following will return the values for all 100 members of NNPDF4.0 for
+the gluon and the d-quark, at three values of ``x`` at ``Q=91.2``.
+
+.. code-block:: python
+    from validphys.api import API
+    pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180")
+    l_pdf = pdf.load()
+    alpha_s = l_pdf.central_member.alphasQ(91.2)
+    results = l_pdf.grid_values([21,1], [0.1, 0.2, 0.3], [91.2])

validphys2/src/validphys/app.py

@@ -15,6 +15,7 @@
 import contextlib
 
 
+import lhapdf
 from reportengine import app
 
 from validphys.config import Config, Environment
@@ -131,7 +132,7 @@ def init(self):
             import NNPDF
 
             NNPDF.SetVerbosity(0)
-            NNPDF.SetLHAPDFVerbosity(0)
+            lhapdf.setVerbosity(0)
 
     @staticmethod
     def upload_context(do_upload, output):

validphys2/src/validphys/checks.py

@@ -39,7 +39,7 @@ def check_pdf_is_montecarlo(ns, **kwargs):
 def check_know_errors(ns, **kwargs):
     pdf = ns['pdf']
     try:
-        pdf.nnpdf_error
+        pdf.stats_class
     except NotImplementedError as e:
         raise CheckError(e) from e
 

validphys2/src/validphys/config.py

@@ -161,7 +161,7 @@ def parse_pdf(self, name: str):
 
         # Check that we know how to compute errors
         try:
-            pdf.nnpdf_error
+            pdf.stats_class
         except NotImplementedError as e:
             raise ConfigError(str(e))
         return pdf

validphys2/src/validphys/core.py

@@ -23,7 +23,7 @@
 from reportengine.baseexceptions import AsInputError
 from reportengine.compat import yaml
 
-from NNPDF import (LHAPDFSet,
+from NNPDF import (LHAPDFSet as libNNPDF_LHAPDFSet,
     CommonData,
     FKTable,
     FKSet,
@@ -38,23 +38,10 @@
 from validphys.theorydbutils import fetch_theory
 from validphys.hyperoptplot import HyperoptTrial
 from validphys.utils import experiments_to_dataset_inputs
+from validphys.lhapdfset import LHAPDFSet
 
 log = logging.getLogger(__name__)
 
-
-#TODO: Remove this eventually
-#Bacward compatibility error type names
-#Swig renamed these for no reason whatsoever.
-try:
-    LHAPDFSet.erType_ER_EIG
-except AttributeError:
-    import warnings
-    warnings.warn("libnnpdf out of date. Setting backwards compatible names")
-    LHAPDFSet.erType_ER_MC = LHAPDFSet.ER_MC
-    LHAPDFSet.erType_ER_EIG = LHAPDFSet.ER_EIG
-    LHAPDFSet.erType_ER_EIG90 = LHAPDFSet.ER_EIG90
-    LHAPDFSet.erType_ER_SYMEIG = LHAPDFSet.ER_SYMEIG
-
 class TupleComp:
 
     @classmethod
@@ -97,12 +84,25 @@ class PDF(TupleComp):
     Statistical estimators which depends on the PDF type (MC, Hessian...)
     are exposed as a :py:class:`Stats` object through the :py:attr:`stats_class` attribute
     The LHAPDF metadata can directly be accessed through the :py:attr:`info` attribute
+
+
+    Examples
+    --------
+    >>> from validphys.api import API
+    >>> from validphys.convolution import predictions
+    >>> args = {"dataset_input":{"dataset": "ATLASTTBARTOT"}, "theoryid":162, "use_cuts":"internal"}
+    >>> ds = API.dataset(**args)
+    >>> pdf = API.pdf(pdf="NNPDF40_nnlo_as_01180")
+    >>> preds = predictions(ds, pdf)
+    >>> preds.shape
+    (3, 100)
     """
 
     def __init__(self, name):
         self.name = name
         self._plotname = name
         self._info = None
+        self._stats_class = None
         super().__init__(name)
 
     @property
@@ -116,11 +116,15 @@ def label(self, label):
     @property
     def stats_class(self):
         """Return the stats calculator for this error type"""
-        error = self.error_type
-        klass = STAT_TYPES[error]
-        if self.error_conf_level is not None:
-            klass = functools.partial(klass, rescale_factor=self._rescale_factor())
-        return klass
+        if self._stats_class is None:
+            try:
+                klass = STAT_TYPES[self.error_type]
+            except KeyError:
+                raise NotImplementedError(f"No Stats class for error type {self.error_type}")
+            if self.error_conf_level is not None:
+                klass = functools.partial(klass, rescale_factor=self._rescale_factor())
+            self._stats_class = klass
+        return self._stats_class
 
     @property
     def infopath(self):
@@ -181,43 +185,46 @@ def _rescale_factor(self):
 
     @functools.lru_cache(maxsize=16)
     def load(self):
-        return LHAPDFSet(self.name, self.nnpdf_error)
+        return LHAPDFSet(self.name, self.error_type)
 
     @functools.lru_cache(maxsize=2)
     def load_t0(self):
         """Load the PDF as a t0 set"""
-        return LHAPDFSet(self.name, LHAPDFSet.erType_ER_MCT0)
+        return LHAPDFSet(self.name, "t0")
 
     def __str__(self):
         return self.label
 
     def __len__(self):
         return self.info["NumMembers"]
 
-    @property
-    def nnpdf_error(self):
-        """Return the NNPDF error tag, used to build the `LHAPDFSet` objeect"""
+    def legacy_load(self):
+        """Returns an libNNPDF LHAPDFSet object
+        Deprecated function used only in the `filter.py` module
+        """
         error = self.error_type
-        if error == "replicas":
-            return LHAPDFSet.erType_ER_MC
-
         cl = self.error_conf_level
-        if error == "hessian":
+        et = None
+        if error == "replicas":
+            et = libNNPDF_LHAPDFSet.erType_ER_MC
+        elif error == "hessian":
             if cl == 90:
-                return LHAPDFSet.erType_ER_EIG90
+                et = libNNPDF_LHAPDFSet.erType_ER_EIG90
             elif cl == 68:
-                return LHAPDFSet.erType_ER_EIG
+                et = libNNPDF_LHAPDFSet.erType_ER_EIG
             else:
                 raise NotImplementedError(f"No hessian errors with confidence interval {cl}")
-        if error == "symmhessian":
+        elif error == "symmhessian":
             if cl == 68:
-                return LHAPDFSet.erType_ER_SYMEIG
+                et = libNNPDF_LHAPDFSet.erType_ER_SYMEIG
             else:
                 raise NotImplementedError(
                     f"No symmetric hessian errors with confidence interval {cl}"
                 )
+        else:
+            raise NotImplementedError(f"Error type for {self}: '{error}' is not implemented")
 
-        raise NotImplementedError(f"Error type for {self}: '{error}' is not implemented")
+        return libNNPDF_LHAPDFSet(self.name, et)
 
     @property
     def grid_values_index(self):
@@ -236,17 +243,12 @@ def grid_values_index(self):
         -----
         The range object can be used efficiently as a Pandas index.
         """
-        err = self.nnpdf_error
-        if err is LHAPDFSet.erType_ER_MC:
+        if self.error_type == "replicas":
             return range(1, len(self))
-        elif err in (
-            LHAPDFSet.erType_ER_SYMEIG,
-            LHAPDFSet.erType_ER_EIG,
-            LHAPDFSet.erType_ER_EIG90,
-        ):
+        elif self.error_type in ("hessian", "symmhessian"):
             return range(0, len(self))
         else:
-            raise RuntimeError("Unknown error type")
+            raise RuntimeError(f"Unknown error type: {self.stats_class}")
 
     def get_members(self):
         """Return the number of members selected in ``pdf.load().grid_values``
@@ -822,7 +824,7 @@ def error_members(self):
 
 
 class SymmHessianStats(Stats):
-    """Compute stats in the 'assymetric' hessian format: The first index (0)
+    """Compute stats in the 'symetric' hessian format: The first index (0)
     is the
     central value. The rest of the indexes are results for each eigenvector.
     A 'rescale_factor is allowed in case the eigenvector confidence interval

validphys2/src/validphys/filters.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 Filters for NNPDF fits
 """
@@ -173,7 +172,7 @@ def _filter_closure_data(filter_path, data, fakepdfset, fakenoise, errorsize):
     """Filter closure test data."""
     total_data_points = 0
     total_cut_data_points = 0
-    fakeset = fakepdfset.load()
+    fakeset = fakepdfset.legacy_load()
     # Load data, don't cache result
     loaded_data = data.load.__wrapped__(data)
     # generate level 1 shift if fakenoise

validphys2/src/validphys/gridvalues.py

@@ -10,14 +10,8 @@
 
 import numpy as np
 
-from NNPDF import REALDOUBLE, LHAPDFSet
-
 from validphys.core import PDF
-
-#Numpy is unhappy about downcasting double to float implicitly, so we have
-#to manually make sure all input arrays correspond to NNPDF::real.
-FLTYPE = np.int32
-REALTYPE = np.double if REALDOUBLE else np.float32
+from validphys.lhapdfset import LHAPDFSet
 
 # Canonical ordering of PDG quark flavour codes
 QUARK_FLAVOURS = (-6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6)
@@ -54,9 +48,9 @@
 
 def _grid_values(lpdf, flmat, xmat, qmat):
     """Compute lpdf.grid_values with more forgiving argument types"""
-    flmat = np.atleast_1d(np.asanyarray(flmat, dtype=FLTYPE))
-    xmat = np.atleast_1d(np.asarray(xmat, dtype=REALTYPE))
-    qmat = np.atleast_1d(np.asarray(qmat, dtype=REALTYPE))
+    flmat = np.atleast_1d(np.asanyarray(flmat))
+    xmat = np.atleast_1d(np.asarray(xmat))
+    qmat = np.atleast_1d(np.asarray(qmat))
     return lpdf.grid_values(flmat, xmat, qmat)
 
 def grid_values(pdf:PDF, flmat, xmat, qmat):