diff --git a/docs/docs/tutorials/simulation/resolution_functions.ipynb b/docs/docs/tutorials/simulation/resolution_functions.ipynb index b5bd017e..d46a76cb 100644 --- a/docs/docs/tutorials/simulation/resolution_functions.ipynb +++ b/docs/docs/tutorials/simulation/resolution_functions.ipynb @@ -340,7 +340,10 @@ "## Afterthoughts\n", "As a last task we will compare the reflectivity determined using a percentage resolution function and a point-wise function.\n", "We should recall that the \"experimental\" data was generated using `Refnx`.\n", - "By comparing the reflectivities determined using a resolution function with a FWHM of 1.0% and the point-wise FHWN constructed from data in a `.ort` file it is apparent that this reference data also was constructed using a resolution function of 1.0%." + "\n", + "The `Pointwise` resolution function derives a per-point resolution width directly from the data: it takes the `[Qz, R, sQz]` triple, where `sQz` is the variance of `Qz` (`Qz_0.variances`), and uses `sqrt(sQz)` as the resolution width at each measured point, linearly interpolating onto the requested `q` (exactly as `LinearSpline` does for explicitly provided widths).\n", + "\n", + "By comparing the reflectivities determined using a resolution function with a FWHM of 1.0% and the point-wise width constructed from the data in a `.ort` file, it is apparent that this reference data also was constructed using a resolution function of 1.0%." ] }, { @@ -404,16 +407,20 @@ " model.unique_name,\n", ")\n", "plt.plot(model_coords, model_data, 'k-', label='Variable', linewidth=5)\n", + "\n", + "# The Pointwise resolution is built from the data triple [Qz, R, sQz],\n", + "# where sQz is the variance of Qz. The width sqrt(sQz) is interpolated onto q.\n", "data_points = []\n", "data_points.append(reference_coords) # Qz\n", "data_points.append(reference_data) # R\n", - "data_points.append(reference_variances) # sQz\n", + "data_points.append(reference_variances) # sQz (variance of Qz)\n", "model.resolution_function = Pointwise(q_data_points=data_points)\n", "model_data = model.interface().reflectity_profile(\n", " model_coords,\n", " model.unique_name,\n", ")\n", "plt.plot(model_coords, model_data, 'r-', label='Pointwise')\n", + "plt.plot(reference_coords, reference_data, 'bx', label='Reference', markersize=4)\n", "\n", "ax = plt.gca()\n", "ax.set_xlim([-0.01, 0.45])\n", diff --git a/src/easyreflectometry/model/resolution_functions.py b/src/easyreflectometry/model/resolution_functions.py index ee0933e2..fe5d2254 100644 --- a/src/easyreflectometry/model/resolution_functions.py +++ b/src/easyreflectometry/model/resolution_functions.py @@ -82,27 +82,45 @@ def as_dict( } -# add pointwise smearing funtion class Pointwise(ResolutionFunction): - def __init__(self, q_data_points: list[np.ndarray]): - """Init function.""" + """Pointwise resolution defined by a per-point resolution provided with the data. + + The resolution is supplied as the variance of the Qz values (``sQz``) at the + measured Qz data points, which is the form produced by data reduction (e.g. + ``Qz_0.variances``). The resolution width at each point is ``sqrt(sQz)``. + For a requested ``q`` the width is obtained by linearly interpolating onto + ``q``, exactly as :class:`LinearSpline` does for explicitly provided widths. + + This is a convenience wrapper around :class:`LinearSpline` that derives the + widths from the ``[Qz, R, sQz]`` triple loaded from a data file; the returned + widths are consumed by the calculators (refnx ``x_err`` / refl1d ``dq``), + which perform the actual convolution against the model. + """ + + def __init__(self, q_data_points: List[np.ndarray]): + """Init function. + + Parameters + ---------- + q_data_points : List[np.ndarray] + ``[Qz, R, sQz]`` where ``Qz`` are the measured Qz values, ``R`` the + measured reflectivity (kept only for serialization round-trips) and + ``sQz`` the variance of ``Qz`` at each point. + """ self.q_data_points = q_data_points - self.q = None - def smearing(self, q: Union[np.ndarray, float] = None) -> np.ndarray: - """Smearing function.""" - Qz = self.q_data_points[0] - R = self.q_data_points[1] - sQz = self.q_data_points[2] - if q is None: - q = self.q_data_points[0] - self.q = q - sQzs = np.sqrt(sQz) - if isinstance(Qz, float): - Qz = np.array(Qz) - - smeared = self.apply_smooth_smearing(Qz, R, sQzs) - return smeared + def smearing(self, q: Optional[Union[np.ndarray, float]] = None) -> np.ndarray: + """Return the resolution width interpolated onto ``q``. + + The width at each data point is ``sqrt(sQz)``; values are linearly + interpolated onto the requested ``q``. When ``q`` is ``None`` the widths + are returned at the stored data points. + """ + Qz = np.asarray(self.q_data_points[0], dtype=float) + sQz = np.asarray(self.q_data_points[2], dtype=float) + q_eval = Qz if q is None else np.asarray(q, dtype=float) + widths = np.sqrt(sQz) + return np.asarray(np.interp(q_eval, Qz, widths)) def as_dict( self, skip: Optional[List[str]] = None @@ -114,29 +132,3 @@ def as_dict( 'R_data_points': list(self.q_data_points[1]), 'sQz_data_points': list(self.q_data_points[2]), } - - def gaussian_smearing(self, qt, Qz, R, sQz): - """Gaussian smearing.""" - weights = np.exp(-0.5 * ((qt - Qz) / sQz) ** 2) - if np.sum(weights) == 0 or not np.isfinite(np.sum(weights)): - return np.sum(R) - weights /= sQz * np.sqrt(2 * np.pi) - return np.sum(R * weights) / np.sum(weights) - - def apply_smooth_smearing(self, Qz, R, sQzs): - """Apply smooth resolution smearing using convolution with Gaussian kernel.""" - if self.q is None: - R_smeared = np.zeros_like(Qz) - else: - R_smeared = np.zeros_like(self.q) - - if not isinstance(Qz, np.ndarray): - Qz = np.array(Qz) - if not isinstance(R, np.ndarray): - R = np.array(R) - R_smeared = np.zeros_like(self.q) - - for i, qt in enumerate(self.q): - R_smeared[i] = self.gaussian_smearing(qt, Qz, R, sQzs) - - return R_smeared diff --git a/tests/model/test_resolution_functions.py b/tests/model/test_resolution_functions.py index e28a99f2..b8a4ea18 100644 --- a/tests/model/test_resolution_functions.py +++ b/tests/model/test_resolution_functions.py @@ -95,12 +95,23 @@ def test_constructor(self): # When resolution_function = Pointwise(q_data_points=self.data_points) - # Then Expect + # Then Expect: smearing returns the resolution width sqrt(sQz) at the data + # points, since sQz holds the variance of Qz (consistent with LinearSpline). + expected_widths = np.sqrt(self.data_points[2]) assert np.allclose( np.array(resolution_function.smearing()), - np.array([2.51664683, 2.84038734, 3.2460762, 3.6796519, 4.07869271]), + expected_widths, ) + def test_smearing_interpolates_onto_q(self): + # When + resolution_function = Pointwise(q_data_points=self.data_points) + + # Then Expect: requesting points between data points linearly interpolates the width. + widths = np.sqrt(self.data_points[2]) + expected = np.interp([0.15, 0.25], self.data_points[0], widths) + assert np.allclose(resolution_function.smearing([0.15, 0.25]), expected) + def test_as_dict(self): # When resolution_function = Pointwise(q_data_points=self.data_points)