Skip to content

Commit 02ca2b4

Browse files
authored
change: Handle uncertainties as None rather than 1s (#193)
* fix: return None instead of 0 for unavailable dx/dy uncertainties * test: update and add tests for None uncertainty behavior * chore: add news entry for uncertainty None fix * test: add tests to ensure new handling of uncertainties doesn't break refinements * change: update profile.set_observed_profile() to not convert dy to ones if dy is None * chore: update news file * Update tests to make it easier to not use a testdata file for interpretability * add as_list as a fixture * add comment about how the fixture is used in test * add three_column case * use as_list as a fixture
1 parent 99955e4 commit 02ca2b4

9 files changed

Lines changed: 260 additions & 50 deletions

File tree

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
**Added:**
2+
3+
* <news item>
4+
5+
**Changed:**
6+
7+
* Change ``ProfileParser.parse_file`` to set unavailable ``dx``/``dy`` uncertainties to ``None`` instead of ``0``, matching the format-specific parsers.
8+
* 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.
9+
* Change ``Profile._validate`` to accept a ``None`` ``dyobs``, since observed uncertainties are optional. ``x``, ``y``, ``dy``, ``xobs`` and ``yobs`` are still required.
10+
11+
**Deprecated:**
12+
13+
* <news item>
14+
15+
**Removed:**
16+
17+
* <news item>
18+
19+
**Fixed:**
20+
21+
* 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.
22+
* 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.
23+
24+
**Security:**
25+
26+
* <news item>

‎src/diffpy/srfit/fitbase/profile.py‎

Lines changed: 14 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -186,8 +186,9 @@ def set_observed_profile(self, xobs, yobs, dyobs=None):
186186
Numpy array of the observed signal.
187187
dyobs
188188
Numpy array of the uncertainty in the observed signal. If
189-
`dyobs` is None (default), it will be set to 1 at each
190-
observed `xobs`.
189+
`dyobs` is None (default), `dyobs` stays None to indicate
190+
no uncertainty was observed, and the calculated `dy` will
191+
be set to 1 at each calculation point instead.
191192
192193
193194
Raises
@@ -206,7 +207,7 @@ def set_observed_profile(self, xobs, yobs, dyobs=None):
206207
self._yobs = numpy.asarray(yobs, dtype=float)
207208

208209
if dyobs is None:
209-
self._dyobs = numpy.ones_like(xobs)
210+
self._dyobs = None
210211
else:
211212
self._dyobs = numpy.asarray(dyobs, dtype=float)
212213

@@ -328,7 +329,11 @@ def _isobs(a):
328329
indices = (lo - epslo <= self.xobs) & (self.xobs <= hi + epshi)
329330
self.x = self.xobs[indices]
330331
self.y = self.yobs[indices]
331-
self.dy = self.dyobs[indices]
332+
self.dy = (
333+
self.dyobs[indices]
334+
if self.dyobs is not None
335+
else numpy.ones_like(self.x)
336+
)
332337
else:
333338
x1 = numpy.arange(lo, hi + epshi, step)
334339
self.set_calculation_points(x1)
@@ -373,6 +378,8 @@ def set_calculation_points(self, x):
373378
# FIXME - This does not follow error propagation rules and it
374379
# introduces (more) correlation between the data points.
375380
self.dy = _rebin_array(self.dyobs, self.xobs, self.x)
381+
elif self.yobs is not None:
382+
self.dy = numpy.ones_like(self.x)
376383

377384
return
378385

@@ -477,8 +484,9 @@ def _flush(self, other):
477484
def _validate(self):
478485
"""Validate my state.
479486
480-
This validates that x, y, dy, xobx, yobs and dyobs are not None.
481-
This validates that x, y, and dy are the same length.
487+
This validates that x, y, dy, xobs and yobs are not None. dyobs
488+
may be None, since observed uncertainties are optional. This
489+
validates that x, y, and dy are the same length.
482490
483491
Raises SrFitError if validation fails.
484492
"""
@@ -490,7 +498,6 @@ def _validate(self):
490498
self.dy,
491499
self.xobs,
492500
self.yobs,
493-
self.dyobs,
494501
]
495502
)
496503
if datanotset:

