fix: give points with invalid uncertainties no weight instead of dropping all uncertainties or making the residual infinite

news/uncertainty-checks.rst

@@ -0,0 +1,24 @@
+**Added:**
+
+* <news item>
+
+**Changed:**
+
+* ``ProfileParser`` reads the ``dy`` column as it is in the file, instead of dropping the whole column when one uncertainty is zero, negative, infinite or not a number.
+* ``FitResults`` counts only the points with weight in Reduced Chi2, and says how many points of each contribution have no weight.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* A point whose uncertainty is zero, negative, infinite or not a number gets no weight in a fit, instead of making the residual infinite or undefined. When no point in the calculation range has a valid uncertainty, the profile is weighted as if it had none.
+
+**Security:**
+
+* <news item>

src/diffpy/srfit/fitbase/fitcontribution.py

@@ -474,6 +474,9 @@ def residual(self):
         The residual equation can be changed with the
         ``set_residual_equation`` method.
 
+        A point with no weight, whose ``dy`` is infinite, has a ``chiv``
+        residual of zero.
+
         Returns
         -------
         numpy.ndarray

src/diffpy/srfit/fitbase/fitrecipe.py

@@ -37,7 +37,7 @@
 from collections import OrderedDict
 from pathlib import Path
 
-from numpy import array, concatenate, dot, sqrt
+from numpy import array, concatenate, count_nonzero, dot, isinf, sqrt
 
 import diffpy.srfit.util.inpututils as utils
 from diffpy.srfit.fitbase.fithook import PrintFitHook
@@ -532,8 +532,11 @@ def residual(self, p=[]):
             ]
         )
 
-        # Calculate the point-average chi^2
-        w = dot(chiv, chiv) / len(chiv)
+        # Calculate the point-average chi^2 over the points with weight
+        unweighted_count = 0
+        for con in self._contributions.values():
+            unweighted_count += count_nonzero(isinf(con.profile.dy))
+        w = dot(chiv, chiv) / (len(chiv) - unweighted_count)
         # Now we must append the restraints
         penalties = [sqrt(res.penalty(w)) for res in self._restraintlist]
         chiv = concatenate([chiv, penalties])

src/diffpy/srfit/fitbase/fitresults.py

@@ -137,7 +137,8 @@ class FitResults(object):
         The cumulative chi-squared as a function of data index.
 
     rchi2 : float
-        The reduced chi-squared of the fit.
+        The reduced chi-squared of the fit. Only the points with weight
+        in the fit are counted.
 
     rw : float
         The weighted R-factor of the fit.
@@ -255,7 +256,10 @@ def update(self):
         self._calculate_metrics()
 
         # Calculate the restraints penalty
-        w = self.chi2 / len(res)
+        unweighted_count = 0
+        for con in self.conresults.values():
+            unweighted_count += numpy.count_nonzero(numpy.isinf(con.dy))
+        w = self.chi2 / (len(res) - unweighted_count)
         self.penalty = sum([r.penalty(w) for r in recipe._restraintlist])
 
         return
@@ -356,7 +360,7 @@ def _calculate_metrics(self):
             c2last = cumchi2[-1:].sum()
             cumchi2 = numpy.concatenate([cumchi2, c2last + cc2w])
             yw2tot += con.weight * (con.chi2 / con.rw**2)
-            numpoints += len(con.x)
+            numpoints += numpy.count_nonzero(numpy.isfinite(con.dy))
 
         chi2 = cumchi2[-1:].sum()
         cumrw = numpy.sqrt(cumchi2 / yw2tot)
@@ -440,6 +444,19 @@ def get_results_string(self, header="", footer="", update=False):
             )
             if err_msg not in self.messages:
                 self.messages.append(err_msg)
