From 961cffd336b7d844cd8ef73678b4c895f64f7a98 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Fri, 30 Oct 2026 12:24:06 -0400 Subject: [PATCH 1/3] fix: raise on uncertainties that are zero or NaN in FitContribution and warn when profileparse drops dy --- news/uncertainty-checks.rst | 23 ++++++ src/diffpy/srfit/fitbase/fitcontribution.py | 23 +++++- src/diffpy/srfit/fitbase/profile.py | 3 + src/diffpy/srfit/fitbase/profileparser.py | 12 +++- tests/conftest.py | 26 +++++++ tests/test_contribution.py | 62 ++++++++++++++++ tests/test_profileparser.py | 79 ++++++++++++++++++--- 7 files changed, 215 insertions(+), 13 deletions(-) create mode 100644 news/uncertainty-checks.rst diff --git a/news/uncertainty-checks.rst b/news/uncertainty-checks.rst new file mode 100644 index 00000000..f4f217a9 --- /dev/null +++ b/news/uncertainty-checks.rst @@ -0,0 +1,23 @@ +**Added:** + +* + +**Changed:** + +* ``ProfileParser`` warns when it ignores a ``dy`` column whose uncertainties are zero, negative, infinite or not a number, instead of weighting every point equally without saying so. + +**Deprecated:** + +* + +**Removed:** + +* + +**Fixed:** + +* ``FitContribution.residual`` raises an ``SrFitError`` that counts the points whose uncertainties are zero or not a number, instead of returning an infinite or undefined residual, unless ``Profile.set_calculation_range`` excludes them. + +**Security:** + +* diff --git a/src/diffpy/srfit/fitbase/fitcontribution.py b/src/diffpy/srfit/fitbase/fitcontribution.py index 505b97ba..1e40a99e 100644 --- a/src/diffpy/srfit/fitbase/fitcontribution.py +++ b/src/diffpy/srfit/fitbase/fitcontribution.py @@ -24,6 +24,8 @@ __all__ = ["FitContribution"] +import numpy + from diffpy.srfit.exceptions import SrFitError from diffpy.srfit.fitbase.parameter import ParameterProxy from diffpy.srfit.fitbase.parameterset import ParameterSet @@ -478,12 +480,31 @@ def residual(self): ------- numpy.ndarray The array of residual values. + + Raises + ------ + SrFitError + If the residual isn't finite at points whose uncertainties + ``dy`` are zero or not a number. """ # Assign the calculated profile self.profile.ycalc = self._eq() # Note that equations only recompute when their inputs are modified, so # the following will not recompute the equation. - return self._reseq() + residual = self._reseq() + if not numpy.all(numpy.isfinite(residual)): + dy = self.profile.dy + invalid_count = numpy.count_nonzero((dy == 0) | numpy.isnan(dy)) + if invalid_count > 0: + raise SrFitError( + f"The residual of FitContribution '{self.name}' can't " + f"be computed at {invalid_count} of its {len(dy)} " + "points, because their uncertainties are zero or not a " + "number. Please exclude those points with " + "Profile.set_calculation_range, or set valid " + "uncertainties with Profile.set_observed_profile." + ) + return residual def evaluate(self): """Evaluate the contribution equation and update profile.ycalc. diff --git a/src/diffpy/srfit/fitbase/profile.py b/src/diffpy/srfit/fitbase/profile.py index 22452c68..b9583aeb 100644 --- a/src/diffpy/srfit/fitbase/profile.py +++ b/src/diffpy/srfit/fitbase/profile.py @@ -195,6 +195,9 @@ 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. + Points whose uncertainties are zero or not a number can't + weight a fit, so the residual raises an error unless + ``set_calculation_range`` excludes them. Raises ------ diff --git a/src/diffpy/srfit/fitbase/profileparser.py b/src/diffpy/srfit/fitbase/profileparser.py index 39823c58..d63024fb 100644 --- a/src/diffpy/srfit/fitbase/profileparser.py +++ b/src/diffpy/srfit/fitbase/profileparser.py @@ -23,6 +23,7 @@ """ import warnings +from pathlib import Path import numpy as np @@ -297,12 +298,21 @@ def _parse_data(self, filename, column_format=None, **kwargs): ) column_format = self._detect_column_format(data, column_format) columns = self._map_column_labels_to_data(data, column_format) + dy = self._validate_uncertainty(columns.get("dy")) + if "dy" in columns and dy is None: + warnings.warn( + f"The dy column of '{Path(filename).name}' has uncertainties " + "that are zero, negative, infinite or not a number, so all " + "of them are ignored and every point is weighted equally. " + "Please check the uncertainties in the file, or set them " + "with Profile.set_observed_profile." + ) self._banks.append( [ columns["x"], columns["y"], self._validate_uncertainty(columns.get("dx")), - self._validate_uncertainty(columns.get("dy")), + dy, ] ) diff --git a/tests/conftest.py b/tests/conftest.py index 1a5f4494..d36e3bb2 100644 --- a/tests/conftest.py +++ b/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( diff --git a/tests/test_contribution.py b/tests/test_contribution.py index 50d168fb..0d366127 100644 --- a/tests/test_contribution.py +++ b/tests/test_contribution.py @@ -14,6 +14,8 @@ ############################################################################## """Tests for the fitcontribution module.""" +import re + import pytest from numpy import arange, array_equal, dot, sin @@ -275,6 +277,66 @@ def test_residual( assert actual_builders == expected_builders +# The chiv residual divides by the uncertainties dy, so it can't be +# computed at a point whose uncertainty is zero or not a number. The +# profile is observed at the 20 points x = 0, 0.5, ..., 9.5 with y = x, and +# the equation is "A*x" with A = 2. +@pytest.mark.parametrize( + "input_dyobs, expected_invalid_count", + [ + # C1: A user loads data whose last point has an uncertainty of + # zero. + # Expected: An SrFitError says that 1 of the 20 points has an + # uncertainty that is zero or not a number, and suggests excluding + # it or setting valid uncertainties. + ([1.0] * 19 + [0.0], "1"), + # C2: A user loads data whose last three points have uncertainties + # of zero, as at the end of a detector bank. + # Expected: An SrFitError says that 3 of the 20 points do. + ([1.0] * 17 + [0.0] * 3, "3"), + # C3: A user loads data with one uncertainty that is not a number. + # Expected: An SrFitError says that 1 of the 20 points does. + ([1.0] * 10 + [float("nan")] + [1.0] * 9, "1"), + ], +) +def test_residual_with_invalid_uncertainties( + input_dyobs, expected_invalid_count +): + profile = Profile() + profile.set_observed_profile(XOBS, XOBS, input_dyobs) + contribution = FitContribution("test") + contribution.set_profile(profile) + contribution.set_equation("A*x") + contribution.A.set_value(2) + expected_msg = ( + "The residual of FitContribution 'test' can't be computed at " + f"{expected_invalid_count} of its 20 points, because their " + "uncertainties are zero or not a number. Please exclude those " + "points with Profile.set_calculation_range, or set valid " + "uncertainties with Profile.set_observed_profile." + ) + with pytest.raises(SrFitError, match=re.escape(expected_msg)): + contribution.residual() + + +def test_residual_with_invalid_uncertainties_outside_calculation_range(): + # Test computing the residual when the points with zero uncertainties + # are outside the calculation range, as when a user cuts off the end + # of a detector bank + # Expected: The residual is computed on the calculation range, with + # the observed uncertainties of one. + profile = Profile() + profile.set_observed_profile(XOBS, XOBS, [1.0] * 17 + [0.0] * 3) + profile.set_calculation_range(xmax=8.0) + contribution = FitContribution("test") + contribution.set_profile(profile) + contribution.set_equation("A*x") + contribution.A.set_value(2) + actual_residual = contribution.residual().tolist() + expected_residual = arange(0, 8.5, 0.5).tolist() + assert actual_residual == expected_residual + + # 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". diff --git a/tests/test_profileparser.py b/tests/test_profileparser.py index cd821387..44751e58 100644 --- a/tests/test_profileparser.py +++ b/tests/test_profileparser.py @@ -1,4 +1,5 @@ import re +import warnings from pathlib import Path import pytest @@ -75,17 +76,6 @@ [0.1, 0.3, 0.5], [0.2, 0.4, 0.6], ), - # C5: 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 None - ( - Path("four_col_nan_inf.gr"), - None, - [1.0, 1.1, 1.2], - [2.0, 2.1, 2.2], - None, - None, - ), ], ) def test_parse_file( @@ -171,6 +161,73 @@ def test_parse_file_bad(parser_datafiles, input_file, column_order, msg): parser.parse_file(parser_datafiles / input_file, column_order) +@pytest.mark.parametrize( + "input_file, expected_dx, expected_dy", + [ + # Test reading a file whose dy column has uncertainties that can't + # weight a fit + # C1: A user loads a file whose dy column has one zero. + # Expected: A warning says the whole dy column is ignored and + # suggests checking it, and dy is None. + ("three_col_zero_dy.dat", None, None), + # C2: A user loads a file whose dx and dy columns have NaN and inf + # values. + # Expected: Only the dy column gives a warning, since dx doesn't + # weight a fit, and dx and dy are None. + ("four_col_nan_inf.gr", None, None), + ], +) +def test_parse_file_with_invalid_dy( + as_list, parser_datafiles, input_file, expected_dx, expected_dy +): + # Parse the file, recording every warning to check that only one is + # raised + parser = ProfileParser() + with warnings.catch_warnings(record=True) as caught_warnings: + warnings.simplefilter("always") + parser.parse_file(parser_datafiles / input_file) + actual_warning_messages = [] + for caught_warning in caught_warnings: + if issubclass(caught_warning.category, UserWarning): + actual_warning_messages.append(str(caught_warning.message)) + actual_dx = as_list(parser._dx) + actual_dy = as_list(parser._dy) + expected_warning_messages = [ + f"The dy column of '{input_file}' has uncertainties that are " + "zero, negative, infinite or not a number, so all of them are " + "ignored and every point is weighted equally. Please check the " + "uncertainties in the file, or set them with " + "Profile.set_observed_profile." + ] + assert actual_warning_messages == expected_warning_messages + assert actual_dx == expected_dx + assert actual_dy == expected_dy + + +def test_parse_file_with_zero_dx(as_list, parser_datafiles): + # Test reading a PDF file whose dr column is all zeros, as PDFgetX2 + # writes it + # Expected: dx is None and dy is read, with no warning, since dx + # doesn't weight a fit. + parser = ProfileParser() + # Parse the file, recording every warning to check that none is raised + with warnings.catch_warnings(record=True) as caught_warnings: + warnings.simplefilter("always") + parser.parse_file(parser_datafiles / "four_col_zero_dx.gr") + actual_warning_messages = [] + for caught_warning in caught_warnings: + if issubclass(caught_warning.category, UserWarning): + actual_warning_messages.append(str(caught_warning.message)) + actual_dx = as_list(parser._dx) + actual_dy = as_list(parser._dy) + expected_warning_messages = [] + expected_dx = None + expected_dy = [0.2, 0.4, 0.6] + assert actual_warning_messages == expected_warning_messages + assert actual_dx == expected_dx + assert actual_dy == expected_dy + + # ProfileParser is an extension point. A subclass customizes a format by # overriding the _parse_metadata and _parse_data hooks; parse_file itself # is a template method that subclasses are not expected to touch. The From 0d8752e74e7ee305e4468cfad7530455f93f6401 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Fri, 30 Oct 2026 14:58:23 -0400 Subject: [PATCH 2/3] fix: give points with invalid uncertainties no weight in a fit --- news/uncertainty-checks.rst | 5 +- src/diffpy/srfit/fitbase/fitcontribution.py | 26 +----- src/diffpy/srfit/fitbase/fitrecipe.py | 9 +- src/diffpy/srfit/fitbase/fitresults.py | 23 +++++- src/diffpy/srfit/fitbase/profile.py | 57 ++++++++----- src/diffpy/srfit/fitbase/profileparser.py | 23 ++---- tests/test_contribution.py | 81 +++++++----------- tests/test_fitrecipe.py | 35 ++++++++ tests/test_fitresults.py | 61 ++++++++++++++ tests/test_profile.py | 91 ++++++++++++++++++++- tests/test_profileparser.py | 41 ++++------ 11 files changed, 310 insertions(+), 142 deletions(-) diff --git a/news/uncertainty-checks.rst b/news/uncertainty-checks.rst index f4f217a9..54329f79 100644 --- a/news/uncertainty-checks.rst +++ b/news/uncertainty-checks.rst @@ -4,7 +4,8 @@ **Changed:** -* ``ProfileParser`` warns when it ignores a ``dy`` column whose uncertainties are zero, negative, infinite or not a number, instead of weighting every point equally without saying so. +* ``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:** @@ -16,7 +17,7 @@ **Fixed:** -* ``FitContribution.residual`` raises an ``SrFitError`` that counts the points whose uncertainties are zero or not a number, instead of returning an infinite or undefined residual, unless ``Profile.set_calculation_range`` excludes them. +* 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:** diff --git a/src/diffpy/srfit/fitbase/fitcontribution.py b/src/diffpy/srfit/fitbase/fitcontribution.py index 1e40a99e..9c16daca 100644 --- a/src/diffpy/srfit/fitbase/fitcontribution.py +++ b/src/diffpy/srfit/fitbase/fitcontribution.py @@ -24,8 +24,6 @@ __all__ = ["FitContribution"] -import numpy - from diffpy.srfit.exceptions import SrFitError from diffpy.srfit.fitbase.parameter import ParameterProxy from diffpy.srfit.fitbase.parameterset import ParameterSet @@ -476,35 +474,19 @@ 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 The array of residual values. - - Raises - ------ - SrFitError - If the residual isn't finite at points whose uncertainties - ``dy`` are zero or not a number. """ # Assign the calculated profile self.profile.ycalc = self._eq() # Note that equations only recompute when their inputs are modified, so # the following will not recompute the equation. - residual = self._reseq() - if not numpy.all(numpy.isfinite(residual)): - dy = self.profile.dy - invalid_count = numpy.count_nonzero((dy == 0) | numpy.isnan(dy)) - if invalid_count > 0: - raise SrFitError( - f"The residual of FitContribution '{self.name}' can't " - f"be computed at {invalid_count} of its {len(dy)} " - "points, because their uncertainties are zero or not a " - "number. Please exclude those points with " - "Profile.set_calculation_range, or set valid " - "uncertainties with Profile.set_observed_profile." - ) - return residual + return self._reseq() def evaluate(self): """Evaluate the contribution equation and update profile.ycalc. diff --git a/src/diffpy/srfit/fitbase/fitrecipe.py b/src/diffpy/srfit/fitbase/fitrecipe.py index 78c6a516..e7e4934d 100644 --- a/src/diffpy/srfit/fitbase/fitrecipe.py +++ b/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]) diff --git a/src/diffpy/srfit/fitbase/fitresults.py b/src/diffpy/srfit/fitbase/fitresults.py index e315fa50..af89cb73 100644 --- a/src/diffpy/srfit/fitbase/fitresults.py +++ b/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) diff --git a/src/diffpy/srfit/fitbase/profile.py b/src/diffpy/srfit/fitbase/profile.py index b9583aeb..fc57a2fe 100644 --- a/src/diffpy/srfit/fitbase/profile.py +++ b/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,9 +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. - Points whose uncertainties are zero or not a number can't - weight a fit, so the residual raises an error unless - ``set_calculation_range`` excludes them. + ``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 ------ @@ -339,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) @@ -365,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 ---------- @@ -380,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 diff --git a/src/diffpy/srfit/fitbase/profileparser.py b/src/diffpy/srfit/fitbase/profileparser.py index d63024fb..849e4340 100644 --- a/src/diffpy/srfit/fitbase/profileparser.py +++ b/src/diffpy/srfit/fitbase/profileparser.py @@ -23,7 +23,6 @@ """ import warnings -from pathlib import Path import numpy as np @@ -101,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 @@ -188,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. @@ -298,21 +300,12 @@ def _parse_data(self, filename, column_format=None, **kwargs): ) column_format = self._detect_column_format(data, column_format) columns = self._map_column_labels_to_data(data, column_format) - dy = self._validate_uncertainty(columns.get("dy")) - if "dy" in columns and dy is None: - warnings.warn( - f"The dy column of '{Path(filename).name}' has uncertainties " - "that are zero, negative, infinite or not a number, so all " - "of them are ignored and every point is weighted equally. " - "Please check the uncertainties in the file, or set them " - "with Profile.set_observed_profile." - ) self._banks.append( [ columns["x"], columns["y"], self._validate_uncertainty(columns.get("dx")), - dy, + columns.get("dy"), ] ) diff --git a/tests/test_contribution.py b/tests/test_contribution.py index 0d366127..1ecbb062 100644 --- a/tests/test_contribution.py +++ b/tests/test_contribution.py @@ -14,10 +14,10 @@ ############################################################################## """Tests for the fitcontribution module.""" -import re +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 @@ -277,64 +277,43 @@ def test_residual( assert actual_builders == expected_builders -# The chiv residual divides by the uncertainties dy, so it can't be -# computed at a point whose uncertainty is zero or not a number. The -# profile is observed at the 20 points x = 0, 0.5, ..., 9.5 with y = x, and -# the equation is "A*x" with A = 2. +# 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_invalid_count", + "input_dyobs, expected_residual", [ - # C1: A user loads data whose last point has an uncertainty of - # zero. - # Expected: An SrFitError says that 1 of the 20 points has an - # uncertainty that is zero or not a number, and suggests excluding - # it or setting valid uncertainties. - ([1.0] * 19 + [0.0], "1"), - # C2: A user loads data whose last three points have uncertainties - # of zero, as at the end of a detector bank. - # Expected: An SrFitError says that 3 of the 20 points do. - ([1.0] * 17 + [0.0] * 3, "3"), - # C3: A user loads data with one uncertainty that is not a number. - # Expected: An SrFitError says that 1 of the 20 points does. - ([1.0] * 10 + [float("nan")] + [1.0] * 9, "1"), + # 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_invalid_count -): - profile = Profile() - profile.set_observed_profile(XOBS, XOBS, input_dyobs) - contribution = FitContribution("test") - contribution.set_profile(profile) - contribution.set_equation("A*x") - contribution.A.set_value(2) - expected_msg = ( - "The residual of FitContribution 'test' can't be computed at " - f"{expected_invalid_count} of its 20 points, because their " - "uncertainties are zero or not a number. Please exclude those " - "points with Profile.set_calculation_range, or set valid " - "uncertainties with Profile.set_observed_profile." - ) - with pytest.raises(SrFitError, match=re.escape(expected_msg)): - contribution.residual() - - -def test_residual_with_invalid_uncertainties_outside_calculation_range(): - # Test computing the residual when the points with zero uncertainties - # are outside the calculation range, as when a user cuts off the end - # of a detector bank - # Expected: The residual is computed on the calculation range, with - # the observed uncertainties of one. +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(XOBS, XOBS, [1.0] * 17 + [0.0] * 3) - profile.set_calculation_range(xmax=8.0) + 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) - actual_residual = contribution.residual().tolist() - expected_residual = arange(0, 8.5, 0.5).tolist() - assert actual_residual == expected_residual + # 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 diff --git a/tests/test_fitrecipe.py b/tests/test_fitrecipe.py index 9e35d972..d52632b5 100644 --- a/tests/test_fitrecipe.py +++ b/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) diff --git a/tests/test_fitresults.py b/tests/test_fitresults.py index fa00b464..99d1dad4 100644 --- a/tests/test_fitresults.py +++ b/tests/test_fitresults.py @@ -15,6 +15,7 @@ """Tests for the fitresults module.""" import pytest +from numpy import array, dot, inf, isfinite, linspace, pi, sin, sqrt from scipy.optimize import leastsq from diffpy.srfit import FitResults @@ -121,3 +122,63 @@ def test_get_results_dictionary(build_recipes_one_contribution): actual_amplitude = actual_results_dict["amplitude"] expected_amplitude = 1 assert actual_amplitude == pytest.approx(expected_amplitude, abs=1e-5) + + +# Points whose uncertainties are zero, negative, infinite or not a number +# are left out of Chi2, Reduced Chi2 and Rw, as if they weren't in the +# data, and the results say how many there are. The recipe compares +# A*sin(k*x + c), with A = 0.8, k = 1.2 and c = 0.1, to sin(x) on 21 points +# in [0, pi], without refining. The three variables are counted in Reduced +# Chi2. +@pytest.mark.parametrize( + "input_dyobs, expected_dy, expected_messages", + [ + # C1: A user fits data whose uncertainties are all valid. + # Expected: Every point counts, and there are no messages. + ([0.1] * 21, [0.1] * 21, []), + # C2: A user fits data whose last three uncertainties are zero, as + # at the end of a detector bank. + # Expected: Only the other 18 points count, and a message says how + # many points have no weight and how to set their uncertainties. + ( + [0.1] * 18 + [0.0] * 3, + [0.1] * 18 + [inf] * 3, + [ + "3 of the 21 points of contribution 'c1' 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." + ], + ), + # C3: A user fits data whose uncertainties are all zero. + # Expected: Every point counts with an uncertainty of one, and the + # message for missing uncertainties is given. + ( + [0.0] * 21, + [1.0] * 21, + ["Some quantities invalid due to missing profile uncertainty"], + ), + ], +) +def test_results_with_invalid_uncertainties( + build_recipe_with_uncertainty, input_dyobs, expected_dy, expected_messages +): + # Compute the results at the starting values + recipe, profile = build_recipe_with_uncertainty(input_dyobs) + results = FitResults(recipe) + results.get_results_string() + actual_messages = results.messages + actual_metrics = [results.chi2, results.rchi2, results.rw] + # Compute the expected metrics from the points with weight only + x = linspace(0, pi, 21) + y = sin(x) + dy = array(expected_dy) + is_weighted = isfinite(dy) + chiv = (0.8 * sin(1.2 * x + 0.1) - y)[is_weighted] / dy[is_weighted] + weighted_y = y[is_weighted] / dy[is_weighted] + expected_chi2 = dot(chiv, chiv) + expected_rchi2 = expected_chi2 / (is_weighted.sum() - 3) + expected_rw = sqrt(expected_chi2 / dot(weighted_y, weighted_y)) + expected_metrics = [expected_chi2, expected_rchi2, expected_rw] + assert actual_messages == expected_messages + assert actual_metrics == pytest.approx(expected_metrics) diff --git a/tests/test_profile.py b/tests/test_profile.py index f67d594c..b32332da 100644 --- a/tests/test_profile.py +++ b/tests/test_profile.py @@ -18,7 +18,7 @@ import re import pytest -from numpy import allclose, arange, array, array_equal +from numpy import allclose, arange, array, array_equal, inf, nan from diffpy.srfit import FitContribution, Profile, ProfileParser from diffpy.srfit.exceptions import SrFitError @@ -85,6 +85,95 @@ def test_set_observed_profile( assert actual_dy == expected_dy +# A point whose observed uncertainty isn't finite and positive gets no +# weight in a fit, so its calculated dy is infinite, unless no point has a +# valid uncertainty, when the profile is weighted as if it had none. The +# observed uncertainties dyobs are kept as given. The profile is observed +# at x = 1, 2, 3, 4 with y = x. +@pytest.mark.parametrize( + "input_dyobs, expected_dy", + [ + # C1: A user sets data with one uncertainty of zero. + # Expected: dy is infinite at that point, and the other points keep + # their uncertainties. + ([0.1, 0.0, 0.3, 0.4], [0.1, inf, 0.3, 0.4]), + # C2: A user sets data with a negative, an infinite and a NaN + # uncertainty. + # Expected: dy is infinite at those three points. + ([0.1, -0.2, inf, nan], [0.1, inf, inf, inf]), + # C3: A user sets data whose uncertainties are all zero, as in a + # file whose uncertainty column is all zeros. + # Expected: dy is one everywhere, as for data with no uncertainties. + ([0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0]), + # C4: A user sets data whose uncertainties are all invalid in + # different ways. + # Expected: dy is one everywhere. + ([0.0, -0.2, inf, nan], [1.0, 1.0, 1.0, 1.0]), + ], +) +def test_set_observed_profile_with_invalid_uncertainties( + input_dyobs, expected_dy +): + input_xobs = [1.0, 2.0, 3.0, 4.0] + profile = Profile() + profile.set_observed_profile(input_xobs, input_xobs, input_dyobs) + actual_dyobs = profile.dyobs + actual_dy = profile.dy + expected_dyobs = input_dyobs + assert array_equal(actual_dyobs, expected_dyobs, equal_nan=True) + assert array_equal(actual_dy, expected_dy) + + +# Whether a profile is fully, partly or not at all weighted depends only on +# the uncertainties inside its calculation range. The profile is observed +# at x = 1, 2, ..., 6 with uncertainties of 0.1, except 0 at x = 5 and 6, as +# at the end of a detector bank. +@pytest.mark.parametrize( + "input_range, expected_dy", + [ + # C1: A user fits the whole range. + # Expected: dy is infinite at x = 5 and 6, so only those points have + # no weight. + ({}, [0.1, 0.1, 0.1, 0.1, inf, inf]), + # C2: A user cuts off the end of the bank. + # Expected: dy keeps the uncertainties of 0.1. + ({"xmax": 4}, [0.1, 0.1, 0.1, 0.1]), + # C3: A user fits only the end of the bank. + # Expected: dy is one everywhere, as for data with no uncertainties. + ({"xmin": 5}, [1.0, 1.0]), + ], +) +def test_set_calculation_range_with_invalid_uncertainties( + input_range, expected_dy +): + input_xobs = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + profile = Profile() + profile.set_observed_profile( + input_xobs, input_xobs, [0.1, 0.1, 0.1, 0.1, 0.0, 0.0] + ) + profile.set_calculation_range(**input_range) + actual_dy = profile.dy + assert array_equal(actual_dy, expected_dy) + + +def test_set_calculation_points_with_invalid_uncertainties(): + # Test resampling a profile whose uncertainty at x = 2 is zero onto a + # grid twice as fine + # Expected: The calculation points next to x = 2 have no weight, so dy + # is infinite there, and dy at the other points is interpolated from + # the valid uncertainties only. + profile = Profile() + profile.set_observed_profile( + [0.0, 1.0, 2.0, 3.0, 4.0], + [0.0, 1.0, 2.0, 3.0, 4.0], + [0.1, 0.2, 0.0, 0.4, 0.5], + ) + profile.set_calculation_points(arange(0, 4.1, 0.5)) + actual_dy = profile.dy + expected_dy = [0.1, 0.15, 0.2, inf, inf, inf, 0.4, 0.45, 0.5] + assert allclose(actual_dy, expected_dy) + + # set_calculation_range needs an observed profile to work on. @pytest.mark.parametrize( "input_args", diff --git a/tests/test_profileparser.py b/tests/test_profileparser.py index 44751e58..cbb7640e 100644 --- a/tests/test_profileparser.py +++ b/tests/test_profileparser.py @@ -3,6 +3,7 @@ from pathlib import Path import pytest +from numpy import array_equal, inf, nan from diffpy.srfit import ProfileParser from diffpy.srfit.exceptions import ParseError @@ -162,27 +163,24 @@ def test_parse_file_bad(parser_datafiles, input_file, column_order, msg): @pytest.mark.parametrize( - "input_file, expected_dx, expected_dy", + "input_file, expected_dy", [ - # Test reading a file whose dy column has uncertainties that can't - # weight a fit + # Test reading a file whose dy column has uncertainties that are + # zero, negative, infinite or not a number, which a Profile gives + # no weight in a fit # C1: A user loads a file whose dy column has one zero. - # Expected: A warning says the whole dy column is ignored and - # suggests checking it, and dy is None. - ("three_col_zero_dy.dat", None, None), + # Expected: dy is read as it is in the file, with no warning. + ("three_col_zero_dy.dat", [0.2, 0.0, 0.6]), # C2: A user loads a file whose dx and dy columns have NaN and inf # values. - # Expected: Only the dy column gives a warning, since dx doesn't - # weight a fit, and dx and dy are None. - ("four_col_nan_inf.gr", None, None), + # Expected: dy is read as it is in the file, with no warning, and + # dx is None, since dx doesn't weight a fit. + ("four_col_nan_inf.gr", [inf, 1.0, nan]), ], ) -def test_parse_file_with_invalid_dy( - as_list, parser_datafiles, input_file, expected_dx, expected_dy -): - # Parse the file, recording every warning to check that only one is - # raised +def test_parse_file_with_invalid_dy(parser_datafiles, input_file, expected_dy): parser = ProfileParser() + # Parse the file, recording every warning to check that none is raised with warnings.catch_warnings(record=True) as caught_warnings: warnings.simplefilter("always") parser.parse_file(parser_datafiles / input_file) @@ -190,18 +188,13 @@ def test_parse_file_with_invalid_dy( for caught_warning in caught_warnings: if issubclass(caught_warning.category, UserWarning): actual_warning_messages.append(str(caught_warning.message)) - actual_dx = as_list(parser._dx) - actual_dy = as_list(parser._dy) - expected_warning_messages = [ - f"The dy column of '{input_file}' has uncertainties that are " - "zero, negative, infinite or not a number, so all of them are " - "ignored and every point is weighted equally. Please check the " - "uncertainties in the file, or set them with " - "Profile.set_observed_profile." - ] + actual_dx = parser._dx + actual_dy = parser._dy + expected_warning_messages = [] + expected_dx = None assert actual_warning_messages == expected_warning_messages assert actual_dx == expected_dx - assert actual_dy == expected_dy + assert array_equal(actual_dy, expected_dy, equal_nan=True) def test_parse_file_with_zero_dx(as_list, parser_datafiles): From 38210a3dc1ad7cb14db258cbf28aaaaf86e1c012 Mon Sep 17 00:00:00 2001 From: Caden Myers Date: Fri, 30 Oct 2026 15:01:54 -0400 Subject: [PATCH 3/3] trigger ci