Merge pull request #690 from scipp/peak-fit-robust-polynomial

Peak fitting: Do not drop params with value=0

Author: jl-wynen

src/scippneutron/peaks/model.py

@@ -311,9 +311,14 @@ def _call(self, x: sc.Variable, params: dict[str, sc.Variable]) -> sc.Variable:
 
     def _guess(self, x: sc.Variable, y: sc.Variable) -> dict[str, sc.Variable]:
         poly = np.polynomial.Polynomial.fit(x.values, y.values, deg=self.degree)
+        # `fit` returns a polynomial with domain=[x.min, x.max], convert to [-1, 1]
+        # because all models are defined for peaks centered around 0.
+        coef = poly.convert().coef
+        # Pad with 0's because `convert` drops those to simplify the polynomial
+        # but the model requires an exact number of params.
+        coef = np.pad(coef, (0, len(self._param_names) - len(coef)))
         return {
-            f'a{i}': sc.scalar(c, unit=y.unit / x.unit**i)
-            for i, c in enumerate(poly.convert().coef)
+            f'a{i}': sc.scalar(c, unit=y.unit / x.unit**i) for i, c in enumerate(coef)
         }
 
 

tests/peaks/model_test.py

@@ -209,6 +209,26 @@ def test_polynomial_degree_2_guess_params():
     sc.testing.assert_allclose(params['a2'], a2, atol=sc.scalar(0.1, unit='K/m^2'))
 
 
+def test_polynomial_degree_2_guess_params_degenerate():
+    x = sc.linspace('xx', 0.8, 4.3, 10, unit='s')
+    y = sc.zeros(sizes=x.sizes, unit='cm')
+    data = sc.DataArray(y, coords={'xx': x})
+
+    m = model.PolynomialModel(degree=2, prefix='')
+    params = m.guess(data)
+    # Test that no params have been optimized away:
+    assert params.keys() == {'a0', 'a1', 'a2'}
+    sc.testing.assert_allclose(
+        params['a0'], sc.scalar(0.0, unit='cm'), atol=sc.scalar(0.1, unit='cm')
+    )
+    sc.testing.assert_allclose(
+        params['a1'], sc.scalar(0.0, unit='cm/s'), atol=sc.scalar(0.1, unit='cm/s')
+    )
+    sc.testing.assert_allclose(
+        params['a2'], sc.scalar(0.0, unit='cm/s^2'), atol=sc.scalar(0.1, unit='cm/s^2')
+    )
+
+
 def test_gaussian_guess_params_linspace():
     amplitude = sc.scalar(2.8, unit='kg')
     loc = sc.scalar(0.4, unit='m')