‎src/diffpy/srfit/fitbase/profileparser.py‎

Lines changed: 12 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -94,11 +94,11 @@ class ProfileParser(object):
9494
from the file.
9595
dx : np.ndarray
9696
The uncertainties associated with x
97-
read from the file. This is 0 if the
97+
read from the file. This is None if the
9898
uncertainty cannot be read.
9999
dy : np.ndarray
100100
The uncertainties associated with y
101-
read from the file. This is 0 if the
101+
read from the file. This is None if the
102102
uncertainty cannot be read.
103103
_x : np.ndarray
104104
Independent variable from the chosen bank
@@ -197,13 +197,13 @@ def parse_file(self, filename, column_format=None):
197197
"""Parse a data file to extract data and metadata, with
198198
automatic handling of uncertainties.
199199
200-
- For files with 2 columns: assumes (x, y) and sets dx, dy to 0.
201-
- For files with 3 columns: assumes (x, y, dy) and sets dx to 0.
200+
- For files with 2 columns: assumes (x, y) and sets dx, dy to None.
201+
- For files with 3 columns: assumes (x, y, dy) and sets dx to None.
202202
- For files with 4 columns: assumes (x, y, dx, dy).
203203
- For other cases: `column_format` must be explicitly specified.
204204
205205
Uncertainty columns (dx, dy) are only considered valid if all values
206-
are positive and not NaN/Inf. Otherwise they are set to 0.
206+
are positive and not NaN/Inf. Otherwise they are set to None.
207207
208208
This wipes out the currently loaded data and selected bank number.
209209
@@ -242,10 +242,8 @@ def parse_file(self, filename, column_format=None):
242242
# Extract required arrays
243243
x = columns["x"]
244244
y = columns["y"]
245-
x_length = len(x)
246-
y_length = len(y)
247-
dx = self._validate_uncertainty(columns.get("dx"), x_length)
248-
dy = self._validate_uncertainty(columns.get("dy"), y_length)
245+
dx = self._validate_uncertainty(columns.get("dx"))
246+
dy = self._validate_uncertainty(columns.get("dy"))
249247
# Store as single bank
250248
self._banks = [(x, y, dx, dy)]
251249
self._meta["nbanks"] = 1
@@ -306,10 +304,10 @@ def _map_column_labels_to_data(self, data, column_format):
306304
return columns
307305

308306
@staticmethod
309-
def _validate_uncertainty(data, length):
310-
"""Return the uncertainty data if valid, otherwise 0."""
307+
def _validate_uncertainty(data):
308+
"""Return the uncertainty data if valid, otherwise None."""
311309
if data is None or not np.all(np.isfinite(data)) or np.any(data <= 0):
312-
return np.zeros(length)
310+
return None
313311
return data
314312

315313
def get_num_banks(self):
@@ -395,8 +393,8 @@ def get_data(self, index=None):
395393
396394
Returns
397395
----------
398-
This returns (x, y, dx, dy) tuple for the bank. dx is 0 if it cannot
399-
be determined from the data format.
396+
This returns (x, y, dx, dy) tuple for the bank. dx is None if it
397+
cannot be determined from the data format.
400398
"""
401399
self.select_bank(index)
402400

