From 34340df8a9c5da2c643198e18adb803fbdb28116 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Tue, 15 Feb 2022 16:16:23 +0100
Subject: [PATCH 01/13] make MC PDFs return all members (transitional commit
 for testing)

---
 n3fit/src/n3fit/tests/test_vpinterface.py   |  2 -
 n3fit/src/n3fit/vpinterface.py              |  9 +++-
 validphys2/src/validphys/core.py            | 25 +++------
 validphys2/src/validphys/correlations.py    |  7 +--
 validphys2/src/validphys/dataplots.py       |  2 +-
 validphys2/src/validphys/lhapdfset.py       |  5 --
 validphys2/src/validphys/sumrules.py        | 59 +++++++++++----------
 validphys2/src/validphys/tests/test_core.py |  3 +-
 8 files changed, 52 insertions(+), 60 deletions(-)

diff --git a/n3fit/src/n3fit/tests/test_vpinterface.py b/n3fit/src/n3fit/tests/test_vpinterface.py
index 255a2acb82..3664ee2134 100644
--- a/n3fit/src/n3fit/tests/test_vpinterface.py
+++ b/n3fit/src/n3fit/tests/test_vpinterface.py
@@ -6,7 +6,6 @@
 from hypothesis import given, settings, example
 from hypothesis.strategies import integers
 from validphys.pdfgrids import xplotting_grid, distance_grids
-from validphys.core import MCStats
 from n3fit.vpinterface import N3PDF
 from n3fit.model_gen import pdfNN_layer_generator
 
@@ -41,7 +40,6 @@ def test_N3PDF(members, layers):
     xx = np.random.rand(xsize)
     n3pdf = generate_n3pdf(layers=layers, members=members)
     assert len(n3pdf) == members + 1
-    assert n3pdf.stats_class == MCStats
     assert n3pdf.load() is n3pdf
     w = n3pdf.get_nn_weights()
     assert len(w) == members
diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py
index 98ccd92bed..9f53b786f2 100644
--- a/n3fit/src/n3fit/vpinterface.py
+++ b/n3fit/src/n3fit/vpinterface.py
@@ -45,6 +45,10 @@
     "T35",
 ]
 
+class N3Stats(MCStats):
+    def central_value(self):
+        return self.data[0]
+
 
 class N3PDF(PDF):
     """
@@ -76,10 +80,13 @@ def __init__(self, pdf_models, fit_basis=None, name="n3fit"):
         # to the average of all others
         self._info = {"ErrorType": "replicas", "NumMembers": len(self._models) + 1}
 
+    def get_members(self):
+        return max(len(self._models), 2)
+
     @property
     def stats_class(self):
         """The stats class of N3PDF is always MCStats"""
-        return MCStats
+        return N3Stats
 
     def load(self):
         """Many vp functions ask for a LHAPDF pdf
diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index 85536b7e23..bc76edf28a 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -243,19 +243,12 @@ def grid_values_index(self):
         -----
         The range object can be used efficiently as a Pandas index.
         """
-        if self.error_type == "replicas":
-            return range(1, len(self))
-        elif self.error_type in ("hessian", "symmhessian"):
-            return range(0, len(self))
-        else:
-            raise RuntimeError(f"Unknown error type: {self.stats_class}")
+        return range(0, len(self))
 
     def get_members(self):
         """Return the number of members selected in ``pdf.load().grid_values``
-        operation. See :py:meth:`PDF.grid_values_index` for details on differences
-        between types of PDF sets.
         """
-        return len(self.grid_values_index)
+        return len(self)
 
 
 kinlabels_latex = CommonData.kinLabel_latex.asdict()
@@ -780,7 +773,7 @@ def central_value(self):
         raise NotImplementedError()
 
     def error_members(self):
-        raise NotImplementedError()
+        return self.data[1:]
 
     def std_error(self):
         raise NotImplementedError()
@@ -809,14 +802,14 @@ class MCStats(Stats):
     """Result obtained from a Monte Carlo sample"""
 
     def central_value(self):
-        return np.mean(self.data, axis=0)
+        return np.mean(self.error_members(), axis=0)
 
     def std_error(self):
         # ddof == 1 to match libNNPDF behaviour
-        return np.std(self.data, ddof=1, axis=0)
+        return np.std(self.error_members(), ddof=1, axis=0)
 
     def moment(self, order):
-        return np.mean(np.power(self.data-self.central_value(),order), axis=0)
+        return np.mean(np.power(self.error_members()-self.central_value(),order), axis=0)
 
     def errorbar68(self):
         #Use nanpercentile here because we can have e.g. 0/0==nan normalization
@@ -828,9 +821,6 @@ def errorbar68(self):
     def sample_values(self, size):
         return np.random.choice(self, size=size)
 
-    def error_members(self):
-        return self.data
-
 
 class SymmHessianStats(Stats):
     """Compute stats in the 'symetric' hessian format: The first index (0)
@@ -848,9 +838,6 @@ def central_value(self):
     def errorbar68(self):
         return self.errorbarstd()
 
-    def error_members(self):
-        return self.data[1:]
-
     def std_error(self):
         data = self.data
         diffsq = (data[0] - data[1:])**2
diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py
index 7b1f6b8920..edafd89fef 100644
--- a/validphys2/src/validphys/correlations.py
+++ b/validphys2/src/validphys/correlations.py
@@ -9,8 +9,7 @@
 
 from reportengine import collect
 
-from validphys.checks import check_pdf_is_montecarlo
-
+from validphys.core import Stats
 
 #This would be a good candidate to be optimized to calculate everything in one
 #pass over x,
@@ -62,7 +61,9 @@ def obs_pdf_correlations(pdf, results, xplotting_grid):
     """Return the correlations between each point in a dataset and the PDF
     values on a grid of (x,f) points in a format similar to `xplotting_grid`."""
     _, th = results
-    corrs = pdf.stats_class(_basic_obs_pdf_correlation(xplotting_grid.grid_values, th.stats))
+    # Wrap the result in a standard Stats class
+    # since the result is (npoints, flavours, ndata) and has nothing to do with the PDF replicas
+    corrs = Stats(_basic_obs_pdf_correlation(xplotting_grid.grid_values, th.stats))
     return xplotting_grid.copy_grid(grid_values=corrs)
 
 
diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py
index 725fb3ec65..3a7460ab05 100644
--- a/validphys2/src/validphys/dataplots.py
+++ b/validphys2/src/validphys/dataplots.py
@@ -871,7 +871,7 @@ def plot_smpdf(pdf, dataset, obs_pdf_correlations, mark_threshold:float=0.9):
 
     basis = obs_pdf_correlations.basis
 
-    fullgrid = obs_pdf_correlations.grid_values.error_members()
+    fullgrid = obs_pdf_correlations.grid_values.data
 
     fls = obs_pdf_correlations.flavours
     x = obs_pdf_correlations.xgrid
diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py
index ff029dec7a..8b52a24c90 100644
--- a/validphys2/src/validphys/lhapdfset.py
+++ b/validphys2/src/validphys/lhapdfset.py
@@ -43,9 +43,6 @@ class LHAPDFSet:
 
     Once instantiated this class will load the PDF set from LHAPDF.
     If it is a T0 set only the CV will be loaded.
-
-    For Monte Carlo sets the central value (member 0) is by default not included when taking
-    the results for all members (i.e., when using ``grid_values``).
     """
 
     def __init__(self, name, error_type):
@@ -78,8 +75,6 @@ def members(self):
         """
         if self.is_t0:
             return self._lhapdf_set[0:1]
-        if self._error_type == "replicas":
-            return self._lhapdf_set[1:]
         return self._lhapdf_set
 
     @property
diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py
index 702ef350fd..17435a0440 100644
--- a/validphys2/src/validphys/sumrules.py
+++ b/validphys2/src/validphys/sumrules.py
@@ -11,19 +11,19 @@
 
 import numpy as np
 import pandas as pd
-import scipy.integrate as integrate
+from scipy.integrate import quad
 
 from reportengine.table import table
 from reportengine.checks import check_positive
 from reportengine.floatformatting import format_error_value_columns
 
 from validphys.core import PDF
-from validphys.pdfbases import ALL_FLAVOURS, parse_flarr
-from validphys.lhapdfset import LHAPDFSet
+from validphys.pdfbases import parse_flarr
 
 
-def _momentum_sum_rule_integrand(x, lpdf:LHAPDFSet, irep, Q):
-    return sum([lpdf.xfxQ(x, Q=Q, n=irep, fl=fl) for fl in ALL_FLAVOURS])
+def _momentum_sum_rule_integrand(x, lpdf, Q):
+    res = lpdf.xfxQ(x, Q)
+    return sum([res[f] for f in lpdf.flavors()])
 
 def _make_momentum_fraction_integrand(fldict):
     """Make a suitable integrand function, which takes x to be integrated over
@@ -49,9 +49,9 @@ def _make_momentum_fraction_integrand(fldict):
     # Do this outside to aid integration time
     fldict = {parse_flarr([k])[0]: v for k, v in fldict.items()}
 
-    def f(x, lpdf, irep, Q):
+    def f(x, lpdf, Q):
         return sum(
-            multiplier * lpdf.xfxQ(x, Q=Q, n=irep, fl=flavour)
+            multiplier * lpdf.xfxQ(x, Q)[flavour]
             for flavour, multiplier in fldict.items()
         )
 
@@ -81,10 +81,10 @@ def _make_pdf_integrand(fldict):
     # Do this outsde to aid integration time
     fldict = {parse_flarr([k])[0]: v for k, v in fldict.items()}
 
-    def f(x, lpdf, irep, Q):
+    def f(x, lpdf, Q):
         return (
             sum(
-                multiplier * lpdf.xfxQ(x, Q=Q, n=irep, fl=flavour)
+                multiplier * lpdf.xfxQ(x, Q)[flavour]
                 for flavour, multiplier in fldict.items()
             )
             / x
@@ -123,24 +123,23 @@ def f(x, lpdf, irep, Q):
 }
 
 
+def _integral(rule_f, pdf_member, Q, config=None):
+    """Integrate `rule_f` for a given `pdf_member` at a given energy
+    separating the regions of integration. Uses quad.
+    """
+    if config is None:
+        config = {"limit":1000, "epsabs": 1e-4, "epsrel": 1e-4}
+    res = 0.0
+    lims = [(1e-9, 1e-5), (1e-5, 1e-3), (1e-3, 1)]
+    for lim in lims:
+        res += quad(rule_f, *lim, args=(pdf_member, Q), **config)[0]
+    return res
+
+
 def _sum_rules(rules_dict, lpdf, Q):
     """Compute a SumRulesGrid from the loaded PDF, at Q"""
-    nmembers = lpdf.n_members
-    #TODO: Write this in something fast
-    #If nothing else, at least allocate and store the result contiguously
-    res = np.zeros((len(rules_dict), nmembers))
-    integrands = rules_dict.values()
-    def integral(f, a, b, irep):
-        #We increase the limit to capture the log scale fluctuations
-        return integrate.quad(f, a, b, args=(lpdf, irep, Q),
-               limit=1000,
-               epsabs=1e-4, epsrel=1e-4)[0]
-
-    for irep in range(nmembers):
-        for r, f in enumerate(integrands):
-            res[r,irep] =  (integral(f, 1e-9, 1e-5, irep)  +
-               integral(f, 1e-5, 1e-3, irep) + integral(f, 1e-3, 1, irep))
-    return {k: v for k,v in zip(rules_dict, res)}
+    return {k: [_integral(r, m, Q) for m in lpdf.members] for k,r in rules_dict.items()}
+    # I very much doubt that looping over PDF members is the slow piece here...
 
 
 @check_positive('Q')
@@ -150,7 +149,10 @@ def sum_rules(pdf:PDF, Q:numbers.Real):
     SumRulesGrid object with the list of values for each sum rule.  The
     integration is performed with absolute and relative tolerance of 1e-4."""
     lpdf = pdf.load()
-    return _sum_rules(KNOWN_SUM_RULES, lpdf, Q)
+    tmp = _sum_rules(KNOWN_SUM_RULES, lpdf, Q)
+    if pdf.error_type == "replicas":
+        return {k: v[1:] for k,v in tmp.items()}
+    return tmp
 
 @check_positive('Q')
 def central_sum_rules(pdf:PDF, Q:numbers.Real):
@@ -175,7 +177,10 @@ def unknown_sum_rules(pdf: PDF, Q: numbers.Real):
        - T8
     """
     lpdf = pdf.load()
-    return _sum_rules(UNKNOWN_SUM_RULES, lpdf, Q)
+    tmp = _sum_rules(UNKNOWN_SUM_RULES, lpdf, Q)
+    if pdf.error_type == "replicas":
+        return {k: v[1:] for k,v in tmp.items()}
+    return tmp
 
 def _simple_description(d):
     res = {}
diff --git a/validphys2/src/validphys/tests/test_core.py b/validphys2/src/validphys/tests/test_core.py
index 5dfefc705d..17ca8b3c4c 100644
--- a/validphys2/src/validphys/tests/test_core.py
+++ b/validphys2/src/validphys/tests/test_core.py
@@ -16,9 +16,8 @@ def test_pdf(pdf_name):
     assert pdf.isinstalled
     error_type = pdf.error_type
     if error_type == "replicas":
-        assert pdf.get_members() == (len(pdf)-1)
         assert pdf.error_conf_level is None
     else:
-        assert pdf.get_members() == len(pdf)
         assert isinstance(pdf.error_conf_level, (int, float))
+    assert pdf.get_members() == len(pdf)
     assert pdf.name == pdf._plotname == pdf_name == str(pdf)

From 907ae940f072d21fda5fa7841505ae9037db0949 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Tue, 15 Feb 2022 16:40:16 +0100
Subject: [PATCH 02/13] update regression test for sumrules and remove
 unnecesary if condition

---
 validphys2/src/validphys/sumrules.py           | 18 +++++-------------
 .../tests/regressions/test_sum_rules_MC.csv    | 14 +++++++-------
 .../regressions/test_sum_rules_hessian.csv     |  8 ++++----
 3 files changed, 16 insertions(+), 24 deletions(-)

diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py
index 17435a0440..05de4d90c6 100644
--- a/validphys2/src/validphys/sumrules.py
+++ b/validphys2/src/validphys/sumrules.py
@@ -27,8 +27,7 @@ def _momentum_sum_rule_integrand(x, lpdf, Q):
 
 def _make_momentum_fraction_integrand(fldict):
     """Make a suitable integrand function, which takes x to be integrated over
