Pointwise update

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",

src/easyreflectometry/model/resolution_functions.py

@@ -82,27 +82,45 @@
         }
 
 
-# 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 @@
             '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

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)