‎src/diffpy/srfit/pdf/pdfparser.py‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -63,11 +63,11 @@ class PDFParser(ProfileParser):
6363
from the file.
6464
dx
6565
A numpy array containing the uncertainty in x
66-
read from the file. This is 0 if the
66+
read from the file. This is None if the
6767
uncertainty cannot be read.
6868
dy
6969
A numpy array containing the uncertainty read
70-
from the file. This is 0 if the uncertainty
70+
from the file. This is None if the uncertainty
7171
cannot be read.
7272
_x
7373
Independent variable from the chosen bank
@@ -128,7 +128,8 @@ def parseString(self, patstring):
128128
"""Parse a string and set the _x, _y, _dx, _dy and _meta
129129
variables.
130130
131-
When _dx or _dy cannot be obtained in the data format it is set to 0.
131+
When _dx or _dy cannot be obtained in the data format it is set to
132+
None.
132133
133134
This wipes out the currently loaded data and selected bank number.
134135

‎tests/conftest.py‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -219,6 +219,43 @@ def build_recipe_two_contributions():
219219
return recipe
220220

221221

222+
@pytest.fixture(scope="session")
223+
def as_list():
224+
def _as_list(values):
225+
"""Unavailable uncertainties are None rather than an array."""
226+
return None if values is None else values.tolist()
227+
228+
return _as_list
229+
230+
231+
@pytest.fixture(scope="function")
232+
def build_recipe_with_uncertainty():
233+
"""Helper to build a sine recipe whose profile carries the observed
234+
uncertainties returned by ``make_dyobs``.
235+
236+
The observed profile is a noiseless sine, so a refinement recovers
237+
``A=1``, ``k=1`` and ``c=0`` regardless of how it is weighted.
238+
"""
239+
240+
def _build_recipe(make_dyobs):
241+
xobs = linspace(0, pi, 21)
242+
yobs = sin(xobs)
243+
profile = Profile()
244+
profile.set_observed_profile(xobs, yobs, make_dyobs(xobs))
245+
contribution = FitContribution("c1")
246+
contribution.set_profile(profile)
247+
contribution.set_equation("A*sin(k*x + c)")
248+
recipe = FitRecipe()
249+
recipe.fithooks[0].verbose = 0
250+
recipe.add_contribution(contribution)
251+
recipe.add_variable(contribution.A, 0.8)
252+
recipe.add_variable(contribution.k, 1.2)
253+
recipe.add_variable(contribution.c, 0.1)
254+
return recipe, profile
255+
256+
return _build_recipe
257+
258+
222259
@pytest.fixture
223260
def temp_data_files(tmp_path):
224261
"""

‎tests/test_fitrecipe.py‎

Lines changed: 83 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@
2020
import matplotlib.pyplot as plt
2121
import numpy as np
2222
import pytest
23-
from numpy import array_equal, dot, linspace, pi, sin
23+
from numpy import array_equal, dot, linspace, ones_like, pi, sin
2424
from scipy.optimize import leastsq
2525

2626
from diffpy.srfit.fitbase import FitResults, ProfileParser
@@ -516,7 +516,7 @@ def test_initialize_recipe_from_recipe(build_recipes_one_contribution):
516516
assert sorted(list(expected_values)) == sorted(list(actual_values))
517517

518518

519-
def test_initialize_recipe_from_recipe_bad(build_recipe_two_contributions):
519+
def test_initialize_recipe_from_recipe_bad():
520520
# Case: User tries to initialize a FitRecipe from a non recipe object
521521
# expected: raised ValueError with message
522522
recipe_bad = 12345 # not a FitRecipe object
@@ -1060,5 +1060,86 @@ def test_plot_recipe_reset_all_defaults(build_recipes_one_contribution):
10601060
assert actual_legend == expected_legend
10611061

10621062

