diff --git a/news/profileparser-uncertainty-none.rst b/news/profileparser-uncertainty-none.rst new file mode 100644 index 00000000..f408ceb6 --- /dev/null +++ b/news/profileparser-uncertainty-none.rst @@ -0,0 +1,26 @@ +**Added:** + +* + +**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:** + +* + +**Removed:** + +* + +**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:** + +* diff --git a/src/diffpy/srfit/fitbase/profile.py b/src/diffpy/srfit/fitbase/profile.py index 18c8331f..3718e907 100644 --- a/src/diffpy/srfit/fitbase/profile.py +++ b/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: diff --git a/src/diffpy/srfit/fitbase/profileparser.py b/src/diffpy/srfit/fitbase/profileparser.py index 78dc3c27..8be242c0 100644 --- a/src/diffpy/srfit/fitbase/profileparser.py +++ b/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) diff --git a/src/diffpy/srfit/pdf/pdfparser.py b/src/diffpy/srfit/pdf/pdfparser.py index b45e6b08..885cff62 100644 --- a/src/diffpy/srfit/pdf/pdfparser.py +++ b/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. diff --git a/tests/conftest.py b/tests/conftest.py index c741e3cb..2f40a218 100644 --- a/tests/conftest.py +++ b/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): """ diff --git a/tests/test_fitrecipe.py b/tests/test_fitrecipe.py index 5eceb921..8b62f786 100644 --- a/tests/test_fitrecipe.py +++ b/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() diff --git a/tests/test_pdf.py b/tests/test_pdf.py index 315f6dff..aaeb0cd4 100644 --- a/tests/test_pdf.py +++ b/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 diff --git a/tests/test_profile.py b/tests/test_profile.py index ccbef6b2..54b47d20 100644 --- a/tests/test_profile.py +++ b/tests/test_profile.py @@ -22,7 +22,7 @@ from numpy import allclose, arange, array, array_equal, ones_like from diffpy.srfit.exceptions import SrFitError -from diffpy.srfit.fitbase import ProfileParser +from diffpy.srfit.fitbase import FitContribution, ProfileParser from diffpy.srfit.fitbase.profile import Profile @@ -68,7 +68,7 @@ def test_set_observed_profile(self): self.assertTrue(array_equal(x, prof.xobs)) self.assertTrue(array_equal(y, prof.yobs)) - self.assertTrue(array_equal(ones_like(prof.xobs), prof.dyobs)) + self.assertTrue(prof.dyobs is None) # Get the ranged profile to make sure its the same self.assertTrue(array_equal(x, prof.x)) @@ -104,7 +104,7 @@ def testSetObservedProfile(self): self.assertTrue(array_equal(x, prof.xobs)) self.assertTrue(array_equal(y, prof.yobs)) - self.assertTrue(array_equal(ones_like(prof.xobs), prof.dyobs)) + self.assertTrue(prof.dyobs is None) # Get the ranged profile to make sure its the same self.assertTrue(array_equal(x, prof.x)) @@ -307,19 +307,76 @@ def _test(p): return -def test_load_parsed_data(parser_datafiles): - """Test the load_parsed_data method.""" - prof = Profile() +# The parsed x, y and dy arrays are copied onto the observed profile. +# Uncertainties on x are dropped, since srfit treats the independent +# variable as having no uncertainty. +@pytest.mark.parametrize( + "input_filename, expected_xobs, expected_yobs, expected_dyobs", + [ + # C1: File has four columns, so uncertainties are present. + # Expected: xobs, yobs and dyobs are loaded and dx is dropped. + ( + "four_col.gr", + [1.0, 1.1, 1.2], + [2.0, 2.1, 2.2], + [0.2, 0.4, 0.6], + ), + # C3: File has three columns, so uncertainties are present. + # Expected: xobs, yobs and dyobs are loaded. + ( + "three_col.dat", + [1.0, 1.1, 1.2], + [2.0, 2.1, 2.2], + [0.2, 0.4, 0.6], + ), + # C2: File has two columns, so no uncertainties are present. + # Expected: xobs and yobs are loaded and dyobs stays None, + # marking the profile as unweighted. + ( + "two_col.txt", + [1.0, 1.1, 1.2], + [2.0, 2.1, 2.2], + None, + ), + ], +) +def test_load_parsed_data( + as_list, + parser_datafiles, + input_filename, + expected_xobs, + expected_yobs, + expected_dyobs, +): + """Load a parsed profile onto the observed arrays.""" parser = ProfileParser() - datafile = parser_datafiles / "four_col.gr" - parser.parse_file(datafile) + parser.parse_file(parser_datafiles / input_filename) + prof = Profile() prof.load_parsed_data(parser) - expected_xobs = [1.0, 1.1, 1.2] - expected_yobs = [2.0, 2.1, 2.2] - expected_dyobs = [0.2, 0.4, 0.6] - assert prof.xobs.tolist() == expected_xobs - assert prof.yobs.tolist() == expected_yobs - assert prof.dyobs.tolist() == expected_dyobs + actual_xobs = prof.xobs.tolist() + actual_yobs = prof.yobs.tolist() + # Unavailable uncertainties are None rather than an array. + actual_dyobs = as_list(prof.dyobs) + assert actual_xobs == expected_xobs + assert actual_yobs == expected_yobs + assert actual_dyobs == expected_dyobs + + +# Case: A profile is set without uncertainties, so dyobs is None, and a +# residual is computed from the equation y = A*x with A = 2. +# Expected: The residual is (A*x - y) with unit weights, since missing +# uncertainties default to ones rather than being used as divisors. +def test_residual_without_uncertainties(): + """Compute the residual for a profile with no uncertainties.""" + prof = Profile() + prof.set_observed_profile(array([1.0, 2.0, 3.0]), array([1.0, 2.0, 3.0])) + contribution = FitContribution("test") + contribution.set_profile(prof) + contribution.set_equation("A*x") + contribution.A.set_value(2) + actual_residual = contribution.residual().tolist() + expected_residual = [1.0, 2.0, 3.0] + assert actual_residual == expected_residual if __name__ == "__main__": diff --git a/tests/test_profileparser.py b/tests/test_profileparser.py index e1676b15..2e078561 100644 --- a/tests/test_profileparser.py +++ b/tests/test_profileparser.py @@ -17,7 +17,7 @@ # expected: x, y, dx, dy, and metadata are all read correctly # UC5: User loads file with dy and dx values containing NaN and inf values # expected: x, y, and metadata are all read correctly and dx and dy are set to -# 0 for all values +# None # UC6: User loads file with only one column # expected: ParseError is raised @@ -36,6 +36,7 @@ # duplicate values # expected: ParseError is raised + EXPECTED_META = { "wavelength": 0.1, "dataformat": "QA", @@ -81,7 +82,7 @@ None, [1.0, 1.1, 1.2], [2.0, 2.1, 2.2], - [0.0, 0.0, 0.0], + None, [0.2, 0.4, 0.6], ), # UC3: 2-column file (x, y) — dx and dy are missing @@ -91,8 +92,8 @@ None, [1.0, 1.1, 1.2], [2.0, 2.1, 2.2], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], + None, + None, ), # UC4: 4-column file in (x, dx, y, dy) order with explicit # column_format @@ -107,18 +108,19 @@ ), # UC5: 4-column file where dx/dy contain NaN and inf values # expected: x, y, and metadata are read correctly; dx and dy - # are set to 0 + # are set to None ( Path("four_col_nan_inf.gr"), None, [1.0, 1.1, 1.2], [2.0, 2.1, 2.2], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], + None, + None, ), ], ) def test_parse_file( + as_list, parser_datafiles, input_file, column_order, @@ -131,8 +133,8 @@ def test_parse_file( parser.parse_file(parser_datafiles / input_file, column_order) actual_x = parser._x.tolist() actual_y = parser._y.tolist() - actual_dx = parser._dx.tolist() - actual_dy = parser._dy.tolist() + actual_dx = as_list(parser._dx) + actual_dy = as_list(parser._dy) actual_metadata = parser._meta actual_metadata["filename"] = actual_metadata["filename"].split("/")[-1]