+        for con_name, con in self.conresults.items():
+            unweighted_count = numpy.count_nonzero(numpy.isinf(con.dy))
+            if unweighted_count == 0:
+                continue
+            unweighted_msg = (
+                f"{unweighted_count} of the {len(con.dy)} points of "
+                f"contribution '{con_name}' have no weight, because their "
+                "uncertainties are zero, negative, infinite or not a number. "
+                "To include them in the fit, please set valid uncertainties "
+                "with Profile.set_observed_profile."
+            )
+            if unweighted_msg not in self.messages:
+                self.messages.append(unweighted_msg)
 
         lines.extend(self.messages)
 

src/diffpy/srfit/fitbase/profile.py

@@ -87,7 +87,9 @@ class Profile(Observable, Validatable):
         property gets and sets the value of ``ypar``.
     dy : numpy.ndarray or None
         The uncertainty in the profile over the calculation range
-        (default None). This property gets and sets the value of
+        (default None). It is infinite at points with no weight in a fit,
+        and one everywhere when no point in the range has a valid
+        uncertainty. This property gets and sets the value of
         ``dypar``.
     ycalc : numpy.ndarray or None
         The calculated signal (default None). This property gets and
@@ -195,6 +197,11 @@ def set_observed_profile(self, xobs, yobs, dyobs=None):
             ``dyobs`` is None (default), ``dyobs`` stays None to
             indicate no uncertainty was observed, and the calculated
             ``dy`` will be set to 1 at each calculation point instead.
+            ``dyobs`` is kept as given. A point whose uncertainty is
+            zero, negative, infinite or not a number gets no weight in a
+            fit, so its calculated ``dy`` is infinite. If no point in the
+            calculation range has a valid uncertainty, ``dy`` is set to 1
+            at each calculation point, as when ``dyobs`` is None.
 
         Raises
         ------
@@ -336,11 +343,7 @@ def _isobs(a):
             indices = (lo - epslo <= self.xobs) & (self.xobs <= hi + epshi)
             self.x = self.xobs[indices]
             self.y = self.yobs[indices]
-            self.dy = (
-                self.dyobs[indices]
-                if self.dyobs is not None
-                else numpy.ones_like(self.x)
-            )
+            self.dy = self._compute_calculation_uncertainty()
         else:
             x1 = numpy.arange(lo, hi + epshi, step)
             self.set_calculation_points(x1)
@@ -362,7 +365,10 @@ def set_calculation_points(self, x):
         """Set the calculation points.
 
         This creates ``y`` and ``dy`` on the specified grid if
-        ``xobs``, ``yobs`` and ``dyobs`` exist.
+        ``xobs``, ``yobs`` and ``dyobs`` exist. On a grid that differs
+        from ``xobs``, ``dy`` is interpolated from the valid observed
+        uncertainties only, and a calculation point next to an observed
+        point with no weight gets no weight either.
 
         Parameters
         ----------
@@ -377,19 +383,31 @@ def set_calculation_points(self, x):
         self.x = x
         if self.yobs is not None:
             self.y = _rebin_array(self.yobs, self.xobs, self.x)
-        if self.dyobs is not None:
-            # work around for interpolation issue making some of these non-1
-            if (self.dyobs == 1).all():
-                self.dy = numpy.ones_like(self.x)
-            else:
-                # FIXME - This does not follow error propagation rules and it
-                # introduces (more) correlation between the data points.
-                self.dy = _rebin_array(self.dyobs, self.xobs, self.x)
-        elif self.yobs is not None:
-            self.dy = numpy.ones_like(self.x)
-
+            self.dy = self._compute_calculation_uncertainty()
         return
 
+    def _compute_calculation_uncertainty(self):
+        """Return dy on the calculation points, infinite where a point
+        has no weight, or ones when no point has a valid uncertainty."""
+        unweighted_dy = numpy.ones_like(self.x, dtype=float)
+        if self.dyobs is None:
+            return unweighted_dy
+        is_valid = numpy.isfinite(self.dyobs) & (self.dyobs > 0)
+        # A calculation point between an observed point with no weight and
+        # its neighbour interpolates to less than 1, so it has no weight.
+        is_weighted = (
+            _rebin_array(is_valid.astype(float), self.xobs, self.x) == 1
+        )
+        if not is_weighted.any():
+            return unweighted_dy
+        dy = numpy.full_like(self.x, numpy.inf, dtype=float)
+        # FIXME - This does not follow error propagation rules and it
+        # introduces (more) correlation between the data points.
+        dy[is_weighted] = numpy.interp(
+            self.x[is_weighted], self.xobs[is_valid], self.dyobs[is_valid]
+        )
+        return dy
+
     @deprecated(setCalculationPoints_dep_msg)
     def setCalculationPoints(self, x):
         """This function has been deprecated and will be removed in version

