change: Handle uncertainties as None rather than 1s

news/profileparser-uncertainty-none.rst

@@ -0,0 +1,26 @@
+**Added:**
+
+* <news item>
+
+**Changed:**
+
+* Change ``ProfileParser.parse_file`` to set unavailable ``dx``/``dy`` uncertainties to ``None`` instead of ``0``, matching the format-specific parsers.
+* Change ``Profile.set_observed_profile`` to leave ``dyobs`` as ``None`` when no uncertainties are observed, instead of silently replacing them with an array of ones. The calculated ``dy`` still defaults to 1 at every calculation point, so unweighted fits refine exactly as before.
+* Change ``Profile._validate`` to accept a ``None`` ``dyobs``, since observed uncertainties are optional. ``x``, ``y``, ``dy``, ``xobs`` and ``yobs`` are still required.
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* Fix ``ProfileParser.parse_file`` producing zero-valued uncertainties for absent or invalid ``dx``/``dy`` columns, which gave non-finite residuals for any file without usable uncertainties instead of falling back to an unweighted fit.
+* Fix ``Profile.dyobs`` reporting an array of ones for data that carries no uncertainties, which made an unweighted profile indistinguishable from one whose uncertainties were genuinely all 1.
+
+**Security:**
+
+* <news item>

src/diffpy/srfit/fitbase/profile.py

@@ -186,8 +186,9 @@ def set_observed_profile(self, xobs, yobs, dyobs=None):
             Numpy array of the observed signal.
         dyobs
             Numpy array of the uncertainty in the observed signal. If