1063+
# Observed uncertainties are optional. A profile loaded without an
1064+
# uncertainty column keeps dyobs as None and is refined unweighted, with
1065+
# dy falling back to one at every calculation point. The cases below
1066+
# refine the same noiseless sine profile so that the known solution
1067+
# (A=1, k=1, c=0) is recovered no matter how the fit is weighted.
1068+
1069+
1070+
# make_input_dyobs returns a function that evaluates dyobs on the array
1071+
# of xobs for insertion into the build_recipe_with_uncertainty fixture.
1072+
@pytest.mark.parametrize(
1073+
"make_input_dyobs, expected_dyobs_is_set",
1074+
[
1075+
# C1: No uncertainties are observed, as for a file with no
1076+
# uncertainty column.
1077+
# Expected: The refinement converges and dyobs stays None.
1078+
(lambda xobs: None, False),
1079+
# C2: Uniform uncertainties of one are observed explicitly.
1080+
# Expected: The refinement converges to the same values as C1,
1081+
# since an unweighted fit is weighted by one everywhere.
1082+
(lambda xobs: ones_like(xobs), True),
1083+
# C3: Uncertainties vary across the profile, so the refinement
1084+
# is weighted point by point.
1085+
# Expected: The refinement still converges to the known values.
1086+
(lambda xobs: linspace(0.1, 1.0, len(xobs)), True),
1087+
],
1088+
)
1089+
def test_refine_with_and_without_uncertainty(
1090+
build_recipe_with_uncertainty, make_input_dyobs, expected_dyobs_is_set
1091+
):
1092+
recipe, profile = build_recipe_with_uncertainty(make_input_dyobs)
1093+
optimize_recipe(recipe)
1094+
actual_dyobs_is_set = profile.dyobs is not None
1095+
actual_names = recipe.get_names()
1096+
actual_values = recipe.get_values()
1097+
expected_names = ["A", "k", "c"]
1098+
expected_values = [1.0, 1.0, 0.0]
1099+
assert actual_dyobs_is_set == expected_dyobs_is_set
1100+
assert actual_names == expected_names
1101+
assert actual_values == pytest.approx(expected_values, abs=1e-5)
1102+
1103+
1104+
# The residual that is optimized is (ycalc - y)/dy, so dy is what
1105+
# actually weights a refinement. These cases check that a profile with
1106+
# no observed uncertainties is weighted identically to one whose
1107+
# uncertainties are all one, and that observed uncertainties are carried
1108+
# through to the residual unchanged.
1109+
1110+
1111+
# make_input_dyobs returns a function that evaluates dyobs on the array
1112+
# of xobs for insertion into the build_recipe_with_uncertainty fixture.
1113+
@pytest.mark.parametrize(
1114+
"make_input_dyobs, make_expected_dy",
1115+
[
1116+
# C1: No uncertainties are observed.
1117+
# Expected: dy falls back to one everywhere, so the residual is
1118+
# unweighted.
1119+
(lambda xobs: None, lambda xobs: ones_like(xobs)),
1120+
# C2: Uniform uncertainties of one are observed explicitly.
1121+
# Expected: dy is one everywhere, matching C1.
1122+
(lambda xobs: ones_like(xobs), lambda xobs: ones_like(xobs)),
1123+
# C3: Uncertainties vary across the profile.
1124+
# Expected: dy keeps the observed values, so the residual is
1125+
# weighted point by point.
1126+
(
1127+
lambda xobs: linspace(0.1, 1.0, len(xobs)),
1128+
lambda xobs: linspace(0.1, 1.0, len(xobs)),
1129+
),
1130+
],
1131+
)
1132+
def test_residual_is_weighted_by_uncertainty(
1133+
build_recipe_with_uncertainty, make_input_dyobs, make_expected_dy
1134+
):
1135+
recipe, profile = build_recipe_with_uncertainty(make_input_dyobs)
1136+
actual_residual = recipe.residual(recipe.values)
1137+
actual_dy = profile.dy
1138+
expected_dy = make_expected_dy(profile.xobs)
1139+
expected_residual = (profile.ycalc - profile.y) / expected_dy
1140+
assert actual_dy == pytest.approx(expected_dy)
1141+
assert actual_residual == pytest.approx(expected_residual)
1142+
1143+
10631144
if __name__ == "__main__":
10641145
unittest.main()

‎tests/test_pdf.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -139,7 +139,8 @@ def testParser2(datafile):
139139
res = numpy.dot(diff, diff)
140140
assert 0 == pytest.approx(res)
141141

142-
assert dx.tolist() == [0] * len(dx)
142+
# si-q27r60-xray.gr has a negative dx column, so it is invalid.
143+
assert dx is None
143144
return
144145

145146

0 commit comments

Comments
 (0)