-    and fixed loaded PDF, replica number and Q that computes the momentum
-    fraction based on ``fldict``.
+    and a PDF member and Q that computes the momentum fraction based on ``fldict``.
 
     The keys of ``fldict`` are free form values corresponding to PDG parton ids
     (that end up being passed by :py:func:`validphys.pdfbases.parse_flarr` and
@@ -59,8 +58,7 @@ def f(x, lpdf, Q):
 
 def _make_pdf_integrand(fldict):
     """Make a suitable integrand function, which takes x to be integrated over
-    and fixed loaded PDF, replica number and Q that computes the integrand of
-    the PDFs based on ``fldict``.
+    and a PDF member and Q that computes the integrand of the PDFs based on ``fldict``.
 
     The keys of ``fldict`` are free form values corresponfing to PDG parton ids
     (that end up being passed :py:func:`validphys.pdfbases.parse_flarr` and
@@ -149,10 +147,8 @@ def sum_rules(pdf:PDF, Q:numbers.Real):
     SumRulesGrid object with the list of values for each sum rule.  The
     integration is performed with absolute and relative tolerance of 1e-4."""
     lpdf = pdf.load()
-    tmp = _sum_rules(KNOWN_SUM_RULES, lpdf, Q)
-    if pdf.error_type == "replicas":
-        return {k: v[1:] for k,v in tmp.items()}
-    return tmp
+    return _sum_rules(KNOWN_SUM_RULES, lpdf, Q)
+
 
 @check_positive('Q')
 def central_sum_rules(pdf:PDF, Q:numbers.Real):
@@ -161,7 +157,6 @@ def central_sum_rules(pdf:PDF, Q:numbers.Real):
     return _sum_rules(KNOWN_SUM_RULES, lpdf, Q)
 
 
-
 @check_positive('Q')
 def unknown_sum_rules(pdf: PDF, Q: numbers.Real):
     """Compute the following integrals
@@ -177,10 +172,7 @@ def unknown_sum_rules(pdf: PDF, Q: numbers.Real):
        - T8
     """
     lpdf = pdf.load()
-    tmp = _sum_rules(UNKNOWN_SUM_RULES, lpdf, Q)
-    if pdf.error_type == "replicas":
-        return {k: v[1:] for k,v in tmp.items()}
-    return tmp
+    return _sum_rules(UNKNOWN_SUM_RULES, lpdf, Q)
 
 def _simple_description(d):
     res = {}
diff --git a/validphys2/src/validphys/tests/regressions/test_sum_rules_MC.csv b/validphys2/src/validphys/tests/regressions/test_sum_rules_MC.csv
index 1493b61173..e90c460b2e 100644
--- a/validphys2/src/validphys/tests/regressions/test_sum_rules_MC.csv
+++ b/validphys2/src/validphys/tests/regressions/test_sum_rules_MC.csv
@@ -1,15 +1,15 @@
 	mean	std	min	max	Central value
-momentum	1.0000359397024803	5.384464994542894e-07	1.000034619691424	1.0000372097117414	1.0000359463817983
-uvalence	1.9997929722405798	0.00028226750040615813	1.9987658750095758	2.000123119074607	1.9997934800108916
-dvalence	0.9997171906346781	0.00028320184172691307	0.9986900055519239	1.0000492249686646	0.9997173761604233
-svalence	0.00024972763167990374	0.0003322524099365173	-0.00042330789669951796	0.00144823410093127	0.0002510958414509573
+momentum	1.0000359402582584	5.358395010408105e-07	1.0000346238340416	1.0000372086565177	1.0000359458599606
+uvalence	1.9997926747574069	0.0002809483175437137	1.998765120611357	2.0001202779127008	1.9997937994497308
+dvalence	0.9997164295681845	0.00028118151759189136	0.9986902415578234	1.000046111974624	0.9997173113369415
+svalence	0.00025025303862949585	0.0003303956500365186	-0.00042417491773433236	0.0014490597353980783	0.0002507254062854847
 u momentum fraction	0.2555	0.0021	-1.0	-1.0	-1.0
-ubar momentum fraction	3.387E-2	7.5E-4	-1.0	-1.0	-1.0
+ubar momentum fraction	0.03387	0.00074	-1.0	-1.0	-1.0
 d momentum fraction	0.1278	0.0021	-1.0	-1.0	-1.0
 dbar momentum fraction	0.0422	0.0012	-1.0	-1.0	-1.0
-s momentum fraction	0.0316	0.0020	-1.0	-1.0	-1.0
+s momentum fraction	0.0316	0.002	-1.0	-1.0	-1.0
 sbar momentum fraction	0.0253	0.0012	-1.0	-1.0	-1.0
 cp momentum fraction	0.0286	0.0017	-1.0	-1.0	-1.0
 g momentum fraction	0.4466	0.0017	-1.0	-1.0	-1.0
 T3	0.79	0.23	-1.0	-1.0	-1.0
-T8	4.73	0.70	-1.0	-1.0	-1.0
+T8	4.73	0.7	-1.0	-1.0	-1.0
diff --git a/validphys2/src/validphys/tests/regressions/test_sum_rules_hessian.csv b/validphys2/src/validphys/tests/regressions/test_sum_rules_hessian.csv
index 2cbd36e8e6..04168d5838 100644
--- a/validphys2/src/validphys/tests/regressions/test_sum_rules_hessian.csv
+++ b/validphys2/src/validphys/tests/regressions/test_sum_rules_hessian.csv
@@ -1,8 +1,8 @@
 	mean	std	min	max	Central value
-momentum	1.0000359343288223	3.262768320141462e-08	1.0000357928883208	1.0000360043344905	1.0000359428250625
-uvalence	1.999794587352472	3.771059740749569e-05	1.9996273397296627	1.9998792759145556	1.999804641339146
-dvalence	0.999719265014927	3.893093441227586e-05	0.9995458470781201	0.9998028245378152	0.9997279363904001
-svalence	0.00021669265424246442	4.864440927527348e-05	3.0026081505308416e-05	0.000370776050729299	0.0002126701429466457
+momentum	1.0000359346945835	3.249002390618474e-08	1.0000357949386856	1.000036007099919	1.000035943610909
+uvalence	1.9997947975241677	3.774745222364767e-05	1.9996255039155293	1.9998795897064394	1.9998052948252174
+dvalence	0.9997193904389039	3.913237872100671e-05	0.9995455646219639	0.9998033688013618	0.9997281952294002
+svalence	0.00021668009787772745	4.881864767300481e-05	2.8212616664191614e-05	0.00037062725233145966	0.00021200815810076823
 u momentum fraction	0.2554	0.00025	-1.0	-1.0	-1.0
 ubar momentum fraction	0.03386	0.0001	-1.0	-1.0	-1.0
 d momentum fraction	0.1277	0.00029	-1.0	-1.0	-1.0

From 178b8b5e6f233172cdfa8c30af1ae606cba94823 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Tue, 15 Feb 2022 17:06:04 +0100
Subject: [PATCH 03/13] use directly central member in MCStats -and reupdate
 regressions-

---
 validphys2/src/validphys/core.py              |   9 +--------
 .../tests/baseline/test_dataspecschi2.png     | Bin 23213 -> 23215 bytes
 .../tests/regressions/test_datasetchi2.csv    |   6 +++---
 .../regressions/test_replicachi2data.csv      |   2 +-
 4 files changed, 5 insertions(+), 12 deletions(-)

diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index bc76edf28a..adaa264fd9 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -770,7 +770,7 @@ def __init__(self, data):
         self.data = np.atleast_2d(data)
 
     def central_value(self):
-        raise NotImplementedError()
+        return self.data[0]
 
     def error_members(self):
         return self.data[1:]
@@ -800,10 +800,6 @@ def errorbarstd(self):
 
 class MCStats(Stats):
     """Result obtained from a Monte Carlo sample"""
-
-    def central_value(self):
-        return np.mean(self.error_members(), axis=0)
-
     def std_error(self):
         # ddof == 1 to match libNNPDF behaviour
         return np.std(self.error_members(), ddof=1, axis=0)
@@ -832,9 +828,6 @@ def __init__(self, data, rescale_factor=1):
         super().__init__(data)
         self.rescale_factor = rescale_factor
 
-    def central_value(self):
-        return self.data[0]
-
     def errorbar68(self):
         return self.errorbarstd()
 
diff --git a/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png b/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png
index 0115b2daa34ff91bd168412dc1607937cb407fbd..103835fd4ddd8938a0b4be7bac5ed91b1936dd14 100644
GIT binary patch
delta 13595
zcmZX*2~>^i8#nw^$lN)GA|w<V<QPiQM97ZjNohc$LGwhaCu7E<kp?9icB3W@JB1KZ
zsdhrMO6{Er6`J4Q-T!;P^{w@tb<SF6IqmlI-1l|;rt3QCc|4hUJh%OLR%X#mLBlv<
zY0DBbKK_LUHv~7G$oV|fm)|=1tsyitwCD4wSi{4I?<y)OEfo>T;M67Uirp+FRcAT*
zODNB-Uc<t|!mhRCsi5>>=jEk0A_M-cqU=Nh{V7#h`H<Bs#Od-W75+aX-LVZSQzP2G
zf7=7>nv3%M#(%0e`;GOPKK0I1@ce$Jup~_Q^7JR8tWCmG11SG5FHh#!_BRy@Z@2ld
zGvDO#(ewR9y)p?&Nm8PsZE1}&`1yP9ow{FBCZO`&=X=$zfzbT?eC3IuI<5%ob%$>N
zBTCLeue&Nn{IEkP`)z#(drK@sSn}4Qq@(w8Ewh!A@1Jg4B277;A1E2qGG}?#hf3%#
zS+pqr?%icd!)$#?Nl8~neKRxZ&6_hdEW7*q(!->TBSrE*w)yf?4IR^;3=vtmjAoxz
zdU0z<OEObGiCfg0p`U7$o@-OHw({OK^`%Rf%G@vc_I3Ibvu9HU{$=gAD^J#)iuFxQ
z%_!2?(dv_vc<5fdNzTcgD!1s~^!z2OR;7PTKT=Uvre^i}X3gVcWu4!?ef!F>muc~$
zoRxTIP-dbVHy(KY{JGty#}5Xx{{E%U>|D9f+S-~C9J$>#-l=rLXoL2BgZ=xj=;`UP
z<974y+q<{G?z?ahPSz}l;-Q3kTAs3g542`pxWL<OB_Sb^-PkvG&KzcLdPGFTg^L#r
zEi5v2x_#?;(&)HN$ua%;*g`5pom}O{jT?22ACF(Oc(FKx!{KNgI%IJ2WTF^t`{L4^
z(X!Q(ZB3Ic=LQTrUf+`S`Rm-*RH$DOvALm6TP6M9b^G5mHXaGib!=9^v6P=}eo?uj
zHSPNKxwlsAsQUZID*-3u9~nUL_}l%c*j$=yRUTO7vG#S~s?@+$ySSHyEN+Xfl6QVK
zaxhD5WOUT7{l(>X6&3msGUo2BzWe{!wmm!VJ=T{W)3`)PsJF2DzP^u7@vH^HkK*D)
z)6EL?ZVAh79ery4_Dm~J+_z(^{-G%Uc`p@2Gw_eha|KZff7z`KX<xIFcJUgfC}-x@
zZQH^^LOQrxxY5DEi_Us`7kS@Q*zn(fIu9N^$PnS9WQ2Km<Mvm)e{XR5boBKb^URer
znyo@!-fFrv`0KacuZ$CGXQ<vXOyxZ>@~g{xVz?pCrDH#yU)8O<sy=S_#i}t`adGj{
z3d_lOZIuZFegP2!9i7`2XWHqRx0Y?$6TJOY)yj!JC%W#w(m4Yoqr{#b8>NxXkfN!{
zu?0Trvit)|X5t$*bgTUQ$$TVl9HbvGekA*TME~t2JzNv+?Ul?kT3N2HuCxE4WR242
z-VWKIo#H$8TUgnrpuc#qO47IT_8K+!o)2q$$A6Yzzka>;JRg6jvwU<l{|tJ@bfcS@
zMRV&@lLktsEXo1|jB`#_vP<cyis5v&eDNRWr}A0L?lKrjlS?-W^72t;w@sId^39t6
z!JmgOrFi1P%z3%)9KHYbsA}*M5fKa9sFtVaj7JBpf|L?b=(ew}FCAYrE>-*o=bn6B
z!B;5l|281~|LZ`7-rwQYwl(qlG_<sgs_!c;&_2TY(ZWYjY|T9G>^GM+4j!y~7bex!
z)s=vn-MC?cI_a94hw8Gvqwd|F<w2{Bot+<RXlS5|5^d`ef6f^f|N4|4UmV4C_m^L}
z@bhUo-`}54k!BLJrQd~02yZ=c6=ke+_T0IH!2q7PmoN7&U9!Xm&mg?ZGi~j<b=C6w
zW*e5D6QC%a9rYY)&fK|r_>MuLqm1i8K@nx`Iz~p}@ed!0Gdx-5y1Kf-j~1>{*eD^v
z2wAWBHKOg{v;X@^d(|SZfsBc<KdEtX&qfYwXk1_Q-IKLv&z`Cgk8o>dLBXe&SMvl1
zlqgEI*&vnoQHcNYU39zR%$|B<%la3W=Y@rbCjz!?xoqZWXNdN8_a2>&*4D6_e?)Mw
zhPtU*&quo=EZIrgtGi>%16N5gyxY$6Z+H84Xu6f={lBf`JyfG!pKwq}*0JfZ+Lie=
zQklhHo4$3v4H4dU@@95+b|;_@kkpG#t~WOKYPRlQk;-#sfCultRIVah+E@UfVqcf2
z$!#jkd35Xenmr-PA=?dv#Rb=Hw=I1%6(;A{l&lr=cc3}!@8=GU@o|@H6T!o~hfk=b
zE6kjSXMGlSE94luFQree_};yH@7}&;W<;nh4GRr@Q&(44Z^dSuKYxDif(0JK%ka~(
z8hJAh{IBm)d%JWD<rqx|K8qI@KMgIGx;Gn7RY27XiY9&S&}X>R4+sC%oC;6tPS_<j
zv&cJKVZFNKef_lP=#>Ccs*TrRP@~7?_!g0nV66$M`8v`sQCa@N=Wl%A$#D^#J3mKs
zO%9JzQ|h0t;FLdIKUUOgPe_&L>9ZNd-NvGRlXqWLua0P&DB=B?Q?BSYS+G~Oq!kB}
z#yz(a&q>Me(0a}r@>5irj`;ra+Tsim9!f{;A5qHxpddirqxWN>_k7CZYFT@7tyRG<
z@+kiu`_tc?{rU9ut5<JFMMb5S^XQaYpE$8`#|~ovrRVTQiu(_%IYTo_PC`k^9Hh&w
zOA7Mw@i9{pOz~eQD*GOQvE8-vz{L1qwWO5Pu3kA1(3K>bS0^IfveW{f^O$vD-@ePz
z6|O%GD7$SH+1u$!p9RbQQ!SmE7?5>tQw3S{0D#?<=`Yt7Zy3VA7T6g?&*af7D<1z<
zAt@vCvDvG26BseaF~@hzQ|V8Wn}W~Z$L$=7y|1?+)7ZcubbMk$oZ;x`i1(=K-c%7>
z)=uBHHn}v%-|o)cyU7t2-i)6=e~wO$4Gg~$vfQ|RyTL>CppSqJ$%*-<1<ns%-(24D
z!uDsY{LIu@JDuD7iAwNAh1#)`Yj)<_H$F+9wQyx}a`MInE9BmP@}anij~<DE{y3mu
z9UGfuFxb+SE7P7mQw3AJJDFoRYE4s=OubjD?}i;ac2q{n9t4U9KeEq!YV@bsD^|pJ
zNOu!t&+gs3gSC7Qxc=&77zf^%+T+g_Egk=L107l9&T(z&1Jcx&()1>a;=2@&^JhYD
z=l=B97r1nYG4Ln`72%tfE?>?xN;j2QyLQj9W5?`1XD9IT@!?_UL$|2{u@!Qsbs$D!
z?i@1-H+^!Nncw1lkl^OupB=+*H^e3v7n?OcanJ|InJaGbyG1c~KXGVGyn1yu_~*XY
zE-x=H+0&nN8uOiYY+vqQloB7mhkQSP8qZ+m&{AyD`!O9&7?pdnQrH)77AzK3nJBo)
zG<i3Dah7S3J2O(wf&29>5BtQvefz>3eOyM>ZV`6Pq$g5EynEIR542k1kn*26wPrlO
z$j|!Uy<a>Wo{*ir6$tpIwbe4*I%ZgZZ$R0D(9k890|KH}FH()kc4(A4d-g2Tx*~k;
z!i7fe?zvayFaPQ@2>|(0=$0gR`V$>D+}nVgb9l**+7}c|92xBH1*yw-6!Ne3(NUQi
z&%v!f#Iur02c&sC50w{^*^{|@bhtMO*SpEQ$e?yAbD~bj^L~lc$EuiJVZp)iy&0@3
z!-u@d>x?oj(!W(k?aDkTwRyAg*|UX%T8HtC#tj+fvM1j!tE8iM>aW=0=-&5?vQw<c
z9vW+%vWVH`72fQq(p$Ny#c#}E=8PFVs7Pa5TjrNs+l~ti7a6%5M}pbQ5)u-M+OMq0
z^%}6Kysvb`*wQ!e`|oWlT;JSYL#QWur(3kpNuFn2|M`udL8tcI*jxgA=O)AZyLVnS
z?JFeBYWDj`#~z`a*FaB(Mnr5-QBkQ+H6BQ7ls3uUJ25(}kHQJNbLZHZJ<$SavowpR
zrhi)L(%QN@*%NR68Sbf3TeD`(uOB~jcJDqv`7?;N`J8PXx$RV-Bdbr&IA~s6cw>n#
zi_rC_Plt_zc&RXun(tuc4#|xhKm2GfGgn&ZM_;Cx=fgK0N$ifLEuc<=txa6LyvRX`
zNo?J!*HIqSpeCTy`Bj34ElY&m!Gj?|K|!QOA)pEi3wLS<wN)8YDHVko7blzmDUzEu
zg+XecXX+TvII!T!xj&BS<}9PJ-+lUiekCu;$~+n!@^?Bh_^Ajec~chcnSS4!E78}?
zG6o?Nh6G;{?mJd?h3PfWk|DB~a-O+H)$Z!eoAYrv+J1h1c;nkF)~;?{KIg*CrQ>=4
znjl_YUKjz<k+Cr;z~7O4M?+K7NFQI{;g^fGdjKZx!^xB#532L&xj&MkqN=EY2o`H>
zCmnQ&94ka+J#eOq<-hHHhg`AjONadC^?(EvvJ{TN;Y+SAiq|yP=E0Q(E8c;5manId
zzWDm;#&GSy&7Aw^`fFGNB~ngKPN?>sv6d{n6+L-YGg@BP#l<Dev3edQWQ7JSnjGz|
z&$3d#eEBjHh<p9kE!yzRnF1otp-JVn^QXJQ;mckhv-9%2dVEU+1>ZF^7zeX6%y)s(
zCGl#gx?g|(7{hRx<~ycex-{cedAXFB*j^OvHy?W{qW;aDbsj(3c0V2dc*L$gB}~k)
zs4)z-g6JL`5F_crWdR0db{&F;b4Oi-*3EZp4vme?eE<FsX=Hed%BUR#s~0U=v}DDK
z)XVdPBzNvK+2#4&dHB=c&EuLnd?;1>p`Mz^U0xZ&a;K9?p8_Ej0La0?!J*O7X?>0P
z_RlZOgmJq&K0b~g@8)LluwA8$GuU{#n4fEywaf?M))Ke}8Y;5k5hEn`SlKMlUgc6W
zt@y8BzlhfX%cS69G_<uF06JljktsydpmQU4INgVebq5|pzrU+r*C9Row(zs+YVFLD
z(Uj34&ci<~=d_`zYI|xP%6bpk6tY+lQmUv|;tim|pl!{Rw)!i?(aFg@#i^08%5xj@
z?2Ws+yDM<HY#Yi>LO?*k`sB$Js8YM8f-K=3j?x@<O#?JPjw$)^W1VZ)t`XH;|HR>l
zkoL`-cFAwTFYp8x`2{`_3m4^B#DLu)cH?AKw1lOlCAoBn?V_&mBiqmXI6B>+yJ(lM
zmmmo90L%>63?Dx^0I-7h&_?^qRKW>AKwYx#LkMLVMpf)CXvUV77Q5O<*IvAMfj`Uu
z_<DT*RCBrRs~O(^MYgn|fkE=V8-kGv9x24(JyaJ^8hpQvJSMRj-KL_Hb0MlgRmpqk
zgOH`qgP2{k*SpHfzO;szu!xA>MyFE46DJ;V8!`tjFnHMO9U7n1L!H2v@G-hAEPFmY
zP?MC9&}Dh{9SRJPT&rp~JF1d&2_s3_cK+PCwWx$R-G}=*In~N_AU~AV009p1{8m2}
zdeNif)(-AP0p$<S{Dx*`X|=U_3_3~c+FyWf20b|#D;+pD&a{BVGUj9$p-+wU^nwT1
zD;H!t3#AFi<H)6q(sVpMpZxA?OtSPF*$3PXO-f3_<7Srlj&TR*Db_YU>e(ewi2bL>
z#}r&UuPS(OAAu@Ua41BJfPoWGLmYato8IohH85IB&vIt8c=vAV>g=@n`tn-ZWuvlx
zMU4v}?1|lib|CP|U|f@~&;y2f^yM{3LYAO+4RNt3ZZjnre}&`CMII!zAl-n7q~~Ff
za&3^xWO8oq_N9V?T(BKFy*@-+nYl*oiu34F|AQY{UQT%qjYb1>@nGO!mSL(96AgJi
zFwnF8M1jtXH%t5T8Z7EA%~`xeSU82G0ptiB3%dl5^Z8d8;Fj?Oqoci<VYhFKGH@p~
zlM`be@NKaBD4@5Ejg4pT%TZ=2Co3a4PV`umh|iDR_-~*4JTG6lV)y-1mNyKZxsr$s
zP35;geL)Opuz{&*3e<`XUIu^^Jh)Kc-+%u_qv!z)+zj2Rf@mT-Iy*ZxT-rGtCe#Os
z51hABd0}PkAMUJEhr>4J1fJWubLaINH{9?d%kU!av!$Poe2entGLv$l3LTn@40n0;
zv-_J_dDazzo9H7?4#A=Q8X7V{zn(A&#}mKWV@b)R0->eWuYX-Q>kY?p*2;wer{z|3
zpX5vS1QOlGt;fg5BZk7??f^U3fsKIaT4wCN?jk9{3~7Av_o1n4Fia%1K`*`?F*7l_
z_xt<jhXkCUH@@WAM|{k*r0JG@&#5?tG-zPFZdI&q(p079f0UU_)DEX`7K=p$#D@<$
zI8|G)D``HcNb#7l@Zy4B|A-nt9vX5V=GGrvwae20twTUlNGJ(0k*sa)USK~?q7H?5
z&Nygp+(HI}p^N%$In%KKzfS_FT}emE>Ob6n%b?WH7h2AB$qd!JFFB`zwC*c1$NE{?
zlarIAba1;}IP%QWNmHU8@eojtq%Hf4ye@aF9g_M4S{i!BCvZT)BhU)v+SJrU@*P>H
zmZq~4>d3m5tX`dorvg0{p0T6y%#R+8Q&dzmFr&@V!0{wfA~lUSIBsQCU*us=k~@eC
zDC;|W)YNt=@U{hac6ZwYe96Hf=gB)ecoIjc4O1owWd*1LE<VCnfP#Co%0);5W2Pim
zb)G7G0a&t`nw)@m(1(-t?E8{{&LAIF52RP#X?8RwSA*0E9M${uPm9~|8_F)LxL6U=
z%Z5{%U{@1=8D29NT7pPx6u<|Vndolov5}FtphxK5?pe#J6gcWTu5NC5_&Fe#CIA}d
zzP8AtcL?^BgKqH}zADeJ^}Of>h=6EP$Z$Hrdzfp;DTK7JCuX@lV6_;7L^;d;Rb3T)
z_wGiTrp3j?l3hAplU5jLd1|Pu`v5Ees*pmYuzdM)e2Jd7x57`Jw&2k4aPHVZX#?^O
z$r#hFzP`R<ttZ;C%2nm_=FaWx=;*kvtGf@f$a&PFVPSyU=g*(%FIiRrAP6Q9XO)Vd
z87Y~dpoUdmP3>G=a3&;~o{r9k<#Zvdu<&hpxw^8lvRzAYA*ATTD5nw=k|IVfQFX{V
z{%Wq$k7s-!ORrf6a1e9n&SisFw>|kYJS?lUR^}fbJw9@SRna>k9Acn}lh&yRZln<y
z5`*V4TyLC8PfgC9J==NsF2(i7e@u32Eio*5R(CSzZF5ZNNy;n|P7_i!rS^;3;J3=U
z^3`<OYQMj`MZizd3xzJ0F3%u+LNePHddfH5BxgUcV6K2b!-5@6aYQDdOF#;-gW+yX
zW5WJ(Q24YS+8#lc=l-KV+-5)PdhFR}it{Gs{;bg(2rR=)3#G3QudG$g&{<87XH5<z
zsSv0v;v`17S8U2pwk-8KwUQ%0A*4)m`2Jg?tZcFER23VU1o<H8Op>EM(pOj1+PuHO
zzCq^c#P3G$_RB&`mMrP~%H*NhCyqQhT?3ua<t&fqSg)<EO$-yF7q4NPf&%~&$z8jS
z1iSzJaV#A|Q7($qRaKtKg$oxt2L`Ggym*%BA3mH6Hr3TnJfwSZ)&c`i8!>R&$|H{^
zM;l6^!VDpp?*Z8o($b`HW5ka_O;L6QzVz5{CNvImRzP5=SVCfA3X;HJt>!I5L;TsE
ztQY}ndy2Pn?r?z0`0hD#=De$_^33p5hBfc{_U&D}GjU$+?d^8(zC@G*g<XZDp#ljW
zz&R4<LIk`RBU08rTzT@Z6-s}a5kP3@_SbaK>dY061-pT^UO(BCUA9`lLMH5{SF=Zi
zX`WpIxHc8KjL0TrePY)<uGoow<SBe%5}xqumn(UfPgz#)-o2~MFlU9rc)|Q3-?N9`
zG=&JxDB?;Iog;Ye+r3jH_@3ysoYI337LQTT7z1AQ&?q8+-^6y9By)OVxN0!BP-e!n
zgirZSu@-(K)}+m3WerVCOcbD$PRm6RuU&lB?bXYd$@)oJWU^4?*_VPC)Hvu?tExX+
z2k}^LuEETiGYR1zKYsl04C<M6J>gsX1}#M3PsGUk_4oG+E1r4$n+EJSm5l6%fUW8g
zRi+Fgph{Nlblr>I{KE3egx3nYcaJmvr@bryx$i?FG+>oHhYS8Lh`wN`hKvxn4M(6q
zox;eYb8rL;Rw=MYy6arBtjZUXEaGg|2&ep-09F4AK#P1HW;nzBwK{AYjR5e`_wT0>
z`btPjzPqd+pCZr8Oa0XZJm49NzD3%z^PMd7Z0q(TFT8&9rU7z{$$_cK>a#nfsbm>A
zM1A}T!RWe&f@@SMs+5Cb9<7U=d_=4c&JI23D!=mA;9%X)j&cvo5mHW7-v`JT<NhOQ
zIfr*OCRD`BZ0S?(+wItI&i?D2q!q({pvI3Q-uT~tA71HyJx-FTR**9?fyDUuwJ1Gq
z+w+U>5OUkurDW&i5Hf&Pf+ti(NpNc(283_6>?flC=&M(+B9+cPA+M=)?zas>taoqT
zd_b-Yxgp7FrZeCp^BkK^S1Gvb;EUa>mHrfXpEls(;iU>}@L_W8cy@C3C{0_4u(|W+
zv!RO3ltlbZi7HeN5~;`S5YuCQ>Xg*Fb!z8+f7eAY4>y7?f7ID|!Z_%<q$tkD@T5uj
zV84;siG^*hyIPxvq$_AW*i!h7(3qHXFaZ+`cV)xNSJOYNr>pBybwd@8oQk4=LhKo=
zR%t-JBh0Zqvxmyl1)F}%Fi(Z|A~ua&(6wt2z_-klgk=`t4U$1yTS+wmHV}a|z>DW0
zE+AqPNjGlAm6l)mrIap(qv_|m@UWd#CWb;t?1aDo^3kA{dBmsu<x86%ZO?^O{E8u>
zV{XL8pCpgh4e7syxOozP!_xSTbO`}I6k?q|#GTiQmVW@xWRK}b(I_`{>Egv5FRt;h
zXW4&D(+v_)BF8c{5<8_^M<;6L(Nlg7pR@N7s&RD8$SC|K><czPjiSm6s2Q`VwqoVc
zo^`NV)5AQNN`MKOJPdQ{)cd~V)qVe*Gt5bfwL@d$s@#SPD7*7OI`RPs4wBp<h?F#+
zhXZ{ep`X}OO;6@b{w-Ir`;trB%0vqaT@U-fLv^FkBslD|p9uC$j(4k6G&h?=8zPTr
zfE8S_Y*~-uM7xjBEnZbaX{0`=>qvQ*c*y_F71;*5y6rwKm6e*rb6Hs2S04YNHqo`2
zKBjqF_)<t0j%)b+wrU#)KyGcq4a`$!E)Y>X%<`oBf>~qC>c=3ABnT5I=Oe*a1*nMc
zNwIzeATyr+Htn`TWLZUP#XhUcv{WI&4o=p+KDl7_^w|^_7jH)#^K`Ikm+(&4M<h&v
zQ0Bs1*PzS^C~Czh#0)Mpp`0I(uR%nWigY8mYW!%&5b_^{3dpijlat^4y;*PuzcNlk
z!K*)QYNB<jzQK}CCut;3Rv!w=W2i?L$X4jBNaY>CfK9=@=SNf{#t~B6w&{ZgB7D~y
ztVVpTBM$fg5P}?$S=>>c8h3*Rd(-Z0^X>h;u>PAfDIqDa8bF8sQsDBCXx#tD|F=6w
zXY~Om!VodbIeaNNs%IwMe_<v7DQWYh;#1zD%9A11Cgekj+OZkcDpO{79&c3P@SAl(
z1!qht+0Qrt_d3lT9M~QNO*&viVqy?M5yJxtgUd0(Ov`U_<do7l{iBS;1gYh(XST^4
z2FkcU>D>jGN3a?%V_wwbRXP^4|CTTjj6#R_<K(X0yvarembrL?)<!9*gDkHBBS=b^
z-CUm$*Hud12EpuHTYd5`_$HbA%Y?4419lCPtb%mv<H!Fs`wUiAxuv5EQV=Q8<T@dA
ziH1FPY>oYSD%l7~CMhj#0APTc^k{p=LvnTiDICE^z}|4V%4AKw)mL7Eif%|wBT0hA
z2!@VDma^TcMUmt_mh?m(T+7$5Ux(QG+Du3(nfrpO4nL3=F%L$_6d|6eW|_OtE?H=q
zu8xj%NC?#G9w?AsnnUtKGB$?HwL{YK<?*qLgbY0fTFh0(f8a!g#SmpB9E{+?o@;1o
zb~zrVWUe8)b8`!UHzFlSOyk0UvceA>T1r|v6ro+agbTA^o*>nhW|Uzjfms2G{Ak$W
z`ny}Krbh#RJAy`klnIc?*#;@fD15mLGD)PHjrHsJlH1WB;vDv<Na<Mf0CFA#m=?B6
zhw9{<TGC)_?~43e?~fsJbbYpU#gfH~AKt$G4~*G?U<l;U+qW^?oW-<w>YZX)T0#b}
zj?G09J9pBbUoGZK=8y^mpiFeEo+560mViiFGs~+2(ZtQ!CJ++He>>wDE)p1AJ~JWu
zEJndY@2U52Co)N&G$^KD!^3c)&MgsBNSJ^tOsVpjI`h&XLDNn1YKi_NwTg22i2qIy
zkLaKV<At=CxHuDGIYFTt<lzlL9K=MNc)N(02($=M=x!dzS<N!|yCOUHRJDjlf79nY
zKEC9iXpu)C#S|zClGm$D{uCgxS3qU_F(7yi;t>Wa%gaynLS%pkK>O52X=ygvOV;7b
z_R70k)GBF&ZF~1#S}3BJjIZPLW|kUZ!GhE<NLCUC!@S7dWm+nNJ(IzJK*I=Y6C@3H
zQ8&Z5v0kz}k(yZ|iS#i>-enKSuLc}tuC!H#7ZfCa5ozq!ty@Fy+<CC+$P*g~cM>wf
z#18kf2Dmh4Hilp<cDD6F60ag-f<DO*S)sZh9H}749x~uS@SsgODAygm@bvc=dt%65
z-&#QoL=SRgpC7tv|DtTSoo+2L2W1of_a5qL_bG-RAQNolI8Pj#Q__#*6N=z|&Pke=
zG*f5A=`Xnq8hXL*W12}q`P1U!p|3YKBR|$8S&$(du?Lq;rbT30awG9o!cmOppDlMR
z;;f%JYZe!tZlcm*Q-$8KV-LUqPDc?$)W`3;F@!hb3^!Ot1Ox;Ox>RH%A|UxDpjdLt
z781u2MdA7VQ~Qf+E?EWi>{XQW6`%{)G)trL%W1i<e9y#LkX=LI3xn$d9Cr0xkod^R
z2nPX<?w&mtkQ!v7cnfXx%p`%e&Q|~<M4th2@iHmN$^Y=b=P}cQqE_h`l45}6Q;ak3
zi%CdJADRN`v2#vo6Q!G|zE`$=2`x>jibLIEgh?t7!CuvwmTYNcH74QWsJE!fZ^eg~
zP%<A06QCt}nMwCuXO(m6GCc=cHlc|jpGN-vC0Gc$4bIHWypi;JTE2vhW75ZMEiHvp
zAz>k!mx)H<e#&|X16av1qEec(tmtxr)o;PHpnq3o%ndU(+&{C2ByCz+ff*&EcAEDU
zP0&NG?YRX8O_<5*y#3Wx(be0VN*Wk}G=672ae}Y|XV07=@gouoLrk-fO~qjj6L{nY
zZ-E$Dw;{Crw~epc3X~D4cs84i(rKdoAihU`H#(7_HIgv6-15tdJM;(=YCuP3K-}5?
z>U>M`Jm{&=d-s?KUAk3Zy$KH?2;T14B(E8{{cXJ!U-I*C8S@lO`0|`uTSh*iN=ccg
z0Decm=Tx)POtvBC(5zM*ya`Bobe@mqFLFL^HkWWF17&(;;VMZuz{{5|ZA3cLVCHsf
zaAbVk2jcVqF$6#Y;38&_kWxrNNhC<>lyEDgho$du7!h7Lji$3MqXcxY4J9QZagAa=
z2Xgs`V7!4-{o2tXrL_L-*{s6YSUMIVA_=>ZiLz|;sJv^y?}e+Cr7_b(T1SF8mC3&k
z@vZut@HWIjge{3$d|?YTU$B{aW{aZ+C}g@z9bz^27dO`1CJkKjP|%{<^qO4CSCOo|
zWc=4b61Nb83c(I<NQQs%`M_7&5RwlmX(#!8#CbRq-v%o=5)~`s+G|eLlPL_Db8SEV
zayH3!z<{K|`w`S=G;T7@J%EQoTi?2|exC=?rCRWQ`uK;WA6#9Ey@rQ}hiu@;TD_i+
z;RPqsKBI<DKv9&#wvF>-(dJAha})uM?9peGkoJv(`}`>?8iT<&gwf_eBdD18_!nn$
z%N)DY%&7=%gptHl8+m^pn#l<!15H@41*?^tmUcW41`1J0Xh%ZowamJ-rsnE#dnAP~
z3Kf2GEay{g#~nUB-i<~@5riR#@bRU65v)lt!u(#(U!`CG<te1%WyfQ*ho3qSOhPy?
z)f>lrqSuJ`!kx33`Wpx?VNgvLLsH9U=v3sQ`ym<#ERtr!)C63H8cIe`L54LT?9TRf
zmyZJbbOuq?6r(o9FqPnd(l7eWI*7R2#Ad>QG&Px_2^0nbrz1GXB~pQ8WIAo7FhC$Q
zCkCj}g*Vh-Fls<0P{;|07RY3I%a$!9h6c8K;J_EJm0bw}xGhE?8zvw#I!SG%;u#%&
z%*b=+j-ueuBeXKsg;Ft=AfZrd`7E7_OR*;f2(4GaI-aZ@`#_N5ud9ePLm>!Bbc$^e
zlBz+Oz%pUlzZtO~44c_nw2&xcbYg(UpKbR2^>#yd@vHu~Asal9oZ#M>eiNf%MhE!)
zZ|Lgj#V?gI>PgRUsZFFcqlkE(7^$W`Gt^0Rf}JET1`S9Y>X3o=e!b?xB)ucXgt77-
z<s=QGk2MtXKDBY^oPK`)iuAO!Ls&;DS+9LR4Qp#kqdgBvxJ@SC#z9x?UMnr-p@JWs
z>n~I%na<=$Lw>7Lqht@IvkDQ232kzZ3`kWbINC&Xk&GTQB#m|ZZ+-|CKWJcK(a+tY
zduzHe$Ye7hpnPYnm;B0=3oz%8Ry?bR?2B|HS=0kUn-ez#UV4OCdy%)I!D7cX6!if{
zY`8{!HAx(+>I1ZQA(xmtZyu73B!7RJVoXitN>c>}Z{NNp$r&-qY_<WoQVZIK&_brg
zQ+E219a?26*5n`rWz3S7Z#&ro!xcU>WxUz~(b9dtQQL`M6`K(w3ClajS<n+62;vIQ
z9H(TP@p_~hQn6!1Y7y#+Y%VoAvNXw@n=UOS<}+h@LVVqBKe}%SdG!$t0qPJ!uHmq`
zWtgVcuM1L)$Io-I`Bw)sMpW)Tb#(&E5ZS7l`L10W#OI=6QG~(@?#ZKn1_sQ8X5<N>
z1V11%CQ?NTal*(aG9^L~Nao=m%xSM?aYhlAymFn}_mK1p1dt&tq`C>4Qi%|d;5M7G
zD+`E=!GbACelWn0T)%!l8Of4d``>>btk;Zm9zC_``Mw@(06Bay*mm;$hfVoSYdLIL
z5i-n#z7pJQA-#I_YDR5Sp&J<p5|M_uKVv1m9AAa}4xk~ve!Zl)_&(qf34VwoLj;X(
zz(i0V^iQH;cP#9&yiXYLk^R+GUZ`x0KsG{F<TttAfZCrN&ZIdIn8d4ZQx<!RaE`2*
zVZ_@zYnF~wXIIyB_SRq6Ez;~(tyEPmMN3NNLXG0hC1K%k<EDkK2}rHX<Ts@(R#&Cs
z0^Q_SdVkL<pI-bU=^0p*tYk21iOOEI=xY-};yPr(Oq7jr;QY8N71h-Z`0YurZV1_9
zW`Riec3Hc;Rr3SD3rW7wNEs^x>jN*(=82;Na~k`KFfN_eTb$Z`sN$!A@J6TZpWP2%
zR5kM%iQet%C9&P*@xQx<(;{0jNNS?VSTvw)4Q7C+e4;0P#xVu2dv=}=gxLAcMYmEO
zdVwonkxl9sySLhrDs%}MO<<!8iYc+yGnV%Gb&ZH<!~z9`{$7#pfu<pgN43md12LMX
zPp5$o(_-mK6U^4O?_4Xh4p76Dhqi_X`tig0mWWbOty{P?5itN7i1#`qU^{s+vOtkA
zAetg#V_QvPReZ(p*ItN?F9;Xk{&c<vx}YhrZ1+e*Fj`#(oZG#b2Q3}>-?0l)>N9k1
zgocF?;KIUG2x<NN!Qa1)z@XqqTW19`M>9TN%EJs>Q(`+i)PkFuwSAR|h=Ek^0e#hD
ztMQvz@;{V}5q{*$6Nl89^91EPU{j(!8c}e>rP7cV|H-eB*)d&F-Pq?XZuNtrCF%qt
z0R&nEal`u0ctEBVI^%@SWBZFt0ySfm%}LaaHDF0!>R0c>EY!UYR2&uuD!?*20c6XV
zggOX^h$zBlTX*v*l9rNeA6A0AC5_hfsv6AWY1=dO{j>g#mNO5@tPqoLyJ!A9ST@ta
zqzH>nfvZ(Au;F5eYbZp&b#igLhxu_2kFkOPlk^r%(qM5e9n+RX5hd?8h{ND4QB)Ib
z49K5imepN+Z~<j@lbj8}u%66EAw=jDP%*(Jyt8M26O&j9*iSK4yEW4%;z=g1ums4~
zNK*k-EIA{%QhRLKU}eZoXz9`f0p-8b>$7lO)4&3SVuuXO`%i!R1GNC|$0JT2!*>hj
z6A<_Kgi+cs+(GZepH^d6*KB(9cV8;z^+vv<LnO`v$mAV;&fDx<wje`fsp=+#dJiD~
zN!X5N#_%OWWKP5!V44opIu7CyzElSaha^83oskrT;6Y7IO@;_2NLI)jjoOmU0o=)D
z(iPaQsoGU_8Z*$v>koYx9~m^HF{J>n+y4Ce3X9gd`}WODI?hAMEXOw_4>3CM{6koW
zCb7+zB9BzE^Mor>gNyN5I|myAB!EhLa{6;Th}OQP_$+b@HspCdCjZ7kGpXeE!*f<a
zv2&^uQ1hA9%43Ia2yHgP6ptXOkoGK!PF_G%WL$7v(&3mA<(IfOuO^oOXeKko<y4$p
z8#dG1(h$f1-S7|((Vp!34lG5_Ax*e3psZc-or?@{$xj@ju!J*TYZjGUNnZCMs^+xZ
zkpz}0oD4b@8{G^&B;N2L37Gz}><&??o62ZeR4PBw<Tyl(ij0vAVYW#MXWNBN?7xw;
z*wEO#3!0nCo0H!@(_bMGp>Ih^Nw64Zb9tvf>)p9?hitm6UAq?hGKp<%NAqlI_K>L#
ziBO0+JN>!#RR|wt$48zy4*5L61xVDo&875sJnB>eyJfXZ{T#xCh={NYLAnm<mGjS+
z_+Je!VEJ6(Bu*_``XG)zOk^__s9^ia(h3ojgdmZ*5aEZ7G$GA-RI&p(0xL+jcz=J2
zls{G;@oBi(u#Abnd5Eu)k%wUaat6gVBDRWq7l%MRB9tA8aagQ#39tC2BqfOkDPAo2
zlZQ&y`ri*CTE(n?Hx}+(*ku^bRaRE+3*n)-ywe>M3>}Gz#SkreI5f4@^e!~R!B_v{
zanyoCkllro@CA@d#wXkl^n;6Jj|k(y*g*wszo*PY<$1ihzLbzv0t(y)TalyRqMWdF
zAjTlU3K<ZW_@46sS3m=kwPs{XH<D&i?Df-+cN=)g1dDr*bVZO95vq%1jbtGaqLHL0
z02jDWE3z^Am&S%8q-l}&X;j0D>8B(xaBkqKF}|40oG`f_@PW-$#f{X#3X<>!^~QO<
z!pXfo7q@d237<|tL17^o&=YGic)=~I@^#*7in5yS`MqdWGRwjk|CG<CzEdspw`|_5
zjj;j7O7uOvkNJNRDD@7MI-4gc$X$MAphnG|s3<80iROvk0^pE$fGfTY*MKy|Om-_3
z7yfnn6Ug@QJOHVBv?uJW0ue$aE5xCCU_}XD9;_=NBct{IQS(&ZKhx|^7eM=L_v$wz
z;mZcCX!84shd^RV&HxktA_6O=Q)8uM4JoyJCdJlU$%A?%OD8x&BkZAw(Bm!`(rFT#
zt~>=Dn2LmvOmf>paX$U;Q?Z4e7ieul03S>0c^<t+o6HJy;Aj=d{*f3Nm?49}@S;5s
z8EnOMQt!`B=5U_#uV$0^d6Qc;2MR+M5wxuJyCo1s<TyaK=vXgo*7=(T+7qpUm;4CX
zJOpRYX`PzP2aK>2HCK@70$It!l0b}CB`+0M<T3-DxC};tCPd%rH<nJ8aL723>>xp-
zeyfhv_Q7th)O1ZF-6tA@kt8VML#!lLcgT5ZEB?Gh-W?*b;EV>P>adr{GS5*PGFc;W
zFOp(Ql22e>;f|383%R-KCj1Z;4FN8TQ*RbJ+Z}`PkL$2JUg9bOK#Dp#U+4i}NS3Il
d(6=d`tj|AA?duxfh(D(o`?U6E>^XYj{{W?!{=xtN

delta 13612
zcmZv@c{G*l8#n%->ExhU$&k`yJe^V^L!-7bPZ4QSnP(zA8dZuibuyP}vkhg)OeZ8c
zRM?bxNJ2<R=HKV;`~KE?|9H<@XRTAm-uv12eO;gHGhEWMIa0GZ9(Z#IJg1Y^&=f_j
zJ?KNJ$et+n%v;1O$Xp|3{N_StU$b4wY<2qf?b|<OzEfg7d-j*RhllvSeI|8@nr=#F
z1_lWwGd(hv?BtNb!ovDuuPag~o^|}~%`eM0tU%T8i@3IsVt0J546<x5F1p1hlBO_R
zn5kyA^Zc7-Lh{Z>Md+F5PoF+@?s!Ra>=m4!-}tgOnm#o@)2lq*=_``#GSDzJQX%3p
zJ(y!P-p}s;_`JJk%D{3}n;XS+w+=NoN)zU!ELVnCC*(QxTE36j-`jHG-)jr=2hs|H
zqoUG$e0(x(R|*QIT#Rx*!lyj*Z1g!kUbn5SZECDelXc~rK){617MrFlMx+Am7Mi84
zv(+#_x59Ia1#_303D?<L@5-qdkBqCj&#0FRi>HSaTn3x;brX(jA3b{e?Afzp^}3o#
zOy-+Bbq9y%`1aynkvTaz25BZ19xcY#Cxj@IfZwmIt*y;+>``-fcekjC*!#A!^20wn
z&N+{L)$Sc%vZw8GY=TB)@&KJyppj~vUHkN~->!=_g8F)TxtFKH{C>a7rVUgFC$I0p
z-QnudUOykF<@bxXdX83y{cV2tLUmZr%&xiAoHiY2=PctKTPZ$1zWeKT9x1)W6BHgE
z?w6RJce$P8<@l0il>OD1m<<aTFJ_%7y;W1u)!p4=FBcpVvYBC78+XT~c=qSX$bIfQ
z>(;GPV=%O@FXCCXe0gW`bE|zgZF@z|1(Zq5m5HII61rsR(xnjr@o{l)Y;z>U#NujS
zoH})CDKD>qyL;Z3{~pIS+ZI%h<z_EJ0ptpEa8m4;pd~^Iw5Xk(9iukiWq^FYUKj2?
z?Vg$9;^Mn(TjaXi+e7Z+t6B=<W0L=sDVcfdFQg#fcxSDGYl^(<Hzxn7%liXSB^<pn
z<$<3E23R+_`Nd^qG)?oJKX^^FCSK)Q+&`Gp?rN)RVv^Dubwy)jWJF-sr3crpT_Z;?
zYgzwB^y`G{=5DXWK4HEtDkyg8WdlCBEAa1v0U!L4V%`Z76%nLgx3amYm+^1^`%h0p
z!=F{bbW&4O``p*p*T17b;^j+~o40P|d;Ce|Udk3%36(KRYc!)x*YB~d#}n*w8#3V$
z*!8v0Ojk$eDQOE*W-i}5_Agwxu<s`SsrR>+ZxUlTefjT0``phxgP$^0&CSj4`ujHo
zX(}a650|vRbla}y^u4WA|L>hUSFc~6*wuBwFhFPGqnN`cx9-BiqToV_$4fZ0&asoJ
zy+L05n`z(DNlz6OpD3@NM&;$@NvDcDrsk%HWj%k{%UU%YpPim4F=<!3UwVblu(Gmp
zQq$8@e0vSMyIyoYG^ICQL2O~ZN2t7*s(6%cN<@Sp=K_`c`Yg`#FMT)!{yMD^l4JJl
z^wm&l)8w%+=Y<Ov91VO`S}I{vSKv0RMtgsIUa`xHZjekp@%)<EHTvd+a2o&CHC#Kl
z_tlqhmLE7x`Ectb${d&<Dj2TMwlY+I{P$r|QPI11?tIYIRE$T@XWH5sg#u&ZK7T$x
z{_RuU_pEI0VMU5kY1L!-`<O7ZvUa7K<msTa!ggNtcXV{zb>Y9YLW&-G1bA%SW_~#=
zTud)!#l;=6wY3e~Z6gr<(rjbf6^io7PtsU?yR`E`^@X~N(OwC+c-E^0CLDkKRP^Aj
zuSFAWL1<F3r2W?ysU1JAw}1at@(j<<Rfzxn_m6wJH?&*|Cr}M8x{<sO{?{%FG)HLM
zuV0t<3i9H*w#&;K2cFMqUcGwtW(M8={>k;gz`$J=wPN0X6fgMyUs}8<?*G{;pCimJ
z8in5V^HaZf!fK?Ai=yh)QAg;f;D7#MmX#f1=y-VKuUNS<<niO}jIh7XT$}sut1OpO
zGSPO0o$GLpM?}d)!ELDiM^Du_O3R8ayS=<d*XNYixW1H0u9k&`MN_ubBZ4QIN|T2U
z`tr)o&6-wT9_FAZt<~B6t5VFoCl7l~jXBJ@fBX16FmAE%WtE*5Ye?}%lq@?In4X@#
zM%sj_D&nPw&eHCFeP0$pUG(FvM2(F`)V`!0=U#7<kT?bawQSALR+>FArrYY$X!@o2
zxq#1vFpf?COs1rti#136iSXSvI=wZKn{8GN`1~{WIPOX9w&%RJ8Pl6*)3Y49kGS-g
z(ZgwF?}QDzdIOGqdz8YLymh@Ld;1B}ifPtCox*~&u1!v7_3u8^a48Z@X5y-&fa=1r
z8cYBDQ}5CxCaVHJ%+yYnQaO9BThGMggsF^gat5R8$Xn^op27Ifrj>`5P=Idc3u^Cb
zk!ZEMth7gEc2ha$6=hmQx4&m+*6-M>GZJyO=kw=f8}L^AF=2)W0#^Es{1fclZ}#Wa
zv!diX)h{u*!JBzibY^Mhvkm`pDt4bLl6j<8B;(uuz`NgUZ(8J2jV{F!&XLgyN$<Jb
z1;Z1oZ$~SB!5=p#kIcR>3N_bnufhdTFY_Fht_kAUC`sp@->MN|KfZ+W5x7QiP$n5)
zzg|i`nN!0)M==iy@~Vaimxd(bQGEYEuZ31$&^>hMMoUYJi6XC$y!56`2lMjsEL(D)
zvs5{>d_NAlk^u&Q`tjq3PPCHO=kD(A{^wc?1v!$={QUMQ_TeVgehgIuBco47IT%-e
z?W0%w<Xg{9jaxQlCGAuJ1fMqw^-JW`3O5RkpB(MALjO<A4%6wV3h|8_51?k^M<?qv
zBIKUnKg;!FyCj!95ctyb$jkt9Z>V_M7LzT~C48_sXCSXPs;;xtmq&APjI?oycM<yh
z^YE|?!_mnpsQQAF5?Wr6m-dyc+IXDO^FzTGjny8z8xYVq(N3q!6#qK#<jE5OMUSjC
zk|(tKF1~+tlUo*D&#&P(|9blg**AEivF-rn{@4H|o!BGyYw}E+U)adncc>gbeAwdS
zvoo*1aB1z?ym|9>6scX|$gKkh53a`e@8~$aaM7Z-jg5x8!&7TzXfDb=iHvicqeVlC
z(NZ3s5HZs$gW37R#nWS7V*@NFfApXy1D*R11_o`_%y%@I;pdE7G>NZ_&S%FTdo*fm
zXl$H3H#==L(pA1-iBFX+YQ3qzE&cWV^;;R1AD=xjeqo*T=FQP;OSanb<;$0HbL&=x
zO08bImf6ukpT4%><>pV%&Rmyw8~S$d&#H|24D(|};|(VPyFRZ2MAAw3_Vn~Ly|4-Q
z`1!4`>&}SPf}=bcHqE;aA3F4=u5PrjIj^9A%<-mTuYxC)RJ}g%h<v_5`l(brgVpGl
zYLJ1*-8M~em;_%7hKqy5K3VV1zVzV$4=u3keQo?Pc9w-Mzb0?Y7720jWan<bIuT3f
zZy&Yk*#Wch_Xqis==r()>7jybuW45v;e8+f&%KvrNQQ)jB*nz2%+1YZNDib#Z<H70
z1sT}3Z(qyI1FApRruju>Rn;-6lP~n_?Vml3jU^b;JGxafs=g*tA^S`zy(rbJ<P>gh
zlM-qljIugf*-WPJ`>S>D6n}MGpPKBACbv$W6*tu1l+EH!D9f{7?!i-=Zl;Xuwdkvj
z51gY1ql=%7j+s7p9coDm4qj#I*6NZbuv$gvxyQJrtX*4+yQwi4$=j-`NViwXW+lb>
zu7k-yuuC*ZN>i>wtkvR0THUDm_W8Mq{+7JVZ?4lucD5rwo9!evY*0Z(mZ%*&#!gUQ
z&wu=`({rPiz_{1jq@>t9ds<?Ou1r{>x$>v_S>`#>pB@~%^{2M6aTwmXyYMe=O4@g%
zyFvn;YUVj|(|UKf)ADC#W|!v%3Qm6-{FzD5wQ%}yYinyiDJ<MKIXOwj_Qlo<m4Tc4
z8`Gmn;oQ6TXX(l=Pqf(u(BAy{x#<u9ALbWRH9&ajAAdab_4P%8js?&&2T@ja5qq7S
zn{9IPJWQzUmY;*o`st?>%U{3F_gF<S%R249hT9cuW72_t5A{2IrzKEKS=OKaR-B(1
zHO{sSk*o?VC@8=|ZDu@t_)tLEyHMEUt0=$bA2DeRO-)pKn5U;FsZsEl%IfNBcaH^b
z#hhBEH7-KmXcm)+!IM;2RK#=qX7WZA%V(_xZdQX$*-bfi(LY8<R|^ZL5k}~KOrFYZ
zk*un1qO+~WL7gpz+lu79XL1F3eXm`=&Mvy*#<PC+ah$B{_&=17#1amjgXhnmXX0@1
zE<yOlEv>D40|b1kywHgGAbG2nEsMbzuWWBO1N=4RyRg8t2|3&vb@>3m<S9#ap%%BC
zAn#ynfguM62dM$P$2;0x>^zQD0NfHdW8bzg;*4-kXR4@e!QgHb+;)f?9D`+D!hIAk
z`{`kBS?iBlm^?LqaZoGQ$ec=fbjpJ<(c%=)<~3pJe&tFP{+XP?<Q<<GeP`8Gb_WBV
z<$UE7mAe*Guy(xBEFElCTvGBlkXqZ?noa?cV#hI7%3gZ^(mkI|y27GvOIP>OvgHp1
zVxB$IDlRS#jB_1qO3KR8!&eiRk_ss)DS7qkReipT)n+lVZAwa~z@|{NJ9|%Z<)$3*
z7fiaon6ID60G2_jSW<0k8p_UJSTPv~#E@St>3Nj1L+5z-ZuNkDKj_;a;R!P{p5_-W
zB%_gqg@rBK3iE)~!4Dq9eEfKVa0&3mrswB;SrxD8PP=OYF@=RDv8n+LEqPAtXJ-zQ
zK1Bi60+54)g9(MRX)iYAUL*M)OzLxQZz3w&Fn}i}^kkl6Jf5!a9v>q_gC|?^818`<
zQ1Z;(q8=89`(~cZu}xdMT1D#NqeqDtgX$4-+CmDhdZ0nP!uy!@^?InJV+F2*wT14Z
z!9hW>z+=mLrgoa}x{LdL>M3t;Z;x@dAv;l11k~<I*e;8}BLTvxA7WLjT3by~uWzfX
z^-rCO%E-vj)6@Ih|B$1egT}PzHD%9?X@5*JEf7~$HpAs69O2SB=6?C|ww*h5tHWe<
z!F>qhVkG((hXw~TaZERF-UMT;pPgt=e__*{fB9?#kKvBSWq1OCeV3mpO@G_j>RhKT
zu*Wt6(n{9p%kjynDRSwUi;Yq_El>ZwzPR9bgVJ|aZe1lLlnP)z8tBwBnz>zEycSPv
zB!RIo00=m8Um*VbcYB88T+8$3!5=?=8U$@oU$tgU0{)VK@3{(p>4-JG2Yl+6g_<Wl
zeR@F3IO~0Vl5R*~;8S5mkFSp}MtvZUd4K&L`9@9-pQ=<;LH~zX-{N0g+=DH2o-Qc!
z#iF4UTxz@bbn3=UoAgV(z4(tWjoEF{{33OWMg(&akHAs+iKzp=yh1t*#vN4mMaO-5
z%h8k{tx1<K<i{E3w6q>5{p#SlFCeQ%geeA3BFc&YhjiZT&*mQV;(EJ6Z4uA!SAPz*
zvbws?hJ=Qy(ef^8WdL2b+0p3v^A36x`&Daek(yg}e``Sw*mT66k$*#O7Mx};!;u?5
zKYKGHGc((+U1|15l?=h{eTs_YamNRn4e9ym_Q;2v)RzA+LLhnTz;<vDTeEGN4alK@
zi(mqalNy>IRi2Mc(oM`*LGSQ6_M_pX-S(|pRab7+I4_kY;7{ef+Pr0pCHS<sgaiQ_
zs6NmlU|70UdXWKW31Le_GxzrP_B(LYFCS<~)y2Q08-(QSRYbhT9(vD=sK&>~C!vEe
z?ewgzQ=&Y7FnHGeeWNZNldi4k3Jb&0gYJ!ViwJ<&S5pt)tq<AGy|wh_MwK-C^rD|(
zTaNz|^Zvb-sS8^_A4lj?C6%*mJ!O)BTLyZ3fZR;Iabp34+{w@3_H=Ok%1--j;^Jx`
z<u#*O9347D`O(wO?PDF+`QCQg6Oq1ammX?{4eB)1KVnaF<nj!!9(s9|D=@CMR>$V!
zv;FuYb#aWK2)E0(Z{OC^(+dU|M6jN6z3j;X>Dat=tKnU4O-)U(k8K+_zPcI2!Au2g
z`wQ-Unwy)e=GIbGS*dMe5`iyLzCk7Jbf|e>^`3TMykb16-wN;El;@P_I^32DR&a^_
zlwl6>{wOgq5&gQcDr1*bgXEeiPJ1K53wG?-Q6sCOG5>dmz1;cCPjY+Kt7JC=iE!&g
zbA?xLf!zbJONJZ(Ob2SrO3F8yKDRd)#Rn(OYpz<c;^~72%SbQQ^uB%jHtXU$VFFIz
z8+B;+Ebm!*YE*jZB8q(m9Qf|-+xBz!<Y(<{uW@uZekyd&#DvHOfq>2jGgiZN)g{d*
zC@9(2lJmJd#Fxsu7N3~NGE6@OJ1Q2f69Ul8&6O)JD`PT(wgCI#K}2|?Fbh2ND3_yZ
zYIjk;<FD7-ky`?4Tha|Fj|F5c<EF+3s_N>3VmDC9>?p-6&ObggTHV^ceoRalqI7V(
zGB|Rt+3`lg9?3(1J)$k`8;c|a`Jb5myW9F>3{h~kNlt>tU{ZYfl7RosbRYdvGoaj*
zWf2R}il=&m@l>wM#T?}1;_B(_e9}dCofT5{HX*Q&n#MPP)xc^<J;~B%8Gu}9gITW<
z5P17xIrT^k;7Gt12S*fYb@(2u_F^w$j(}N#t)xRDlaIg%4p990Prh?sSO0Sd-UV+~
zr=eh{Cx!`m(9zM!vTBq@XIS9H@{aW9+uvLZ$kovE-uHzSjfaaSqFAbnwEVz%w}W~e
zdnB#zw?+!yR6OAlgj%Bjz_+^h$uz%vr-|_h-m9?Z4W%`B_wL=`@NgFSIovRj^g23$
z_%eyn-ct;8OA$-8Y_I*&hSfFI)#`+kR)x}1Y+Rq(yVO%f#=40bY&e4(7!Uz2<+rK+
z2M{Hb)6){7qWb}#251G3iJ{cHcbA3)2XDh^?A^Qf?VC4R_>kJVx;uQ-qd3gyc6z>O
z(VuIz`s&4Da`n1o5RI*LN=s{lh=_>cy}RO|MLu;C`L1$;Z)<Ak2GD_B7!W|5{$C^A
zq+~oj3x{&Mc?|=WVrGJp!MLjVRn>}TR#H;>F*2eJLQkkD@%@sSO#2UUl%?Nl34^a!
z7Yl1)SngE8LFj2|#bB;RJTHQf?_2J06$Ekvo|!mK#CNl6+b6Ww@esX#|Nh?skCk3*
z43EeV4ic+kf!1O{{7<jxpQ*JAqhK`-4b8}Mb<|q7<btYdhcpO>?Xinen)`Kb_IoMc
z4~(Kt=iIe3)%G{hC$;@-);OFeD~wJUoZ6Q!=i26WE+RC6bji%jOhisZcUvQMe66e{
z9gMbriiy0F^Lq3_kmYQvMsb`sp;(y*2*L6g?_Vw=yNbi{=G3q5fbm4-Sq7kzt-hx<
z?xbzrVf3v+_tA>JlJAqpM11~$*B7?eDdB}}YcgJ#gwD*L{1U!7QohT<W5+Ag0@prR
zgALp=g)yjBeoY=aW=nI1xgU6flyB*ijEukGXj|fZ``b%OSgNIm0mWLMK79&|EA{0g
zwsXcOYt)r<+U)XGtHfxUwk_a(`Om9Y2a0}m`4JBmMq>;{8xe3GKPtD5H<rwylTvFl
z%&Q2p!8C&GK1%c`t74H>=^)*lqX~`^+-gf8P;9%btbU||Tc~rt|C2=J2ieVop4(1x
zE!$woRG#VNQ#*Q8%hmOTQ}!sNxx}75hVX%*yq5Cu8KU0VAm#T3_XZZ6=EzMUctDgR
zQ7(kQGip<eGRNSbSD}o^M8H5px69JrGnZTKOL+(%_Dfu(6-1nPv>vg|>yz~*Fl$f2
zjxYdi3X`?k!bsVv$Mu}|&FUIRzekV2_m5Wi3+Wmdgo2^QLwFJKr>JNOCdUR%z~D0c
zEkr)<e77(5u-}?Kv-zn<FzCel>b08)UeU^PzhVJ?i4aKxQv>O4QHj+qqLg%bvDfq%
zX1oMSHYq6y08tC3gc+a*t^K7sTw=?XBa!khI@3P~*{~E3tEv`wSkEltjJuT~c?{?x
zzHJ*B_^Vg0EZqu!mqq3+Zr=mx!Da@&pQD>wW5>VOvuVJN*J$+~5Lm|}i$^$O@`FP|
z``<j=gx+kxi-^|x@4w@{C43@E#1|z*LwSBClT}e}*Rd)^;O&LX%DOsTl<RL7rh-Fj
zaa>&7rLXUw)F?;$`{8OY!HcNztJYN3@T8d109xXX+uPg6?^MEv(#~Jsi<)1&D6xGz
z1Dz!uu|mIm`7-LIouix}c}CerlUN)@GD;pzt#kf-9DIy%dVYegDLBKS;$rl|`Y~}f
zfYFbIvzNeO|D_IH$c~;H+e)MkIlIcrhj2inqoe&rl+#Yy73nKp`Kbqx!Tqyo2!H<?
z&Rnrj^S!aLfcqN1_X;lkalAtE+i=8JuU?goTETxE&eY77stkAnJ1}^+qx2T@3F}?B
zkf5Wh8xk3rieVrDp9LHnfTR^5teE9F=}1WbA_-~fhHsxTUGO!D$qf_B;K-3%yn=Eu
zpc|hCX;-`^3E2*DPv!B!SBCN`*ym!GBBOYDn7kO*9vmH=j<YAF*_>l%c=~h<SaDDP
z2vxrVtWe%#jLD1J0S^#YP&kRqqABxY01z1{9fvLlW5M**<G|T4|E|i2aLei1=tMng
zo*YWJXherXP@<9(zkRdDB+x{S$nF!BUPOL45GU?*YC*o+@F`GLaWKSadM@v^zYq&%
zPU%Wc=6%TIY}Y~K!InH-%pr0?$B*B`IRD)L!Tuq>feILUNFS;Rut6*^;>bYMlDCpD
zNO^!;p-Vg8!#ZLYx;q&LaL3%5?n{g&wv#vv1C()~Mr)S>?g8Je0n-5_dcb66l`Iz?
zPc9uau8*jBM@Qgi##oh1$u?P8_M|7`8t1OEKVa0H;9AUE&o(qm{pTMTXMe`*K9PCh
zzQ@?tCpebbck?p=Vf1dBWWc@^xD8SN@a9~qWgEX@8tA~VL6tAuph8i-;~!K-EFiUr
zuLwA+PchnwmSbr~D_J(CMP7Ipm8z*cs|SS}h%#eUKx&^iHzzN^3Dv>5>*?;+?K6IW
z(MimIdTum&yh3?SyDouNck0Z_)k?X+i!X#mV<gs;76W_6dX?vmhK7a+H-sNS3`Ybm
zy?uPW1(#M)$!FpEknpVR&Pqc6Vcg)Ydq*K<h~~mmjv+Gam!T)Qt*+H?h>%>2cPGJL
z=EZjgkEBh<UUAD@`rh$xQ#Nh&sHIiY>sGdHt2{=J2AFf;Ax|pzyz=Z=5E{f~FuPG+
z#0me*IZZ8y@}8Ya0<klOl4Tt!Km7Y_fF|)IPKb1C`+%#?hM7Vn2gLPlQsw|fm;U;`
z2)8#dddSDP47C_R?FKu4S-d`xo^OGvz)m+agAjdvv8;R(j$6Mr&*|dq@cdkQ@!XUZ
zaU%hlVN(;s>5$z(wm0l<4(4IxLhL68ecM}Vnww8Pe*73?pyb=mCsVMmeZ(E)1ObE?
zj&%qT$Kax|z`-ys;P9t7UzDVBPhhHvZ{JRTfKVeG_y6dBUf=F}HUYp1Sj;dYUM~#l
zI_tnM2|$wf8tLs@{oZDGID8wtNY-o0*=ugHmmn>wFiBO|y{x>BlUez~`qR^J+k6I*
zJusTcM1+We1(giJLq-@>4k*po=sV?sPTE(b0XBh&;%;fm?O#r6CLz<b{Y>Vk`uJn~
zE`7C!OQu_gcUjhNLoei=mZ!2AK)~Ypxhs)!4hM{kjGBtB<bYBJL+ticJsSA>AvQ2B
zK|Q=5|A9}kA1TY#66D~ZA{yaUkWK|$^nD4pbB(g{Lgs3I{ya^1Iw88?u%OmP+zV8)
zq55>soOwck0c;X2v<=M&AZ>yoSi~osqC7VekS4s6-v_7|QI};=N5Wr5Ab2eCQb4&a
z;QLHrm?2YL#m}FZ80Db`D&;bgMOj`aei2D3Vb6H94c?{w{6eCmWk3z#5D0|vvagai
z)dy3yEgU(5EKHX;Er_kiFk~S}nnGdky7Ybn5-2Itd_9B*T9%gLZl#>|3|Mzad7+I+
z2@=tW8B%joWBw$*3{1N}HSZiNH+@D?7qbzAej0QomLE+l)z|md{f-7IMCf{75pO0T
z4oKD-i7tmHX#utkk41$i2h5$3*f8;{b=W%*J`W5`HC9aY3-{3r5VmDwx0#C;ETKA5
zLZwb70C8^KxIuE7K(O$q(5Z++kAgwzq9|c@Z2;!0S60%M(ot(w)N@dQ0F(id`CmkB
z@4Ec`^k7@jNn%YLdIKT>|92z5`rTvueJHBtH$SUS3IfmE3RTp5jxZ*$T&U2|;u%$c
zLEwtI>0W6c<I|v^WF3L<Pg1L-T!=pg#3MYYq}Y+(fWZsF;06>r@%!P);J~$SJ{Ecj
zqzWza0SrBYwX%886=w#d!SrN0_12^btyQ_71Fy~nMmaEm`1L3Z`Vv0V<5#Dz3m_OS
zJ3|MRHOjTTeqRa~K>Hx>f{bB{x&#XxyCn^?ZR5s`&IpFl9qMNEY_@UMr6;IiaST}q
zjQ-kqzyDE1OixU}2S7tK;IlKOWGdpuOx`Ka7n0Nu?&ARPkPJV}sDb{64$H|9T?aQd
z<l)07@R4mkWgHIQW7hzakwhWTmxeTAjx-`bac9Y}S2LOi!B1SrKT?*@VBe5okMR>2
zp{hWPDD=YD?8Z~N5M(5(Cj!EkPvpz>0D-%tY|ns=iAzfxVE)&~s`B-AW46m*{vHD#
zC(C8Pz@{}{pG1N0q7+T#q<orA<C<u~7cwtYiz^NRM)}P=wGmuS!3#<ODP%#BVIq;V
z#B`UG^g<HO*Vxyo@7ey>Us}*~eRj51Mxkh1G6T4X5kc}>d?StYlDYhS`}RcyoD79a
zK>SS}K{>5EcZeNZ2S$-~u7a0e<DOqRefTTovoj7x4W_AnQ^SSK<sy{yCqlb0FA{V-
zk0yG*hWM<6)QtB0ac9Mf6>tqaP`uTLgSvE<18a8@XBzwr-A2C5;lnrOPV;s+f>F<{
zOEhBymUe0td5nuO&YU?jg`xKZ^1dJ7ut#xY@*j!xX%0%f^6gt4M3|_&)oa$cz8(!V
zhmZ3nLlX6d3hFI?O4&arlYq2@rs8<;PC1Qe@9ZlRw)g=Yab-n?DzFe%=tEC0FVlUp
z|D!LTNzykzH@kY(syHY#GTqQ92H@*#IL8_}oLYnaLbO5x!D>ueGX6(Le*AD$94T8#
zoHoSl&ahLhetfPZJ;a?@9^WN>g7Q(rs)MwdnHgzd7}8-1ZfQgbf<#r7mL7!DrUl0W
zX;un)p5!ns*Z7feWmyD~vHb615jhbFj0~&+jkjGHrU~~W?r~1LH;cslaKfN+r6blp
zB5|l*26!Y%x?%M;z|SK#09XP{7lrxCH*el!4u!x7cIbX>Q5(ld6JD>v^TN4L2MJ#?
z5Vxjcw`~iPlsW#uCiy>BWGvu@czk($F$<8wBS_1=zqztohs;bh65cz1)gn7hVq%Bj
z&g6dyvx;~JlW*YtlY>MM014{hF7s<}$L#)ybe&j&q)stg0aXZq(h2cWocq;X<5y<~
z&{2<8GG#Cr$XJe`b|ZjLIMr1mB4&L%7!U=sGxQ8`$YFMm54AQQ4&r301HWB<eLsUt
zugW%O0LC2IdV;4Nt(LQ3Vq&5Z)3P3H5jp6wa+}K>lm?anf(r_+#GgFVH0avfnx>)P
zwKb*Afs=aF<*zXG9$xtv!b<F1c@cpM!j1^zF97PRJX2d+It>lRhD-YU{h&O2>HxQp
zl2z?9&JRxxb22qRSQ|*rwS3bN;_VRPj|O*!QDe|#`{C(f;!&ZkO?GWhCcM<1WQ}C`
z2=0er0A~#EA4`_J$KqJ(x}6sgnW;CLINcY%+`+*})Zof)rU%v3)!8s;k~z2*<z6;6
z5y#d~N7k9rMXlb`gFqv&7=TbstdOh1#X-(oPneNo5Y^eW=MrU(J|cl87Sk+#eeX<L
z%g4Trl99tyCqq3@Bl1GekoOm}mE`s8JDT%nk?LE)7)GNKvkOMV=K*>^0waWxJht-!
zG#gR9zw&gspRJImR04?!he@U4T3`ykrU&MMqe;|=l!tSyfuPzD?`5jE#m5+Ey_=8*
zz~U-28&VU@GHgQ~!rzmy1{7!mA74PG_qJ<czIeo=fd;9@@dO7R(yGEAK;p=T8Wc#+
zmoLd#S@rI2Da3-qKp+*k3z6ynq?Dv&Fvgr1fhuT2L%4|;CDJl7GQ{RZvW$%yZTPDM
z+i~EL%aS@V0DQ#y0NHhhW?nt}gu0-4<&lxc>%X6X3GX#c|J=TJuMvqQz+J$?d<0T*
zb5D?=FgN**PKF@#`20|nWI-T=#5z@Q{YDHgP=Z)IFw%a+zc1P6x>o=SDMsxuJ%0s%
z*@rEry}eSdnphNpb0UFzA4i&zbwHjnxfd9?0-Hhy?25)ichd-^33+};&t*C(K81Bc
z_RdJq0LLM%=%+Ei=DlNw(K})cBn!25n#!r@V}shH@<p2S;V#OQi7h-8kNMdNS@$33
z0d-%YREW8q252bo*rm5~C}|1xK&;j2og(pcM&Feb&klVu*uU7P|Flu)Tl(!=Jp`n)
zeVVimNEh+wy#RUwJI>w+l{9(^I5Z5nV_8NTLl?;=W=hRuVQFW_CFJ~HfMM{yNX)g%
z{Oziuxv2r7hA@}5%Gxy5v%8;!*>iJHJ|w&fP?|al$I&@W$+HQ*g#L*WFc6g&E}$12
zotjTH=Xe>frlHXQMni;hd_n?tYFyNiw_+^7$BMUY&36zmkit$Tm1W;~p!%O32O+q!
z;^s@HzN1&iumCjh$jtj;t8F6?V;tiJLH^tFm0Zj|qB&6wrs`q49-<b(uIgj<bG5oQ
z+mM`_o|&B`;`4Vxtf#jV-4{cGl>Ua)aAHC-5)DM!kfsW|e_l>xZceLfI1yuDxn%nZ
zu<SZgMp@>ENEi~=1Px2px&W8#5}NMK-@q&)CD@Q=$|6(+c0x@BNr?tOKL4GE57GGp
z(u`WNwnw)3;C^8MII?po%d7C>Kp@{>wyoc7{T?a?i+4%H`9XkzT!2Kf#Mi!kdx@+?
z?TtERt3RX<ET)(sR$)8N-LlQ=S;I(78@cp9ewB53<m}|HUq=JK)I_onuMuyO*!?~9
z8+;V}JAlTv9XqtK+Dh;U#t&K51FK6yH^AI^`|e#Mtl^k}^pZCArdR>{1q-PDTvFLs
ziUF<2?+HEm*<Nleot`}g0{Q`3J#s{p0dtO|+K6~f_x$B!EQ&2f;=MV)Q=Ir3b!7;)
zFpcCNY*Yf@91;#U-uN~|0<Kl}D+}#NIZjT>B?1n?xd6q4y-r`g68DTuQLJPz0tw4T
zxY6rMEb}(cQFW9}zT<gH*2LSp81KGU>Qvzgl3Bpg738g0H^2*cK?JIfMsjR@#&ng$
zMilK5fto>F$|7?NW_v!W*ag~1arD2P&_xc?R48UHJj_0TqtoxZhSIXaR*c6`UT-~O
zQ792)fP1^U$9lWly(S{bzVOqh8x8&aEU~CVX;pwTqkds*c-PKMj^ZH_q)qb)V+wgO
zO@9v3go4@vTcDhbKPl2Z;4~IkJc9Y*K5TvV+&MB2XOX81$l1qKE_YZ5sL`qbGr}j|
z^!2qG`}!dyDpT?pu?vunLj%Fitwzf}vIfKKZ!0PQ)YnYS^lidS7bRo^j@og5->Q2Q
z2TdUpQ+QVr8i!WL%vGTO_>Pvo@yAM&y&HXu<5E(N1Go&0j2@HL&m&8{n4rXJy<XSp
z*e$S(Vs^vx7OPw*68L$@yMGih{Vrm3kt|jC%C{Z82X2DB;YaFU*qp>Z^viW3(AX)>
z7tWvm3YGeMYw+Jc9P%7UMak{iv!{3VCpR}Y5j=<lh@l0<Wja@|K+`I3;)Hmx{pE}#
z@Y2*+AF;Zj4GsFZ9+-*Tk3`*r#g#9AV0??>@R7I1hB0AvFo&{0im=%xH>^|#x0LX9
zND1;S>3YODPEqMR6#vtU@1pWYUaucVRw$0-n`_aZT39wiQe>Xf?1CN3A`qM@Tm$+|
z%rF4AE025dB`XNfP-v3;pG{Wm)RTM0zkWCfJB%nx;_M<I??;xU*lTJ@b)1D4s0uJ-
zl6(|Igw{fmNpJ})KHI5RgFG7Krx=G0_wN&VMiN&@0(fh+Xeyw}yxKb~i=4!kji?f(
zSHNU$^ZNepV!}ebe|0V<z(N=bJh1u#H#;67c|8bi7OSGn{3R-aSxSOZr~~8SlKJAm
zpdgx^V`ru~{VlV<E-?URRS!UBA@`cEx7OR*k25LDJd*Z=gb>FQ&5ZEHVf<+Bae`@J
z>*OH3AoZ4l;Q;(dbOsmX5Rn=F{;uOlkk-O)%;8iWJb-0-p$kT#l2wjv4i_+^$v)?o
zchTuI(0w}=1N``946sR#I+0t=P1#rC6I9|g1^9&I<LrtjkC4SoPy+Zw%0;*$sF?1R
z4C>`(;uc-ZHm_Qv9<|R1u?6uh$R;w9|6C7AN`<FF4f;k+gO+Av9pCHcr;Emgm!Fg6
zCm?CT<-aLf#pWepks%RE<^5yhqPVyGZhwVUt5^4QRZ>(ohbs95^UO~`H?phkHPx4x
zCcK!^ItI$DrHT&h)%JHHzY&0I#uCnH`@bm_8<1TUAyf@sC*I0_n<0sx!xsC|Kr|nT
z|6GRKTJ_IiU4I>&#(_$`1SH2HY9T6uGeojY6i)PVv5Ei@L%79c$sL@VVyfrz$WZk=
z;Ntuk9}f{jHpeor3cY*po*qm|tf)W10vTCqN7B%QREJn7gbm_nq@x5Wt@Y%Y7bBv$
z1Y96A)q9fe;U^0L05I^yFfs}8=1JIvAUp@|mEYDyRQ(zl1NA$QeFmBv;pkaom+^mG
zzkPfHT6hU0kgOXK*9*iC8)<o$jBy0y*C&L?uQp&?`YXymr20aTxe>@V%XzROo!D3K
zL$H6jnTmNqWR+HhHC8vF5y9*p0m2YHOGqi5z<2;K^+rCs%}1%Y{O<>0ttKJs$HLvc
zSYd>7kRH-ue^Offzo*gz&@xi7qN1HCi6*lRuWnG3iMT&`oWO_(NhItrUtsP1j!%4P
zULmqagcZ4&h{mfFrDclBV;(d}Ie}p%g91Nf)sXT5>-n*ObcPPjjZ`_-%E+uB92hwQ
zyd{%;vQ*6OPVzl`#|9cTF~w`)N2h=^4z`z=!D}Q7iG(!5Jt1o~BsIdu=xZ7qj-aM&
z_Bk%fa{2G4%m56uUKuX3L^O`fDFMlKuM=?+T)+K8e10-60XXEhCAbq8QuRmR@cGD>
zVX%>52Wb+)AJJ^RJ8c%f#8OIneiK@ijVuc>zI=D7t@7xF+n}f<Hh2sEw`%?Rfolv7
zDtjGDJ%(4#E3UNDUhl3y9-UDOYd#x%iwJLg2cn8GN`fCg>`Z=$SFHN|if5m!;S#~O
zNNPwoONX2#hClHNciXk4fFO~*3p|{IgTw#D<|*dQ-_xBgBiskSu|mujd55kWcz>cH
z;Fyv#@Ic~E4t##LeO@1HNT!O?R7_~s64c{zpa?la?4gV!7-b=HBQkw-5<Jj|Y;0rf
zb0;pK>Rm6M@>)VMUEmr4@K>!{mxMsivh&qVqG-wfkr)H+DiSt?7k|2yaZ)A*q~1AD
zS$s;rVn}{Yev6{rbM!v3pttKLFhPoN92mCx=vZvlQT2`{aJO8J<9P^e#`cxw{M<wd
zlvg~~PA}nNNRMDiV7Sm~4Mnl<Z$c+lTy}S-$)KOYnk52LI}(kOxk(&Add{;8ySYZc
zr$7u)7%j?>&;H&)q$C6wIWLc{J6p%_-K#1pD#~>MC!0WD`=8aRDiEB3-CKlhNirF*
z*c*#TqY1t_Whsgu>d-rJf{9b_=dV#14t|6rvx0-FKY;Izf6nDD5eoWO)NRc?U+dI&
VbL)x}7Wr$&VU0s+2Tou6{{Wt=HT3`h

diff --git a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
index 7f00493b6c..73d8e50cc8 100644
--- a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
+++ b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
@@ -1,6 +1,6 @@
 		test	test
 		ndata	$\chi^2/ndata$
 group	dataset		
-NMC	NMC	204	1.6064258130123148
-ATLAS	ATLASTTBARTOT	3	1.9377358549444097
-CMS	CMSZDIFF12	28	1.8519163102572869
+NMC	NMC	204	1.6064258192222478
+ATLAS	ATLASTTBARTOT	3	1.9383178404814991
+CMS	CMSZDIFF12	28	1.8520435691568384
diff --git a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
index a669de4f65..0db9efa215 100644
--- a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
+++ b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
@@ -1,6 +1,6 @@
 name	Total	NMC	ATLAS	CMS
 npoints	236	204	4	28
-0	1.6331414605961554	1.6064258130123148	1.4642155397440857	1.851916310257288
+0	1.6331643952670938	1.6064258192222478	1.4646775563260181	1.8520435691568402
 1	1.7334090365532264	1.734962043323776	1.109948798376749	1.8111600212501475
 2	1.6877884536672854	1.6453401923939792	2.065149379329199	1.9431456535639562
 3	1.6987177907128304	1.6814300483245213	1.2167235538564396	1.8935276619499943

From c5425ffbbd777dcdcbcc50f0afc769d6d522fb0f Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Tue, 15 Feb 2022 18:08:48 +0100
Subject: [PATCH 04/13] make N3PDF closer to vp PDF (and add a N3LHAPDFSet)

---
 n3fit/src/n3fit/io/writer.py              |  10 +-
 n3fit/src/n3fit/tests/test_vpinterface.py |   9 +-
 n3fit/src/n3fit/vpinterface.py            | 301 ++++++++++++----------
 3 files changed, 177 insertions(+), 143 deletions(-)

diff --git a/n3fit/src/n3fit/io/writer.py b/n3fit/src/n3fit/io/writer.py
index a68cc33404..7d0a5fb7e7 100644
--- a/n3fit/src/n3fit/io/writer.py
+++ b/n3fit/src/n3fit/io/writer.py
@@ -11,6 +11,8 @@
 from reportengine.compat import yaml
 import validphys
 import n3fit
+from n3fit import vpinterface
+from n3fit.vpinterface import integrability_numbers, compute_arclength
 
 
 class WriterWrapper:
@@ -57,9 +59,9 @@ def write_data(self, replica_path_set, fitname, tr_chi2, vl_chi2, true_chi2):
                 chi2 of the replica to the central experimental data
         """
         # Compute the arclengths
-        arc_lengths = self.pdf_object.compute_arclength()
+        arc_lengths = vpinterface.compute_arclength(self.pdf_object)
         # Compute the integrability numbers
-        integrability_numbers = self.pdf_object.integrability_numbers()
+        integrability_numbers = vpinterface.integrability_numbers(self.pdf_object)
         # Construct the chi2exp file
         allchi2_lines = self.stopping_object.chi2exps_str()
         # Construct the preproc file (the information is only in the json file)
@@ -140,8 +142,8 @@ def jsonfit(replica_status, pdf_object, tr_chi2, vl_chi2, true_chi2, stop_epoch,
     all_info["erf_vl"] = vl_chi2
     all_info["chi2"] = true_chi2
     all_info["pos_state"] = replica_status.positivity_status
-    all_info["arc_lengths"] = pdf_object.compute_arclength().tolist()
-    all_info["integrability"] = pdf_object.integrability_numbers().tolist()
+    all_info["arc_lengths"] = vpinterface.compute_arclength(pdf_object).tolist()
+    all_info["integrability"] = vpinterface.integrability_numbers(pdf_object).tolist()
     all_info["timing"] = timing
     # Versioning info
     all_info["version"] = version()
diff --git a/n3fit/src/n3fit/tests/test_vpinterface.py b/n3fit/src/n3fit/tests/test_vpinterface.py
index 3664ee2134..6b9edf0d15 100644
--- a/n3fit/src/n3fit/tests/test_vpinterface.py
+++ b/n3fit/src/n3fit/tests/test_vpinterface.py
@@ -6,7 +6,7 @@
 from hypothesis import given, settings, example
 from hypothesis.strategies import integers
 from validphys.pdfgrids import xplotting_grid, distance_grids
-from n3fit.vpinterface import N3PDF
+from n3fit.vpinterface import N3PDF, integrability_numbers, compute_arclength
 from n3fit.model_gen import pdfNN_layer_generator
 
 
@@ -39,19 +39,18 @@ def test_N3PDF(members, layers):
     xsize = np.random.randint(2, 20)
     xx = np.random.rand(xsize)
     n3pdf = generate_n3pdf(layers=layers, members=members)
-    assert len(n3pdf) == members + 1
-    assert n3pdf.load() is n3pdf
+    assert len(n3pdf) == members
     w = n3pdf.get_nn_weights()
     assert len(w) == members
     assert len(w[0]) == 16 + (layers + 1) * 2  # 16=8*2 preprocessing
     ret = n3pdf(xx)
     assert ret.shape == (members, xsize, 14)
-    int_numbers = n3pdf.integrability_numbers()
+    int_numbers = integrability_numbers(n3pdf)
     if members == 1:
         assert int_numbers.shape == (5,)
     else:
         assert int_numbers.shape == (members, 5)
-    assert n3pdf.compute_arclength().shape == (5,)
+    assert compute_arclength(n3pdf).shape == (5,)
     # Try to get a plotting grid
     res = xplotting_grid(n3pdf, 1.6, xx)
     assert res.grid_values.data.shape == (members, 8, xsize)
diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py
index 9f53b786f2..7234be6391 100644
--- a/n3fit/src/n3fit/vpinterface.py
+++ b/n3fit/src/n3fit/vpinterface.py
@@ -24,6 +24,7 @@
 import numpy.linalg as la
 from validphys.core import PDF, MCStats
 from validphys.pdfbases import ALL_FLAVOURS, check_basis
+from validphys.lhapdfset import LHAPDFSet
 from validphys.arclength import integrability_number, arc_lengths
 
 log = logging.getLogger(__name__)
@@ -45,9 +46,111 @@
     "T35",
 ]
 
+
 class N3Stats(MCStats):
+    """The PDFs from n3fit are MC PDFs
+    however, since there is no grid, the CV has to be computed manually"""
+
+    def error_members(self):
+        return self.data
+
     def central_value(self):
-        return self.data[0]
+        return np.mean(self.data, axis=0)
+
+
+class N3LHAPDFSet(LHAPDFSet):
+    """Extension of LHAPDFSet using n3fit models"""
+
+    def __init__(self, name, pdf_models):
+        log.info("Creating LHAPDF-like n3fit PDF")
+        self._error_type = "replicas"
+        self._name = name
+        self._lhapdf_set = pdf_models
+        self._flavors = None
+        self.basis = check_basis("evolution", EVOL_LIST)["basis"]
+
+    def xfxQ(self, x, Q, n, fl):
+        """Return the value of the PDF member for the given value in x"""
+        return self.grid_values([fl], [x]).squeeze()[n]
+
+    def __call__(self, xarr, flavours=None, replica=None):
+        """Uses the internal model to produce pdf values for the grid
+        The output is on the evolution basis.
+
+        Parameters
+        ----------
+            xarr: numpy.ndarray
+                x-points with shape (xgrid_size,) (size-1 dimensions are removed)
+            flavours: list
+                list of flavours to output
+            replica: int
+                replica whose value must be returned (by default return all members)
+                replica 0 corresponds to the central value
+
+        Returns
+        -------
+            numpy.ndarray
+                (xgrid_size, flavours) pdf result
+        """
+        if flavours is None:
+            flavours = EVOL_LIST
+        # Ensures that the input has the shape the model expect, no matter the input
+        # as the scaling is done by the model itself
+        mod_xgrid = xarr.reshape(1, -1, 1)
+
+        if replica is None or replica == 0:
+            # We need generate output values for all replicas
+            result = np.concatenate([m.predict([mod_xgrid]) for m in self._lhapdf_set], axis=0)
+            if replica == 0:
+                # We want _only_ the central value
+                result = np.mean(result, axis=0, keepdims=True)
+        else:
+            result = self._lhapdf_set[replica - 1].predict([mod_xgrid])
+
+        if flavours != "n3fit":
+            # Ensure that the result has its flavour in the basis-defined order
+            ii = self.basis._to_indexes(flavours)
+            result[:, :, ii] = result
+        return result
+
+    def grid_values(self, flavours, xarr, qmat=None):
+        """
+        Parameters
+        ----------
+            flavours: numpy.ndarray
+                flavours to compute
+            xarr: numpy.ndarray
+                x-points to compute, dim: (xgrid_size,)
+            qmat: numpy.ndarray
+                q-points to compute (not used by n3fit, used only for shaping purposes)
+
+        Returns
+        ------
+            numpy.ndarray
+            array of shape (replicas, flavours, xgrid_size, qmat) with the values of
+                the ``pdf_model``(s) evaluated in ``xarr``
+        """
+        n3fit_result = self(xarr.reshape(1, -1, 1))
+
+        # The results of n3fit are always in the 14-evolution basis used in fktables
+        # the calls to grid_values always assume the result will be LHAPDF flavours
+        # we need then to rotate them to the LHAPDF-flavour basis,
+        # we don't care that much for precision here
+        to_flav = la.inv(self.basis.from_flavour_mat)
+        to_flav[np.abs(to_flav) < 1e-12] = 0.0
+        flav_result = np.tensordot(n3fit_result, to_flav, axes=(-1, 1))
+        # Now drop the indices that are not requested
+        requested_flavours = [ALL_FLAVOURS.index(i) for i in flavours]
+        flav_result = flav_result[..., requested_flavours]
+
+        # Swap the flavour and xgrid axis for vp-compatibility
+        ret = flav_result.swapaxes(-2, -1)
+        # If given, insert as many copies of the grid as q values
+        ret = np.expand_dims(ret, -1)
+        if qmat is not None:
+            lq = len(qmat)
+            ret = ret.repeat(lq, -1)
+        return ret
 
 
 class N3PDF(PDF):
@@ -67,7 +170,6 @@ class N3PDF(PDF):
 
     def __init__(self, pdf_models, fit_basis=None, name="n3fit"):
         self.fit_basis = fit_basis
-        self.basis = check_basis("evolution", EVOL_LIST)["basis"]
 
         if isinstance(pdf_models, Iterable):
             self._models = pdf_models
@@ -75,28 +177,14 @@ def __init__(self, pdf_models, fit_basis=None, name="n3fit"):
             self._models = [pdf_models]
 
         super().__init__(name)
-        # The number of members will be minimum 2, needed for legacy compatibility
-        # Note that the member 0 does not exist as an actual model as it corresponds
-        # to the average of all others
-        self._info = {"ErrorType": "replicas", "NumMembers": len(self._models) + 1}
-
-    def get_members(self):
-        return max(len(self._models), 2)
-
-    @property
-    def stats_class(self):
-        """The stats class of N3PDF is always MCStats"""
-        return N3Stats
+        self._stats_class = N3Stats
+        self._lhapdf_set = N3LHAPDFSet(self.name, self._models)
+        # Since there is no info file, create a fake `_info` dictionary
+        self._info = {"ErrorType": "replicas", "NumMembers": len(self._models)}
 
     def load(self):
-        """Many vp functions ask for a LHAPDF pdf
-        from nnpdflib, this class fakes it until a time in which vp is free from C++
-        """
-        return self
-
-    def get_member_model(self, replica):
-        """Return the member corresponding to the given replica"""
-        return self._models[replica - 1]
+        """If the function needs an LHAPDF object, return a N3LHAPDFSet"""
+        return self._lhapdf_set
 
     def get_nn_weights(self):
         """Outputs all weights of the NN as numpy.ndarrays"""
@@ -111,7 +199,8 @@ def get_preprocessing_factors(self, replica=None):
         # If no replica is explicitly requested, get the preprocessing layer for the first model
         if replica is None:
             replica = 1
-        preprocessing_layers = self.get_member_model(replica).get_layer_re(r"pdf_prepro_\d")
+        # Replicas start counting in 1 so:
+        preprocessing_layers = self._models[replica - 1].get_layer_re(r"pdf_prepro_\d")
         if len(preprocessing_layers) != 1:
             # We really don't want to fail at this point, but print a warning at least...
             log.warning("More than one preprocessing layer found within the model!")
@@ -140,7 +229,7 @@ def get_preprocessing_factors(self, replica=None):
         return alphas_and_betas
 
     def __call__(self, xarr, flavours=None, replica=None):
-        """Uses the internal model to produce pdf values.
+        """Uses the internal model to produce pdf values for the grid
         The output is on the evolution basis.
 
         Parameters
@@ -158,122 +247,66 @@ def __call__(self, xarr, flavours=None, replica=None):
             numpy.ndarray
                 (xgrid_size, flavours) pdf result
         """
-        if flavours is None:
-            flavours = EVOL_LIST
-        # Ensures that the input has the shape the model expect, no matter the input
-        # as the scaling is done by the model itself
-        mod_xgrid = xarr.reshape(1, -1, 1)
-
-        if replica is None or replica == 0:
-            # We need generate output values for all replicas
-            result = np.concatenate([m.predict([mod_xgrid]) for m in self._models], axis=0)
-            if replica == 0:
-                # We want _only_ the central value
-                result = np.mean(result, axis=0, keepdims=True)
-        else:
-            result = self.get_member_model(replica).predict([mod_xgrid])
+        return self._lhapdf_set(xarr, flavours=flavours, replica=replica)
 
-        if flavours != "n3fit":
-            # Ensure that the result has its flavour in the basis-defined order
-            ii = self.basis._to_indexes(flavours)
-            result[:, :, ii] = result
-        return result
-
-    def grid_values(self, flavours, xarr, qmat=None):
-        """
-        Parameters
-        ----------
-            flavours: numpy.ndarray
-                flavours to compute
-            xarr: numpy.ndarray
-                x-points to compute, dim: (xgrid_size,)
-            qmat: numpy.ndarray
-                q-points to compute (not used by n3fit, used only for shaping purposes)
-
-        Returns
-        ------
-            numpy.ndarray
-            array of shape (replicas, flavours, xgrid_size, qmat) with the values of
-                the ``pdf_model``(s) evaluated in ``xarr``
-        """
-        n3fit_result = self(xarr.reshape(1, -1, 1))
-
-        # The results of n3fit are always in the 14-evolution basis used in fktables
-        # the calls to grid_values always assume the result will be LHAPDF flavours
-        # we need then to rotate them to the LHAPDF-flavour basis,
-        # we don't care that much for precision here
-        to_flav = la.inv(self.basis.from_flavour_mat)
-        to_flav[np.abs(to_flav) < 1e-12] = 0.0
-        flav_result = np.tensordot(n3fit_result, to_flav, axes=(-1, 1))
-        # Now drop the indices that are not requested
-        requested_flavours = [ALL_FLAVOURS.index(i) for i in flavours]
-        flav_result = flav_result[..., requested_flavours]
 
-        # Swap the flavour and xgrid axis for vp-compatibility
-        ret = flav_result.swapaxes(-2, -1)
-        # If given, insert as many copies of the grid as q values
-        ret = np.expand_dims(ret, -1)
-        if qmat is not None:
-            lq = len(qmat)
-            ret = ret.repeat(lq, -1)
-        return ret
+# Utilities and wrapper to avoid having to pass around unnecessary information
+def integrability_numbers(n3pdf, q0=1.65, flavours=None):
+    """Compute the integrability numbers for the current PDF
+    using the corresponding validphys action
 
-    # Utilities
-    def integrability_numbers(self, q0=1.65, flavours=None):
-        """Compute the integrability numbers for the current PDF
-        using the corresponding validphys action
+    Parameters
+    ----------
+        q0: float
+            energy at which the integrability is computed
+        flavours: list
+            flavours for which the integrability is computed
 
-        Parameters
-        ----------
-            q0: float
-                energy at which the integrability is computed
-            flavours: list
-                flavours for which the integrability is computed
+    Returns
+    -------
+        np.array(float)
+            Value for the integrability for each of the flavours
 
-        Returns
-        -------
-            np.array(float)
-                Value for the integrability for each of the flavours
+    Example
+    -------
+    >>> from n3fit.vpinterface import N3PDF
+    >>> from n3fit.model_gen import pdfNN_layer_generator
+    >>> fake_fl = [{'fl' : i, 'largex' : [0,1], 'smallx': [1,2]} for i in ['u', 'ubar', 'd', 'dbar', 'c', 'cbar', 's', 'sbar']]
+    >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl)
+    >>> n3pdf = N3PDF(pdf_model)
+    >>> res = n3pdf.integrability_numbers()
+    """
+    if flavours is None:
+        flavours = ["V", "T3", "V3", "T8", "V8"]
+    return integrability_number(n3pdf, [q0], flavours=flavours)
 
-        Example
-        -------
-        >>> from n3fit.vpinterface import N3PDF
-        >>> from n3fit.model_gen import pdfNN_layer_generator
-        >>> fake_fl = [{'fl' : i, 'largex' : [0,1], 'smallx': [1,2]} for i in ['u', 'ubar', 'd', 'dbar', 'c', 'cbar', 's', 'sbar']]
-        >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl)
-        >>> n3pdf = N3PDF(pdf_model)
-        >>> res = n3pdf.integrability_numbers()
-        """
-        if flavours is None:
-            flavours = ["V", "T3", "V3", "T8", "V8"]
-        return integrability_number(self, [q0], flavours=flavours)
 
-    def compute_arclength(self, q0=1.65, basis="evolution", flavours=None):
-        """
-        Given the layer with the fit basis computes the arc length
-        using the corresponding validphys action
+def compute_arclength(self, q0=1.65, basis="evolution", flavours=None):
+    """
+    Given the layer with the fit basis computes the arc length
+    using the corresponding validphys action
 
-        Parameters
-        ----------
-            pdf_function: function
-                pdf function has received by the writer or ``pdf_model``
-            q0: float
-                energy at which the arc length is computed
-            basis: str
-                basis in which to compute the arc length
-            flavours: list
-                output flavours
+    Parameters
+    ----------
+        pdf_function: function
+            pdf function has received by the writer or ``pdf_model``
+        q0: float
+            energy at which the arc length is computed
+        basis: str
+            basis in which to compute the arc length
+        flavours: list
+            output flavours
 
-        Example
-        -------
-        >>> from n3fit.vpinterface import N3PDF
-        >>> from n3fit.model_gen import pdfNN_layer_generator
-        >>> fake_fl = [{'fl' : i, 'largex' : [0,1], 'smallx': [1,2]} for i in ['u', 'ubar', 'd', 'dbar', 'c', 'cbar', 's', 'sbar']]
-        >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl)
-        >>> n3pdf = N3PDF(pdf_model)
-        >>> res = n3pdf.compute_arclength()
-        """
-        if flavours is None:
-            flavours = ["sigma", "gluon", "V", "V3", "V8"]
-        ret = arc_lengths(self, [q0], basis, flavours)
-        return ret.stats.central_value()
+    Example
+    -------
+    >>> from n3fit.vpinterface import N3PDF
+    >>> from n3fit.model_gen import pdfNN_layer_generator
+    >>> fake_fl = [{'fl' : i, 'largex' : [0,1], 'smallx': [1,2]} for i in ['u', 'ubar', 'd', 'dbar', 'c', 'cbar', 's', 'sbar']]
+    >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl)
+    >>> n3pdf = N3PDF(pdf_model)
+    >>> res = n3pdf.compute_arclength()
+    """
+    if flavours is None:
+        flavours = ["sigma", "gluon", "V", "V3", "V8"]
+    ret = arc_lengths(self, [q0], basis, flavours)
+    return ret.stats.central_value()

From b0b8a0f272a95a554b7d9a2866f34249e9f4ea25 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Thu, 17 Feb 2022 09:59:38 +0100
Subject: [PATCH 05/13] phi needs actually the mean of the error members

---
 validphys2/src/validphys/results.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index 5060ac9455..ba84c518c4 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -547,8 +547,7 @@ def phi_data(abs_chi2_data):
     1410.8849
     """
     alldata, central, npoints = abs_chi2_data
-    cv = float(alldata.central_value())
-    return (np.sqrt((cv - central) / npoints), npoints)
+    return (np.sqrt((alldata.error_members().mean() - central) / npoints), npoints)
 
 
 def dataset_inputs_phi_data(dataset_inputs_abs_chi2_data):

From bca31e56879ac85f74781d12cc619edb03934233 Mon Sep 17 00:00:00 2001
From: "Juan M. Cruz-Martinez" <juacrumar@lairen.eu>
Date: Thu, 17 Feb 2022 15:58:49 +0100
Subject: [PATCH 06/13] Update n3fit/src/n3fit/vpinterface.py

Co-authored-by: Zaharid <zk261@cam.ac.uk>
---
 n3fit/src/n3fit/vpinterface.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py
index 7234be6391..6f58b1cad6 100644
--- a/n3fit/src/n3fit/vpinterface.py
+++ b/n3fit/src/n3fit/vpinterface.py
@@ -62,7 +62,7 @@ class N3LHAPDFSet(LHAPDFSet):
     """Extension of LHAPDFSet using n3fit models"""
 
     def __init__(self, name, pdf_models):
-        log.info("Creating LHAPDF-like n3fit PDF")
+        log.debug("Creating LHAPDF-like n3fit PDF")
         self._error_type = "replicas"
         self._name = name
         self._lhapdf_set = pdf_models

From 01bd364694adbb4ee1c5bd197e86ceb680a4178a Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Wed, 2 Mar 2022 11:21:31 +0100
Subject: [PATCH 07/13] use the raw data for the correlations

---
 validphys2/src/validphys/correlations.py | 34 +++++++++++++-----------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py
index edafd89fef..2841d2bb4b 100644
--- a/validphys2/src/validphys/correlations.py
+++ b/validphys2/src/validphys/correlations.py
@@ -10,15 +10,16 @@
 from reportengine import collect
 
 from validphys.core import Stats
+from validphys.checks import check_pdf_is_montecarlo
 
 #This would be a good candidate to be optimized to calculate everything in one
 #pass over x,
-def _basic_obs_pdf_correlation(pdf_stats, obs_stats):
+def _basic_obs_pdf_correlation(pdf_val, obs_val):
     """Calculate the correlation between pdfs and observables.
-    The expected format is are Stats instances of:
+    The expected format is two arrays
 
-    obs_stats: (nbins x nreplicas), as returned from thresult.
-    pdf_stats: (nreplicas x nf x nx), as returned from xplotting_grid.grid_values
+    obs_val: (nbin x nreplicas) as returned from thresults.error_members
+    pdf_val: (nreplicas x nf x nf) as returned from xplotting_grid.grid_values.error_members
 
     The returned array contains the PDF correlation between
     the value of the obsevable and the PDF at the corresponding point in (fl,x)
@@ -26,8 +27,8 @@ def _basic_obs_pdf_correlation(pdf_stats, obs_stats):
     (nbins x nf x nx), compatible with grid_values, upon
     changing replicas->bins.
     """
-    x = pdf_stats.error_members() - pdf_stats.central_value()
-    y = obs_stats.error_members() - obs_stats.central_value()
+    x = pdf_val - np.mean(pdf_val, axis=0)
+    y = (obs_val - np.mean(obs_val, axis=-1, keepdims=True)).T
 
     #We want to compute:
     #sum(x*y)/(norm(x)*norm(y))
@@ -44,34 +45,37 @@ def _basic_obs_pdf_correlation(pdf_stats, obs_stats):
 
 def _basic_obs_obs_correlation(obs1, obs2):
     """Calculate the correlation between two observables. The expected format is
-    Stats instances of:
-
-    obs1: (nreplicas, nbins1)
-    obs2: (nreplicas, nbins2)
+    arrays instances of:
+    
+    obs1: (nbins, nreplicas)
+    obs2: (nbins, nreplicas)
 
     The result is (nbins1 , nbins2), a mareix containing the correlation
     coefficients between the two sets.
     """
-    x = (obs1.error_members() - obs1.central_value()).T
-    y = obs2.error_members() - obs2.central_value()
+    x = obs1 - np.mean(obs1, axis=1, keepdims=True)
+    y = (obs2 - np.mean(obs2, axis=1, keepdims=True)).T
 
     return x@y/np.outer(la.norm(x,axis=1),la.norm(y,axis=0))
 
+@check_pdf_is_montecarlo
 def obs_pdf_correlations(pdf, results, xplotting_grid):
     """Return the correlations between each point in a dataset and the PDF
     values on a grid of (x,f) points in a format similar to `xplotting_grid`."""
     _, th = results
     # Wrap the result in a standard Stats class
     # since the result is (npoints, flavours, ndata) and has nothing to do with the PDF replicas
-    corrs = Stats(_basic_obs_pdf_correlation(xplotting_grid.grid_values, th.stats))
+    pdf_val = xplotting_grid.grid_values.error_members()
+    obs_val = th.error_members
+    corrs = Stats(_basic_obs_pdf_correlation(pdf_val, obs_val))
     return xplotting_grid.copy_grid(grid_values=corrs)
 
 
 corrpair_results = collect("results", ["corrpair"])
 corrpair_datasets = collect("dataset", ["corrpair"])
 
-
+@check_pdf_is_montecarlo
 def obs_obs_correlations(pdf, corrpair_results):
     """Return the theoretical correlation matrix between a pair of observables."""
     (_, th1), (_, th2) = corrpair_results
-    return _basic_obs_obs_correlation(th1.stats, th2.stats)
+    return _basic_obs_obs_correlation(th1.error_members, th2.error_members)

From 7e28ee25efa68fcaccef018b13dd9ae62f70d9d2 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Wed, 2 Mar 2022 11:30:44 +0100
Subject: [PATCH 08/13] enforce grid_values being stats

---
 validphys2/src/validphys/pdfgrids.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py
index 97596d2602..53eb6a16db 100644
--- a/validphys2/src/validphys/pdfgrids.py
+++ b/validphys2/src/validphys/pdfgrids.py
@@ -55,6 +55,11 @@ class XPlottingGrid:
     grid_values: Stats
     scale: str
 
+    def __post_init__(self):
+        """Enforce grid_values being a Stats instance"""
+        if not isinstance(self.grid_values, Stats):
+            raise ValueError("`XPlottingGrid` grid_values can only be instances of `Stats`")
+
     def select_flavour(self, flindex):
         """Return a new grid for one single flavour"""
         if isinstance(flindex, str):
@@ -72,11 +77,8 @@ def mask_replicas(self, mask):
         gv = self.grid_values.__class__(new_grid)
         return dataclasses.replace(self, grid_values=gv)
 
-    def copy_grid(self, grid_values=None):
+    def copy_grid(self, grid_values):
         """Create a copy of the grid with potentially a different set of values"""
-        if not isinstance(grid_values, Stats):
-            log.warning("XPlottingGrid being called with a numpy grid, should be using Stats instead!")
-            grid_values = self.grid_values.__class__(grid_values)
         return dataclasses.replace(self, grid_values=grid_values)
 
 

From 7081d4d81868b90c3f20d68c6fada9ae814f1d7d Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Wed, 2 Mar 2022 12:53:08 +0100
Subject: [PATCH 09/13] remove the special case of mask_replicas

---
 validphys2/src/validphys/deltachi2.py | 9 +++++++--
 validphys2/src/validphys/pdfgrids.py  | 6 ------
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/validphys2/src/validphys/deltachi2.py b/validphys2/src/validphys/deltachi2.py
index 2ba7297458..c1e8aff2fd 100644
--- a/validphys2/src/validphys/deltachi2.py
+++ b/validphys2/src/validphys/deltachi2.py
@@ -144,8 +144,13 @@ def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid):
     pos_mask = np.append(True, positive_eigenvalue_mask)
     neg_mask = np.append(True, ~positive_eigenvalue_mask)
 
-    pos_xplotting_grid = xplotting_grid.mask_replicas(pos_mask)
-    neg_xplotting_grid = xplotting_grid.mask_replicas(neg_mask)
+    pos_grid = xplotting_grid.grid_values.data[pos_mask]
+    neg_grid = xplotting_grid.grid_values.data[neg_mask]
+
+    # Wrap everything back into the original stats class
+    stats_class = xplotting_grid.grid_values.__class__
+    pos_xplotting_grid = xplotting_grid.copy_grid(stats_class(pos_grid))
+    neg_xplotting_grid = xplotting_grid.copy_grid(stats_class(neg_grid))
 
     return [xplotting_grid, pos_xplotting_grid, neg_xplotting_grid]
 
diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py
index 53eb6a16db..b89d98e4eb 100644
--- a/validphys2/src/validphys/pdfgrids.py
+++ b/validphys2/src/validphys/pdfgrids.py
@@ -71,12 +71,6 @@ def select_flavour(self, flindex):
         gv = self.grid_values.__class__(new_grid)
         return dataclasses.replace(self, grid_values=gv, flavours=[flstr])
 
-    def mask_replicas(self, mask):
-        """Return a copy of XPlottingGrid with the mask applied to the replicas"""
-        new_grid = self.grid_values.data[mask]
-        gv = self.grid_values.__class__(new_grid)
-        return dataclasses.replace(self, grid_values=gv)
-
     def copy_grid(self, grid_values):
         """Create a copy of the grid with potentially a different set of values"""
         return dataclasses.replace(self, grid_values=grid_values)

From 08d347e7f3452798ab1636430171c24fcdb112a0 Mon Sep 17 00:00:00 2001
From: "Juan M. Cruz-Martinez" <juacrumar@lairen.eu>
Date: Thu, 3 Mar 2022 18:17:11 +0100
Subject: [PATCH 10/13] Apply suggestions from code review

Co-authored-by: Roy Stegeman <roystegeman@live.nl>
---
 n3fit/src/n3fit/io/writer.py         | 1 -
 validphys2/src/validphys/sumrules.py | 1 -
 2 files changed, 2 deletions(-)

diff --git a/n3fit/src/n3fit/io/writer.py b/n3fit/src/n3fit/io/writer.py
index 7d0a5fb7e7..a988bf25a8 100644
--- a/n3fit/src/n3fit/io/writer.py
+++ b/n3fit/src/n3fit/io/writer.py
@@ -12,7 +12,6 @@
 import validphys
 import n3fit
 from n3fit import vpinterface
-from n3fit.vpinterface import integrability_numbers, compute_arclength
 
 
 class WriterWrapper:
diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py
index 05de4d90c6..731c338ec7 100644
--- a/validphys2/src/validphys/sumrules.py
+++ b/validphys2/src/validphys/sumrules.py
@@ -137,7 +137,6 @@ def _integral(rule_f, pdf_member, Q, config=None):
 def _sum_rules(rules_dict, lpdf, Q):
     """Compute a SumRulesGrid from the loaded PDF, at Q"""
     return {k: [_integral(r, m, Q) for m in lpdf.members] for k,r in rules_dict.items()}
-    # I very much doubt that looping over PDF members is the slow piece here...
 
 
 @check_positive('Q')

From 9199df9cccbe2739d88a8e46071b95ddd5ac8449 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Thu, 3 Mar 2022 20:53:41 +0100
Subject: [PATCH 11/13] remove grid_values_index

---
 validphys2/src/validphys/arclength.py   |  2 --
 validphys2/src/validphys/convolution.py | 14 +++++---------
 validphys2/src/validphys/core.py        | 19 -------------------
 validphys2/src/validphys/gridvalues.py  |  6 ------
 validphys2/src/validphys/lhapdfset.py   | 10 +++-------
 5 files changed, 8 insertions(+), 43 deletions(-)

diff --git a/validphys2/src/validphys/arclength.py b/validphys2/src/validphys/arclength.py
index f380d9fb36..8a27efba22 100644
--- a/validphys2/src/validphys/arclength.py
+++ b/validphys2/src/validphys/arclength.py
@@ -41,8 +41,6 @@ def arc_lengths(
     npoints = 199  # 200 intervals
     seg_min = [1e-6, 1e-4, 1e-2]
     seg_max = [1e-4, 1e-2, 1.0]
-    # Note that we discard replica 0 for MC PDFs, see core.PDF.grid_values_index
-    # for more details
     res = np.zeros((pdf.get_members(), len(flavours)))
     # Integrate the separate segments
     for a, b in zip(seg_min, seg_max):
diff --git a/validphys2/src/validphys/convolution.py b/validphys2/src/validphys/convolution.py
index f35e872fda..5203cb8e77 100644
--- a/validphys2/src/validphys/convolution.py
+++ b/validphys2/src/validphys/convolution.py
@@ -147,9 +147,7 @@ def predictions(dataset, pdf):
         A dataframe corresponding to the hadronic prediction for each data
         point for the PDF members. The index of the dataframe corresponds to
         the selected data points, based on the dataset :ref:`cuts <filters>`. The
-        columns correspond to the selected PDF members in the LHAPDF set, which
-        depend on the PDF error type (see
-        :py:meth:`validphys.core.PDF.grid_values_index`)
+        columns correspond to the selected PDF members in the LHAPDF set.
 
     Examples
     --------
@@ -219,9 +217,7 @@ def fk_predictions(loaded_fk, pdf):
         point for the PDF members. The index of the dataframe corresponds to
         the selected data points (use
         :py:meth:`validphys.coredata.FKTableData.with_cuts` to filter out
-        points). The columns correspond to the selected PDF members in the
-        LHAPDF set, which depend on the PDF error type (see
-        :py:meth:`validphys.core.PDF.grid_values_index`)
+        points). The columns correspond to the selected PDF members in the LHAPDF set.
 
     Notes
     -----
@@ -351,7 +347,7 @@ def hadron_predictions(loaded_fk, pdf):
     """Implementation of :py:func:`fk_predictions` for hadronic observables."""
     gv = functools.partial(evolution.grid_values, pdf=pdf)
     res = _gv_hadron_predictions(loaded_fk, gv)
-    res.columns = pdf.grid_values_index
+    res.columns = range(pdf.get_members())
     return res
 
 
@@ -380,7 +376,7 @@ def gv2(*args, **kwargs):
         return 2 * replica_values - central_value
 
     res = _gv_hadron_predictions(loaded_fk, gv1, gv2)
-    res.columns = pdf.grid_values_index
+    res.columns = range(pdf.get_members())
     return res
 
 
@@ -388,7 +384,7 @@ def dis_predictions(loaded_fk, pdf):
     """Implementation of :py:func:`fk_predictions` for DIS observables."""
     gv = functools.partial(evolution.grid_values, pdf=pdf)
     res = _gv_dis_predictions(loaded_fk, gv)
-    res.columns = pdf.grid_values_index
+    res.columns = range(pdf.get_members())
     return res
 
 
diff --git a/validphys2/src/validphys/core.py b/validphys2/src/validphys/core.py
index adaa264fd9..8681856129 100644
--- a/validphys2/src/validphys/core.py
+++ b/validphys2/src/validphys/core.py
@@ -226,25 +226,6 @@ def legacy_load(self):
 
         return libNNPDF_LHAPDFSet(self.name, et)
 
-    @property
-    def grid_values_index(self):
-        """A range object describing which members are selected in a
-        ``pdf.load().grid_values`` operation.  This is ``range(1,
-        len(pdf))`` for Monte Carlo sets, because replica 0 is not selected
-        and ``range(0, len(pdf))`` for hessian sets.
-
-
-        Returns
-        -------
-        index : range
-            A range object describing the proper indexing.
-
-        Notes
-        -----
-        The range object can be used efficiently as a Pandas index.
-        """
-        return range(0, len(self))
-
     def get_members(self):
         """Return the number of members selected in ``pdf.load().grid_values``
         """
diff --git a/validphys2/src/validphys/gridvalues.py b/validphys2/src/validphys/gridvalues.py
index 3b77f8bd3a..d038faa7f4 100644
--- a/validphys2/src/validphys/gridvalues.py
+++ b/validphys2/src/validphys/gridvalues.py
@@ -75,12 +75,6 @@ def grid_values(pdf:PDF, flmat, xmat, qmat):
 
         grid_values[replica][flavour][x][Q]
 
-    Notes
-    ----
-    This uses libnnpdf, and therefore follows the convention to throw away
-    replica 0 for Monte Carlo ensembles (so index 0 corresponds to replica 1).
-    Use ``pdf.grid_values_index`` to index the result properly.
-
     See Also
     --------
     :py:meth:`validphys.pdfbases.Basis.grid_values` offers a higher level
diff --git a/validphys2/src/validphys/lhapdfset.py b/validphys2/src/validphys/lhapdfset.py
index 8b52a24c90..ed087e27db 100644
--- a/validphys2/src/validphys/lhapdfset.py
+++ b/validphys2/src/validphys/lhapdfset.py
@@ -7,15 +7,13 @@
 
     The ``.members`` and ``.central_member`` of the ``LHAPDFSet`` are
     LHAPDF objects (the typical output from ``mkPDFs``) and can be used normally.
-    For MC PDFs the ``central_member`` is the average of all replicas (members 1-100)
-    while for Hessian PDFfs the ``central_member`` is also ``members[0]``
 
     Examples
     --------
     >>> from validphys.lhapdfset import LHAPDFSet
     >>> pdf = LHAPDFSet("NNPDF40_nnlo_as_01180", "replicas")
     >>> len(pdf.members)
-    100
+    101
     >>> pdf.central_member.alphasQ(91.19)
     0.11800
     >>> pdf.members[0].xfxQ2(0.5, 15625)
@@ -68,10 +66,8 @@ def n_members(self):
 
     @property
     def members(self):
-        """Return the members of the set, this depends on the error type:
-        t0: returns only member 0
-        MC: skip member 0
-        Hessian: return all
+        """Return the members of the set
+        the special error type t0 returns only member 0
         """
         if self.is_t0:
             return self._lhapdf_set[0:1]

From 2f0bc16e0f678b75bb05f583fb78d61a38be4712 Mon Sep 17 00:00:00 2001
From: juacrumar <juacrumar@lairen.eu>
Date: Thu, 3 Mar 2022 21:11:50 +0100
Subject: [PATCH 12/13] add a warning for calls with different Q

---
 n3fit/src/n3fit/performfit.py  |  5 +++--
 n3fit/src/n3fit/vpinterface.py | 11 ++++++++---
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/n3fit/src/n3fit/performfit.py b/n3fit/src/n3fit/performfit.py
index c2acd8a712..ff3a846196 100644
--- a/n3fit/src/n3fit/performfit.py
+++ b/n3fit/src/n3fit/performfit.py
@@ -264,14 +264,15 @@ def performfit(
             replica_path_set = replica_path / f"replica_{replica_number}"
 
             # Create a pdf instance
-            pdf_instance = N3PDF(pdf_model, fit_basis=basis)
+            q0 = theoryid.get_description().get("Q0")
+            pdf_instance = N3PDF(pdf_model, fit_basis=basis, Q=q0)
 
             # Generate the writer wrapper
             writer_wrapper = WriterWrapper(
                 replica_number,
                 pdf_instance,
                 stopping_object,
-                theoryid.get_description().get("Q0") ** 2,
+                q0**2,
                 final_time,
             )
 
diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py
index 6f58b1cad6..55e0665c0d 100644
--- a/n3fit/src/n3fit/vpinterface.py
+++ b/n3fit/src/n3fit/vpinterface.py
@@ -61,16 +61,21 @@ def central_value(self):
 class N3LHAPDFSet(LHAPDFSet):
     """Extension of LHAPDFSet using n3fit models"""
 
-    def __init__(self, name, pdf_models):
+    def __init__(self, name, pdf_models, Q=1.65):
         log.debug("Creating LHAPDF-like n3fit PDF")
         self._error_type = "replicas"
         self._name = name
         self._lhapdf_set = pdf_models
         self._flavors = None
+        self._fitting_q = Q
         self.basis = check_basis("evolution", EVOL_LIST)["basis"]
 
     def xfxQ(self, x, Q, n, fl):
         """Return the value of the PDF member for the given value in x"""
+        if Q != self._fitting_q:
+            log.warning(
+                "Querying N3LHAPDFSet at a value of Q=%f different from %f", Q, self._fitting_q
+            )
         return self.grid_values([fl], [x]).squeeze()[n]
 
     def __call__(self, xarr, flavours=None, replica=None):
@@ -168,7 +173,7 @@ class N3PDF(PDF):
             name of the N3PDF object
     """
 
-    def __init__(self, pdf_models, fit_basis=None, name="n3fit"):
+    def __init__(self, pdf_models, fit_basis=None, name="n3fit", Q=1.65):
         self.fit_basis = fit_basis
 
         if isinstance(pdf_models, Iterable):
@@ -178,7 +183,7 @@ def __init__(self, pdf_models, fit_basis=None, name="n3fit"):
 
         super().__init__(name)
         self._stats_class = N3Stats
-        self._lhapdf_set = N3LHAPDFSet(self.name, self._models)
+        self._lhapdf_set = N3LHAPDFSet(self.name, self._models, Q=Q)
         # Since there is no info file, create a fake `_info` dictionary
         self._info = {"ErrorType": "replicas", "NumMembers": len(self._models)}
 

From bf14730b94b4bc5f433bbdca989a15f12c8f242f Mon Sep 17 00:00:00 2001
From: "Juan M. Cruz-Martinez" <juacrumar@lairen.eu>
Date: Wed, 9 Mar 2022 14:35:58 +0100
Subject: [PATCH 13/13] Update validphys2/src/validphys/sumrules.py

Co-authored-by: Zaharid <zk261@cam.ac.uk>
---
 validphys2/src/validphys/sumrules.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/validphys2/src/validphys/sumrules.py b/validphys2/src/validphys/sumrules.py
index 731c338ec7..5196feea16 100644
--- a/validphys2/src/validphys/sumrules.py
+++ b/validphys2/src/validphys/sumrules.py
@@ -22,8 +22,8 @@
 
 
 def _momentum_sum_rule_integrand(x, lpdf, Q):
-    res = lpdf.xfxQ(x, Q)
-    return sum([res[f] for f in lpdf.flavors()])
+    xqvals = lpdf.xfxQ(x, Q)
+    return sum([xqvals[f] for f in lpdf.flavors()])
 
 def _make_momentum_fraction_integrand(fldict):
     """Make a suitable integrand function, which takes x to be integrated over