-            `dyobs` is None (default), it will be set to 1 at each
-            observed `xobs`.
+            `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.
 
 
         Raises
@@ -206,7 +207,7 @@ def set_observed_profile(self, xobs, yobs, dyobs=None):
         self._yobs = numpy.asarray(yobs, dtype=float)
 
         if dyobs is None:
-            self._dyobs = numpy.ones_like(xobs)
+            self._dyobs = None
         else:
             self._dyobs = numpy.asarray(dyobs, dtype=float)
 
@@ -328,7 +329,11 @@ 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]
+            self.dy = (
+                self.dyobs[indices]
+                if self.dyobs is not None
+                else numpy.ones_like(self.x)
+            )
         else:
             x1 = numpy.arange(lo, hi + epshi, step)
             self.set_calculation_points(x1)
@@ -373,6 +378,8 @@ def set_calculation_points(self, x):
                 # 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)
 
         return
 
@@ -477,8 +484,9 @@ def _flush(self, other):
     def _validate(self):
         """Validate my state.
 
-        This validates that x, y, dy, xobx, yobs and dyobs are not None.
-        This validates that x, y, and dy are the same length.
+        This validates that x, y, dy, xobs and yobs are not None. dyobs
+        may be None, since observed uncertainties are optional. This
+        validates that x, y, and dy are the same length.
 
         Raises SrFitError if validation fails.
         """
@@ -490,7 +498,6 @@ def _validate(self):
                 self.dy,
                 self.xobs,
                 self.yobs,
-                self.dyobs,
             ]
         )
         if datanotset:

src/diffpy/srfit/fitbase/profileparser.py

@@ -94,11 +94,11 @@ class ProfileParser(object):
             from the file.
         dx : np.ndarray
             The uncertainties associated with x
-            read from the file. This is 0 if the
+            read from the file. This is None if the
             uncertainty cannot be read.
         dy : np.ndarray
             The uncertainties associated with y
-            read from the file. This is 0 if the
+            read from the file. This is None if the
             uncertainty cannot be read.
     _x : np.ndarray
         Independent variable from the chosen bank
@@ -197,13 +197,13 @@ def parse_file(self, filename, column_format=None):
         """Parse a data file to extract data and metadata, with
         automatic handling of uncertainties.
 
-        - For files with 2 columns: assumes (x, y) and sets dx, dy to 0.
-        - For files with 3 columns: assumes (x, y, dy) and sets dx to 0.
+        - For files with 2 columns: assumes (x, y) and sets dx, dy to None.
+        - For files with 3 columns: assumes (x, y, dy) and sets dx to None.
         - 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 0.
+        are positive and not NaN/Inf. Otherwise they are set to None.
 
         This wipes out the currently loaded data and selected bank number.
 
@@ -242,10 +242,8 @@ def parse_file(self, filename, column_format=None):
         # Extract required arrays
         x = columns["x"]
         y = columns["y"]
-        x_length = len(x)
-        y_length = len(y)
-        dx = self._validate_uncertainty(columns.get("dx"), x_length)
-        dy = self._validate_uncertainty(columns.get("dy"), y_length)
+        dx = self._validate_uncertainty(columns.get("dx"))
+        dy = self._validate_uncertainty(columns.get("dy"))
         # Store as single bank
         self._banks = [(x, y, dx, dy)]
         self._meta["nbanks"] = 1
@@ -306,10 +304,10 @@ def _map_column_labels_to_data(self, data, column_format):
         return columns
 
     @staticmethod
-    def _validate_uncertainty(data, length):
-        """Return the uncertainty data if valid, otherwise 0."""
+    def _validate_uncertainty(data):
+        """Return the uncertainty data if valid, otherwise None."""
         if data is None or not np.all(np.isfinite(data)) or np.any(data <= 0):
-            return np.zeros(length)
+            return None
         return data
 
     def get_num_banks(self):
@@ -395,8 +393,8 @@ def get_data(self, index=None):
 
         Returns
         ----------
-        This returns (x, y, dx, dy) tuple for the bank. dx is 0 if it cannot
-        be determined from the data format.
+        This returns (x, y, dx, dy) tuple for the bank. dx is None if it
+        cannot be determined from the data format.
         """
         self.select_bank(index)
 

src/diffpy/srfit/pdf/pdfparser.py

@@ -63,11 +63,11 @@ class PDFParser(ProfileParser):
             from the file.
         dx
             A numpy array containing the uncertainty in x
-            read from the file. This is 0 if the
+            read from the file. This is None if the
             uncertainty cannot be read.
         dy
             A numpy array containing the uncertainty read
-            from the file. This is 0 if the uncertainty
+            from the file. This is None if the uncertainty
             cannot be read.
     _x
         Independent variable from the chosen bank
@@ -128,7 +128,8 @@ def parseString(self, patstring):
         """Parse a string and set the _x, _y, _dx, _dy and _meta
         variables.
 
-        When _dx or _dy cannot be obtained in the data format it is set to 0.
+        When _dx or _dy cannot be obtained in the data format it is set to
+        None.
 
         This wipes out the currently loaded data and selected bank number.
 

tests/conftest.py

@@ -219,6 +219,43 @@ def build_recipe_two_contributions():
     return recipe
 
 
+@pytest.fixture(scope="session")
+def as_list():
+    def _as_list(values):
+        """Unavailable uncertainties are None rather than an array."""
+        return None if values is None else values.tolist()
+
+    return _as_list
+
+
+@pytest.fixture(scope="function")
+def build_recipe_with_uncertainty():
+    """Helper to build a sine recipe whose profile carries the observed
+    uncertainties returned by ``make_dyobs``.
+
+    The observed profile is a noiseless sine, so a refinement recovers
+    ``A=1``, ``k=1`` and ``c=0`` regardless of how it is weighted.
+    """
+
+    def _build_recipe(make_dyobs):
+        xobs = linspace(0, pi, 21)
+        yobs = sin(xobs)
+        profile = Profile()
+        profile.set_observed_profile(xobs, yobs, make_dyobs(xobs))
+        contribution = FitContribution("c1")
+        contribution.set_profile(profile)
+        contribution.set_equation("A*sin(k*x + c)")
+        recipe = FitRecipe()
+        recipe.fithooks[0].verbose = 0
+        recipe.add_contribution(contribution)
+        recipe.add_variable(contribution.A, 0.8)
+        recipe.add_variable(contribution.k, 1.2)
+        recipe.add_variable(contribution.c, 0.1)
+        return recipe, profile
+
+    return _build_recipe
+
+
 @pytest.fixture
 def temp_data_files(tmp_path):
     """

tests/test_fitrecipe.py

@@ -20,7 +20,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
-from numpy import array_equal, dot, linspace, pi, sin
+from numpy import array_equal, dot, linspace, ones_like, pi, sin
 from scipy.optimize import leastsq
 
 from diffpy.srfit.fitbase import FitResults, ProfileParser
@@ -516,7 +516,7 @@ def test_initialize_recipe_from_recipe(build_recipes_one_contribution):
     assert sorted(list(expected_values)) == sorted(list(actual_values))
 
 
-def test_initialize_recipe_from_recipe_bad(build_recipe_two_contributions):
+def test_initialize_recipe_from_recipe_bad():
     # Case: User tries to initialize a FitRecipe from a non recipe object
     # expected: raised ValueError with message
     recipe_bad = 12345  # not a FitRecipe object
@@ -1060,5 +1060,86 @@ def test_plot_recipe_reset_all_defaults(build_recipes_one_contribution):
     assert actual_legend == expected_legend
 
 
+# Observed uncertainties are optional. A profile loaded without an
+# uncertainty column keeps dyobs as None and is refined unweighted, with
+# dy falling back to one at every calculation point. The cases below
+# refine the same noiseless sine profile so that the known solution
+# (A=1, k=1, c=0) is recovered no matter how the fit is weighted.
+
+
+# make_input_dyobs returns a function that evaluates dyobs on the array
+# of xobs for insertion into the build_recipe_with_uncertainty fixture.
+@pytest.mark.parametrize(
+    "make_input_dyobs, expected_dyobs_is_set",
+    [
+        # C1: No uncertainties are observed, as for a file with no
+        # uncertainty column.
+        # Expected: The refinement converges and dyobs stays None.
+        (lambda xobs: None, False),
+        # C2: Uniform uncertainties of one are observed explicitly.
+        # Expected: The refinement converges to the same values as C1,
+        # since an unweighted fit is weighted by one everywhere.
+        (lambda xobs: ones_like(xobs), True),
+        # C3: Uncertainties vary across the profile, so the refinement
+        # is weighted point by point.
+        # Expected: The refinement still converges to the known values.
+        (lambda xobs: linspace(0.1, 1.0, len(xobs)), True),
+    ],
+)
+def test_refine_with_and_without_uncertainty(
+    build_recipe_with_uncertainty, make_input_dyobs, expected_dyobs_is_set
+):
+    recipe, profile = build_recipe_with_uncertainty(make_input_dyobs)
+    optimize_recipe(recipe)
+    actual_dyobs_is_set = profile.dyobs is not None
+    actual_names = recipe.get_names()
+    actual_values = recipe.get_values()
+    expected_names = ["A", "k", "c"]
+    expected_values = [1.0, 1.0, 0.0]
+    assert actual_dyobs_is_set == expected_dyobs_is_set
+    assert actual_names == expected_names
+    assert actual_values == pytest.approx(expected_values, abs=1e-5)
+
+
+# The residual that is optimized is (ycalc - y)/dy, so dy is what
+# actually weights a refinement. These cases check that a profile with
+# no observed uncertainties is weighted identically to one whose
+# uncertainties are all one, and that observed uncertainties are carried
+# through to the residual unchanged.
+
+
+# make_input_dyobs returns a function that evaluates dyobs on the array
+# of xobs for insertion into the build_recipe_with_uncertainty fixture.
+@pytest.mark.parametrize(
+    "make_input_dyobs, make_expected_dy",
+    [
+        # C1: No uncertainties are observed.
+        # Expected: dy falls back to one everywhere, so the residual is
+        # unweighted.
+        (lambda xobs: None, lambda xobs: ones_like(xobs)),
+        # C2: Uniform uncertainties of one are observed explicitly.
+        # Expected: dy is one everywhere, matching C1.
+        (lambda xobs: ones_like(xobs), lambda xobs: ones_like(xobs)),
+        # C3: Uncertainties vary across the profile.
+        # Expected: dy keeps the observed values, so the residual is
+        # weighted point by point.
+        (
+            lambda xobs: linspace(0.1, 1.0, len(xobs)),
+            lambda xobs: linspace(0.1, 1.0, len(xobs)),
+        ),
+    ],
+)
+def test_residual_is_weighted_by_uncertainty(
+    build_recipe_with_uncertainty, make_input_dyobs, make_expected_dy
+):
+    recipe, profile = build_recipe_with_uncertainty(make_input_dyobs)
+    actual_residual = recipe.residual(recipe.values)
+    actual_dy = profile.dy
+    expected_dy = make_expected_dy(profile.xobs)
+    expected_residual = (profile.ycalc - profile.y) / expected_dy
+    assert actual_dy == pytest.approx(expected_dy)
+    assert actual_residual == pytest.approx(expected_residual)
+
+
 if __name__ == "__main__":
     unittest.main()

tests/test_pdf.py

@@ -139,7 +139,8 @@ def testParser2(datafile):
     res = numpy.dot(diff, diff)
     assert 0 == pytest.approx(res)
 
-    assert dx.tolist() == [0] * len(dx)
+    # si-q27r60-xray.gr has a negative dx column, so it is invalid.
+    assert dx is None
     return