Add in x scaling prior to fit

pyquant/peaks.py

@@ -262,7 +262,7 @@ def findEnvelope(xdata, ydata, measured_mz=None, theo_mz=None, max_mz=None, prec
 
     return {'envelope': env_dict, 'micro_envelopes': micro_dict, 'ppms': ppm_dict}
 
-def findAllPeaks(xdata, ydata_original, min_dist=0, method=None, local_filter_size=0, filter=False, peak_boost=False, bigauss_fit=False,
+def findAllPeaks(xdata_original, ydata_original, min_dist=0, method=None, local_filter_size=0, filter=False, peak_boost=False, bigauss_fit=False,
                  rt_peak=None, mrm=False, max_peaks=4, debug=False, peak_width_start=3, snr=0, zscore=0, amplitude_filter=0,
                  peak_width_end=4, fit_baseline=False, rescale=True, fit_negative=False, percentile_filter=0, micro=False,
                  method_opts=None, smooth=False, r2_cutoff=None, peak_find_method=PEAK_FINDING_REL_MAX, min_slope=None,
@@ -276,10 +276,14 @@ def findAllPeaks(xdata, ydata_original, min_dist=0, method=None, local_filter_si
     if micro:
         fit_baseline = False
 
+    xdata = np.copy(xdata_original)
+
     if not micro and gap_interpolation:
         ydata_original = interpolate_data(xdata, ydata_original, gap_limit=gap_interpolation)
 
     rel_peak_constraint = (0.0 if fit_baseline else 0.5)
+    original_x_max = np.max(xdata)
+    xdata /= original_x_max
     original_max = np.abs(ydata_original).max() if fit_negative else ydata_original.max()
     amplitude_filter /= original_max
     ydata = ydata_original / original_max
@@ -291,6 +295,8 @@ def findAllPeaks(xdata, ydata_original, min_dist=0, method=None, local_filter_si
     if smooth and len(ydata) > 5:
         ydata_peaks = savgol_smooth(ydata_peaks)
 
+    if rt_peak is not None:
+        rt_peak /= original_x_max
 
     if filter or peak_boost:
         if len(ydata) >= 5:
@@ -455,7 +461,12 @@ def findAllPeaks(xdata, ydata_original, min_dist=0, method=None, local_filter_si
                 jacobian = bigauss_jac if bigauss_fit else gauss_jac
 
             if debug:
-                print('left and right segments', xdata[left_break_point], xdata[right_break_point-1])
+                print('left and right segments {}({}), {}({})'.format(
+                    xdata[left_break_point],
+                    xdata[left_break_point]*original_x_max,
+                    xdata[right_break_point-1],
+                    xdata[right_break_point-1]*original_x_max,
+                ))
                 print('guess and bnds', segment_guess, segment_bounds)
             hessian = None  # if bigauss_fit else gauss_hess
 
@@ -618,6 +629,17 @@ def findAllPeaks(xdata, ydata_original, min_dist=0, method=None, local_filter_si
             # Intercept
             best_fit[step_size-1::step_size] *= original_max
 
+    # rescale the x axis back up
+    best_fit[1::step_size] *= original_x_max
+    best_fit[2::step_size] *= original_x_max
+    if bigauss_fit:
+        best_fit[3::step_size] *= original_x_max
+
+    if fit_baseline:
+        # Slope
+        best_fit[step_size - 2::step_size] *= original_max
+
+
     return best_fit, residual
 
 

pyquant/tests/test_peaks.py

@@ -111,15 +111,31 @@ def test_experimental(self):
         # Experimental data
         x, y = self.peak_data['offset_fit']
         params, residual = peaks.findAllPeaks(x, y, bigauss_fit=True, filter=True)
+        expected = np.array([
+          13219262.587656807, 46.821340819991505, 0.06523272363478014, 0.18374422913588656,
+          3347200.309180678, 47.497, 0.6166821402545103, 0.3817876338966981,
+          1880722.1582992678, 48.14756707645761, 0.17537885391522443, 0.4846763157315077,
+          1473766.5256005626, 49.52607160264086, 0.020199999999999108, 0.22250905781532157
+        ])
+        np.testing.assert_allclose(
+          params[::4],
+          expected[::4],
+          atol=10000,
+        )
+        np.testing.assert_allclose(
+          params[1::4],
+          expected[1::4],
+          atol=1,
+        )
+        np.testing.assert_allclose(
+          params[2::4],
+          expected[2::4],
+          atol=0.1,
+        )
         np.testing.assert_allclose(
-            params,
-            np.array([
-                13219262.587656807, 46.821340819991505, 0.06523272363478014, 0.18374422913588656,
-                3347200.309180678, 47.497, 0.6166821402545103, 0.3817876338966981,
-                1880722.1582992678, 48.14756707645761, 0.17537885391522443, 0.4846763157315077,
-                1473766.5256005626, 49.52607160264086, 0.020199999999999108, 0.22250905781532157
-            ]),
-            atol=10,
+          params[3::4],
+          expected[3::4],
+          atol=0.1,
         )
 
 if __name__ == '__main__':

pyquant/tests/test_utils.py

@@ -142,7 +142,7 @@ def test_get_scan_resolution(self):
         x, y = data['low_res_scan']
         scan = pd.Series(y, index=x)
         resolution = utils.get_scan_resolution(scan)
-        self.assertAlmostEqual(int(resolution), 30459)
+        self.assertAlmostEqual(int(resolution), 30459, delta=1000)
 
         x = np.array([
             805.47264226,  805.47265226,  805.47495304,  805.47739824,