improve tests for the peak position error; fine-tune the major peak-detecting function

Author: comcon1

src/fairmd/lipids/analib/formfactor.py

@@ -3,12 +3,11 @@
 import numpy as np
 import scipy.interpolate
 import scipy.signal
-from scipy.optimize import curve_fit
 
 
 def get_mins_from_ffdata(ffdata: np.ndarray) -> list[float]:
     """Find the positions of minimums in form factor data."""
-    sg_window_q = 0.03  # Savitsky-Golay window (in Q)
+    sg_window_q = 0.05  # Savitsky-Golay window (in Q)
     delta_q = ffdata[1, 0] - ffdata[0, 0]  # Q step in FF data
     sg_window_n = int(np.ceil(sg_window_q / delta_q))  # S-G window (in num frames)
     try:
@@ -19,7 +18,8 @@ def get_mins_from_ffdata(ffdata: np.ndarray) -> list[float]:
 
     min_q_distance = 0.01  # Min distance btw peaks (in Q)
     mqd_n = int(np.ceil(min_q_distance / delta_q))  # same in num frames
-    peak_ind = scipy.signal.find_peaks(-filtered, distance=mqd_n)
+    peak_prominence = (filtered.max() - filtered.min()) * 0.1
+    peak_ind = scipy.signal.find_peaks(-filtered, distance=mqd_n, prominence=peak_prominence)
     min_peak_q = 0.1
 
     return [ffdata[i, 0] for i in peak_ind[0] if ffdata[i, 0] > min_peak_q]
@@ -77,5 +77,9 @@ def calc_minpos_with_error(ffdata: np.ndarray) -> (float, float):
     )
     a, b, _c = popt
     min_x = -b / 2 / a
-    delta_minx = (-1/2/a)**2*pcov[0, 0] + (b/2/a**2)**2*pcov[1, 1] + 2*(-1/2/a)*(b/2/a**2)*pcov[0, 1]
+    delta_minx = (
+        (-1 / 2 / a) ** 2 * pcov[0, 0]
+        + (b / 2 / a**2) ** 2 * pcov[1, 1]
+        + 2 * (-1 / 2 / a) * (b / 2 / a**2) * pcov[0, 1]
+    )
     return min_x, np.sqrt(delta_minx)

tests/test_formfactor.py

@@ -57,22 +57,32 @@ def test_estimate_error_of_min():
         abs=1e-3,
         msg="Error of minimum with no error is not estimated as err=0.1 [const]",
     )
-    # Case 1: constant error
-    # sin(8x) = 0
-    # 8x = pi*n => x = pi*n/8
-    #
-    # Error must be estimated as following: we down the curve by the value of error and find intersections
-    # with x axis. The error is the mean distance between the minimum and these intersections. In this
-    # case, the intersections are defined by the equation:
-    #
-    # |sin(8x)|-0.1 = 0 =>> |sin(8x)| = 0.1 =>> sin(8x) = 0.1 or sin(8x) = -0.1, i.e.,
-    #
-    # x = +-arcsin(0.1)/8 + pi*n/8; so error should be arcsin(0.1)/8
-    #
+    # Case 1: 2 constant errors
     pos_comp, err_comp = calc_minpos_with_error(synt_data)
+    synt_data[:, 2] = 0.2
+    pos_berr, err_berr = calc_minpos_with_error(synt_data)
+    check.less(
+        err_comp,
+        err_berr,
+        msg="Error of minimum with constant error 0.2 must be > than with err 0.1",
+    )
     check.almost_equal(
+        pos_comp,
+        pos_berr,
+        abs=1e-5,
+        msg="Minimum position with constant error 0.2 must be close to minimum position with constant error 0.1",
+    )
+    # Case 2: noise
+    synt_data[:, 1] += 0.05 * np.random.rand(1000) - 0.025
+    pos_noised, err_noised = calc_minpos_with_error(synt_data)
+    check.less(
         err_comp,
-        np.arcsin(0.1) / 8,
+        err_noised,
+        msg="Error of minimum with noise is not greater than error of minimum with constant error",
+    )
+    check.almost_equal(
+        pos_comp,
+        pos_noised,
         abs=1e-3,
-        msg="Error of minimum with constant error is not estimated correctly",
+        msg="Minimum position with and w/o noise should be similar",
     )