diff --git a/pyquant/peaks.py b/pyquant/peaks.py index 6b1e70c..9661bd7 100644 --- a/pyquant/peaks.py +++ b/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 diff --git a/pyquant/tests/test_peaks.py b/pyquant/tests/test_peaks.py index ebfe211..d0fe666 100644 --- a/pyquant/tests/test_peaks.py +++ b/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__': diff --git a/pyquant/tests/test_utils.py b/pyquant/tests/test_utils.py index 5f4adaa..d06cc95 100644 --- a/pyquant/tests/test_utils.py +++ b/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,