src/diffpy/srfit/fitbase/profileparser.py

@@ -100,8 +100,8 @@ class ProfileParser(object):
             uncertainty cannot be read.
         dy : np.ndarray
             The uncertainties associated with y
-            read from the file. This is None if the
-            uncertainty cannot be read.
+            read from the file, as they are in the file.
+            This is None if the file has no dy column.
     _x : np.ndarray
         Independent variable from the chosen bank
     _y : np.ndarray
@@ -187,8 +187,11 @@ def parse_file(
         - For files with 4 columns: assumes (x, y, dx, dy).
         - For other cases: ``column_format`` must be explicitly specified.
 
-        Uncertainty columns (dx, dy) are only considered valid if all values
-        are positive and not NaN/Inf. Otherwise they are set to None.
+        The dy column is read as it is, since a ``Profile`` gives the
+        points whose uncertainties are zero, negative, infinite or not a
+        number no weight in a fit. The dx column is only considered valid
+        if all values are positive and not NaN/Inf. Otherwise it is set to
+        None.
 
         This wipes out the currently loaded data and selected bank number.
 
@@ -302,7 +305,7 @@ def _parse_data(self, filename, column_format=None, **kwargs):
                 columns["x"],
                 columns["y"],
                 self._validate_uncertainty(columns.get("dx")),
-                self._validate_uncertainty(columns.get("dy")),
+                columns.get("dy"),
             ]
         )
 

tests/conftest.py

@@ -418,6 +418,32 @@ def parser_datafiles(tmp_path):
         )
     )
 
+    (tmp_path / "three_col_zero_dy.dat").write_text(
+        metadata
+        + "\n".join(
+            [
+                r"#L r($\AA$)  G($\AA^{-2}$) dG($\AA^{-2}$)",
+                "1.0 2.0 0.2",
+                "1.1 2.1 0.0",
+                "1.2 2.2 0.6",
+                "",
+            ]
+        )
+    )
+
+    (tmp_path / "four_col_zero_dx.gr").write_text(
+        metadata
+        + "\n".join(
+            [
+                r"#L r($\AA$)  G($\AA^{-2}$) dr($\AA$) dG($\AA^{-2}$)",
+                "1.0 2.0 0.0 0.2",
+                "1.1 2.1 0.0 0.4",
+                "1.2 2.2 0.0 0.6",
+                "",
+            ]
+        )
+    )
+
     (tmp_path / "one_col.gr").write_text(
         metadata
         + "\n".join(

tests/test_contribution.py

@@ -14,8 +14,10 @@
 ##############################################################################
 """Tests for the fitcontribution module."""
 
+import warnings
+
 import pytest
-from numpy import arange, array_equal, dot, sin
+from numpy import allclose, arange, array_equal, dot, nan, sin
 
 from diffpy.srfit import FitContribution, Profile, ProfileGenerator
 from diffpy.srfit.exceptions import SrFitError
@@ -275,6 +277,45 @@ def test_residual(
     assert actual_builders == expected_builders
 
 
+# A point whose uncertainty is zero or not a number has no weight, so its
+# chiv residual is zero. The profile is observed at x = 1, 2, 3, 4, 5 with
+# y = x, and the equation is "A*x" with A = 2, so the residual of every
+# point with an uncertainty of one is x.
+@pytest.mark.parametrize(
+    "input_dyobs, expected_residual",
+    [
+        # C1: A user loads data whose last two points have uncertainties of
+        # zero, as at the end of a detector bank.
+        # Expected: The residual is zero at those two points, with no
+        # warning.
+        ([1.0, 1.0, 1.0, 0.0, 0.0], [1.0, 2.0, 3.0, 0.0, 0.0]),
+        # C2: A user loads data with one uncertainty that is not a number,
+        # at x = 3.
+        # Expected: The residual is zero at x = 3, with no warning.
+        ([1.0, 1.0, nan, 1.0, 1.0], [1.0, 2.0, 0.0, 4.0, 5.0]),
+    ],
+)
+def test_residual_with_invalid_uncertainties(input_dyobs, expected_residual):
+    input_xobs = [1.0, 2.0, 3.0, 4.0, 5.0]
+    profile = Profile()
+    profile.set_observed_profile(input_xobs, input_xobs, input_dyobs)
+    contribution = FitContribution("test")
+    contribution.set_profile(profile)
+    contribution.set_equation("A*x")
+    contribution.A.set_value(2)
+    # Compute the residual, recording every warning to check that none is
+    # raised
+    with warnings.catch_warnings(record=True) as caught_warnings:
+        warnings.simplefilter("always")
+        actual_residual = contribution.residual()
+    actual_warning_messages = []
+    for caught_warning in caught_warnings:
+        actual_warning_messages.append(str(caught_warning.message))
+    expected_warning_messages = []
+    assert allclose(actual_residual, expected_residual)
+    assert actual_warning_messages == expected_warning_messages
+
+
 # The cached values of the profile and residual equations are cleared
 # whenever their result may change. The contribution has a profile observed
 # at x = 0, 0.5, ..., 9.5 with y = x, and a generator "I".

tests/test_fitrecipe.py

@@ -1302,3 +1302,38 @@ def test_residual_is_weighted_by_uncertainty(
     expected_residual = (profile.ycalc - profile.y) / expected_dy
     assert actual_dy == pytest.approx(expected_dy)
     assert actual_residual == pytest.approx(expected_residual)
+
+
+def test_residual_with_restraint_and_points_with_no_weight(
+    build_recipe_with_uncertainty,
+):
+    # Test the residual of a recipe with scaled soft bounds whose data have
+    # three points with no weight, compared with the same recipe on data
+    # without those points
+    # Expected: The residuals are the same, apart from the zeros of the
+    # points with no weight, because the scaled restraint is weighted by
+    # the average chi2 of the points with weight only.
+    # Build the recipe whose last three points have no weight
+    recipe_with_unweighted_points, _ = build_recipe_with_uncertainty(
+        [0.1] * 18 + [0.0] * 3
+    )
+    recipe_with_unweighted_points.add_soft_bounds(
+        "A", lower_bound=1.0, scaled=True
+    )
+    # Build the same recipe on the first 18 points only
+    xobs = linspace(0, pi, 21)[:18]
+    profile = Profile()
+    profile.set_observed_profile(xobs, sin(xobs), [0.1] * 18)
+    contribution = FitContribution("c1")
+    contribution.set_profile(profile)
+    contribution.set_equation("A*sin(k*x + c)")
+    recipe_without_points = FitRecipe()
+    recipe_without_points.add_contribution(contribution)
+    recipe_without_points.add_variable(contribution.A, 0.8)
+    recipe_without_points.add_variable(contribution.k, 1.2)
+    recipe_without_points.add_variable(contribution.c, 0.1)
+    recipe_without_points.add_soft_bounds("A", lower_bound=1.0, scaled=True)
+    actual_residual = recipe_with_unweighted_points.residual().tolist()
+    expected_residual = recipe_without_points.residual().tolist()
+    expected_residual[18:18] = [0.0, 0.0, 0.0]
+    assert actual_residual == pytest.approx(expected_residual)