From b2b57df0d93b8c6414540cea8703e588d5a25048 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schl=C3=BC=C3=9Fler?= Date: Fri, 23 Sep 2022 13:27:05 +0200 Subject: [PATCH 1/4] Add quadruple fit function --- bmlab/fits.py | 128 ++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 124 insertions(+), 4 deletions(-) diff --git a/bmlab/fits.py b/bmlab/fits.py index 2b70e29..e38918e 100644 --- a/bmlab/fits.py +++ b/bmlab/fits.py @@ -141,6 +141,122 @@ def error(params, xdata, ydata): return w0s, fwhms, intens, offset +def fit_quadruple_lorentz(x, y, bounds_w0=None): + offset_guess = (y[0] + y[-1]) / 2. + # gam_guess = (w0_guess ** 2 * (np.max(y) - offset_guess)) ** -0.5 + fwhm_guess = 10 * (x[-1] - x[0]) / x.shape[0] + intensity_guess = np.max(y) - offset_guess + + # We run peak finding to get a good guess of the peak positions + # and use the four peaks with the highest prominence. + peaks, properties = find_peaks(y, prominence=1) + idx = np.argsort(properties['prominences'])[::-1] + + # Return if we didn't find four peaks + if len(idx) < 4: + return + + idx_sort = np.sort(peaks[idx[0:4]]) + w0_guess = list(x[idx_sort]) + + def error(params, xdata, ydata): + return (ydata + - lorentz(xdata, *params[0:3]) + - lorentz(xdata, *params[3:6]) + - lorentz(xdata, *params[6:9]) + - lorentz(xdata, *params[9:12]) + - params[12]) ** 2 + + # Create the bounds array + if bounds_w0 is not None: + # Lower limits + bounds_lower = -np.Inf * np.ones(13) + + # 1st peak: + # central position + bounds_lower[0] = bounds_w0[0][0] + # full-width-half-maximum + # The VIPA spectrometer has an instrument width of + # approx. 750 MHz for the FOB setup + # and 180 MHz for the 780 nm setup. + # This is far higher than the step size, + # so we limit it to this. + fwhm_lower_bound = (x[-1] - x[0]) / x.shape[0] + + bounds_lower[1] = fwhm_lower_bound + # intensity + bounds_lower[2] = 0 + + # 2nd peak: + # central position + bounds_lower[3] = bounds_w0[1][0] + # full-width-half-maximum + bounds_lower[4] = fwhm_lower_bound + # intensity + bounds_lower[5] = 0 + + # 3rd peak: + # central position + bounds_lower[6] = bounds_w0[2][0] + # full-width-half-maximum + bounds_lower[7] = fwhm_lower_bound + # intensity + bounds_lower[8] = 0 + + # 4th peak: + # central position + bounds_lower[9] = bounds_w0[3][0] + # full-width-half-maximum + bounds_lower[10] = fwhm_lower_bound + # intensity + bounds_lower[11] = 0 + + # offset + bounds_lower[12] = 0 + + # Upper limits + bounds_upper = np.Inf * np.ones(13) + + bounds_upper[0] = bounds_w0[0][1] + bounds_upper[3] = bounds_w0[1][1] + bounds_upper[6] = bounds_w0[2][1] + bounds_upper[9] = bounds_w0[3][1] + bounds = (bounds_lower, bounds_upper) + + # Sort the guesses to the bounds + bounds_w0_center = [np.mean( + np.clip(bound, *x[::len(x) - 1])) for bound in bounds_w0] + w0_guess.sort(reverse=(bounds_w0_center[0] > bounds_w0_center[1])) + + # Check that the initial guesses are within the bounds + w0_guess = [np.clip( + guess, *bounds_w0[idx]) for idx, guess in enumerate(w0_guess)] + else: + bounds = (-np.Inf, np.Inf) + + opt_result = least_squares( + error, + x0=(w0_guess[0], fwhm_guess, intensity_guess, + w0_guess[1], fwhm_guess, intensity_guess, + w0_guess[2], fwhm_guess, intensity_guess, + w0_guess[3], fwhm_guess, intensity_guess, + offset_guess + ), + args=(x, y), + bounds=bounds + ) + + if not opt_result.success: + raise FitError('Lorentz fit failed.') + + res = opt_result.x + w0s, fwhms, intens = (res[0], res[3], res[6], res[9]),\ + (res[1], res[4], res[7], res[10]),\ + (res[2], res[5], res[8], res[11]) + offset = res[12] + return w0s, fwhms, intens, offset + + def fit_circle(points): """ Fits a circle to a given set of points. Returnes the center and the radius @@ -224,13 +340,17 @@ def fit_lorentz_region(region, xdata, ydata, nr_peaks=1, bounds_w0=None): idx_r = np.nanargmin(np.abs(xdata - region[1])) x = xdata[idx_l:idx_r] y = ydata[idx_l:idx_r] - if nr_peaks == 2: + if nr_peaks == 1: + w0s, fwhms, intensities, offset = fit_lorentz( + x, y) + elif nr_peaks == 2: w0s, fwhms, intensities, offset = fit_double_lorentz( x, y, bounds_w0=bounds_w0) - elif nr_peaks == 1: - w0s, fwhms, intensities, offset = fit_lorentz( - x, y) + elif nr_peaks == 4: + w0s, fwhms, intensities, offset = fit_quadruple_lorentz( + x, y, + bounds_w0=bounds_w0) else: return except Exception: From a168083d66b23d8f4d341c52664b1bff2194d343 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schl=C3=BC=C3=9Fler?= Date: Mon, 26 Sep 2022 11:32:52 +0200 Subject: [PATCH 2/4] Sort Rayleigh to Brillouin peaks by minimizing Brillouin shift --- bmlab/controllers.py | 60 ++++++++++++----------------- tests/test_evaluation_controller.py | 43 +++++++-------------- 2 files changed, 38 insertions(+), 65 deletions(-) diff --git a/bmlab/controllers.py b/bmlab/controllers.py index 6e2dd35..7df3215 100644 --- a/bmlab/controllers.py +++ b/bmlab/controllers.py @@ -890,55 +890,43 @@ def get_indices_from_key(resolution, key): def calculate_derived_values(): """ - We calculate the derived parameters here: - - Brillouin shift [GHz] + We calculate the Brillouin shift in GHz here """ session = Session.get_instance() evm = session.evaluation_model() if not evm: return - if len(evm.results['time']) == 0: + if evm.results['brillouin_peak_position_f'].size == 0: + return + + if evm.results['rayleigh_peak_position_f'].size == 0: return shape_brillouin = evm.results['brillouin_peak_position_f'].shape shape_rayleigh = evm.results['rayleigh_peak_position_f'].shape - # If we have the same number of Rayleigh and Brillouin regions, - # we can simply subtract the two arrays (regions are always - # sorted by center in the peak selection model, so corresponding - # regions should be at the same array index) - if shape_brillouin[4] == shape_rayleigh[4]: - evm.results['brillouin_shift_f'] = abs( + + # We calculate every possible combination of + # Brillouin peak and Rayleigh peak position difference + # and then use the smallest absolute value. + # That saves us from sorting Rayleigh peaks to Brillouin peaks, + # because a Brillouin peak always belongs to the Rayleigh peak nearest. + brillouin_shift_f = np.NaN * np.ones((*shape_brillouin, shape_rayleigh[4])) + for idx in range(shape_rayleigh[4]): + brillouin_shift_f[:, :, :, :, :, :, idx] = abs( evm.results['brillouin_peak_position_f'] - - evm.results['rayleigh_peak_position_f'] + np.tile( + evm.results['rayleigh_peak_position_f'][:, :, :, :, [idx], :], + (1, 1, 1, 1, 1, shape_brillouin[5]) + ) ) - # Having a different number of Rayleigh and Brillouin regions - # doesn't really make sense. But in case I am missing something - # here, we assign each Brillouin region the nearest (by center) - # Rayleigh region. - else: - psm = session.peak_selection_model() - if not psm: - return - brillouin_centers = list(map(np.mean, psm.get_brillouin_regions())) - rayleigh_centers = list(map(np.mean, psm.get_rayleigh_regions())) - - for idx in range(len(brillouin_centers)): - # Find corresponding (nearest) Rayleigh region - d = list( - map( - lambda x: abs(x - brillouin_centers[idx]), - rayleigh_centers - ) - ) - idx_r = d.index(min(d)) - evm.results['brillouin_shift_f'][:, :, :, :, idx, :] = abs( - evm.results[ - 'brillouin_peak_position_f'][:, :, :, :, idx, :] - - evm.results[ - 'rayleigh_peak_position_f'][:, :, :, :, idx_r, :] - ) + with warnings.catch_warnings(): + warnings.filterwarnings( + action='ignore', + message='All-NaN slice encountered' + ) + evm.results['brillouin_shift_f'] = np.nanmin(brillouin_shift_f, 6) class Controller(object): diff --git a/tests/test_evaluation_controller.py b/tests/test_evaluation_controller.py index 1f3c5ac..876e586 100644 --- a/tests/test_evaluation_controller.py +++ b/tests/test_evaluation_controller.py @@ -118,7 +118,7 @@ def test_calculate_derived_values_equal_region_count(): evm.results['brillouin_peak_position_f'][:, :, :, :, 0, :] = 1 evm.results['brillouin_peak_position_f'][:, :, :, :, 1, :] = 4 - evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = 3 + evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = -1 evm.results['rayleigh_peak_position_f'][:, :, :, :, 1, :] = 8 calculate_derived_values() @@ -159,14 +159,14 @@ def test_calculate_derived_values_equal_region_count_nr_peaks_2(): evm.results['brillouin_peak_position_f'][:, :, :, :, 1, :] = 4 evm.results['brillouin_peak_position_f'][:, :, :, :, 0, 1] = 2 evm.results['brillouin_peak_position_f'][:, :, :, :, 1, 1] = 5 - evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = 3 + evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = -1 evm.results['rayleigh_peak_position_f'][:, :, :, :, 1, :] = 8 calculate_derived_values() assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 0] == 2).all() assert (evm.results['brillouin_shift_f'][:, :, :, :, 1, 0] == 4).all() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 1] == 1).all() + assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 1] == 3).all() assert (evm.results['brillouin_shift_f'][:, :, :, :, 1, 1] == 3).all() @@ -193,22 +193,13 @@ def test_calculate_derived_values_different_region_count(): evm.results['brillouin_peak_position_f'][:, :, :, :, 0, :] = 1 evm.results['brillouin_peak_position_f'][:, :, :, :, 1, :] = 4 - evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = 5 + evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = -2 evm.results['rayleigh_peak_position_f'][:, :, :, :, 1, :] = 9 evm.results['rayleigh_peak_position_f'][:, :, :, :, 2, :] = 10 - psm = evc.session.peak_selection_model() - - psm.add_brillouin_region((1, 2)) - psm.add_brillouin_region((4, 5)) - - psm.add_rayleigh_region((0, 1)) - psm.add_rayleigh_region((6, 7)) - psm.add_rayleigh_region((8, 9)) - calculate_derived_values() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, :] == 4).all() + assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, :] == 3).all() assert (evm.results['brillouin_shift_f'][:, :, :, :, 1, :] == 5).all() @@ -237,25 +228,19 @@ def test_calculate_derived_values_different_region_count_nr_peaks_2(): evm.results['brillouin_peak_position_f'][:, :, :, :, 1, :] = 4 evm.results['brillouin_peak_position_f'][:, :, :, :, 0, 1] = 2 evm.results['brillouin_peak_position_f'][:, :, :, :, 1, 1] = 5 - evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = 5 + evm.results['rayleigh_peak_position_f'][:, :, :, :, 0, :] = -2 evm.results['rayleigh_peak_position_f'][:, :, :, :, 1, :] = 9 - evm.results['rayleigh_peak_position_f'][:, :, :, :, 2, :] = 10 - - psm = evc.session.peak_selection_model() - - psm.add_brillouin_region((1, 2)) - psm.add_brillouin_region((4, 5)) - - psm.add_rayleigh_region((0, 1)) - psm.add_rayleigh_region((6, 7)) - psm.add_rayleigh_region((8, 9)) + evm.results['rayleigh_peak_position_f'][0:4, :, :, :, 2, :] = 10 + evm.results['rayleigh_peak_position_f'][4, :, :, :, 2, :] = 8 calculate_derived_values() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 0] == 4).all() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 1, 0] == 5).all() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 1] == 3).all() - assert (evm.results['brillouin_shift_f'][:, :, :, :, 1, 1] == 4).all() + assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 0] == 3).all() + assert (evm.results['brillouin_shift_f'][0:4, :, :, :, 1, 0] == 5).all() + assert (evm.results['brillouin_shift_f'][4, :, :, :, 1, 0] == 4).all() + assert (evm.results['brillouin_shift_f'][:, :, :, :, 0, 1] == 4).all() + assert (evm.results['brillouin_shift_f'][0:4, :, :, :, 1, 1] == 4).all() + assert (evm.results['brillouin_shift_f'][4, :, :, :, 1, 1] == 3).all() def test_get_data_0D(mocker): From 41460dd8c52af33adc2cc193d597031ae3f9b494 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schl=C3=BC=C3=9Fler?= Date: Tue, 27 Sep 2022 12:22:21 +0200 Subject: [PATCH 3/4] Adjust bounds calculation for 4-peak-fit --- bmlab/controllers.py | 48 +++++++++++++----- tests/test_evaluation_controller.py | 75 +++++++++++++++++++++++++++++ 2 files changed, 112 insertions(+), 11 deletions(-) diff --git a/bmlab/controllers.py b/bmlab/controllers.py index 7df3215..c5a54be 100644 --- a/bmlab/controllers.py +++ b/bmlab/controllers.py @@ -703,24 +703,49 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks): if bounds is None: return None - # We need a Rayleigh peak position for - # every brillouin_region - if rayleigh_peaks.shape[0] != (len(brillouin_regions)): + # We need two Rayleigh peak positions to determine the FSR + # and decide whether a peak is Stokes or Anti-Stokes + if rayleigh_peaks.shape[0] != 2: return None w0_bounds = [] # We have to create a separate bound for every region for region_idx, region in enumerate(brillouin_regions): local_time = [] - for f_rayleigh in rayleigh_peaks[region_idx, :]: - # In case this is an Anti-Stokes peak, we find the peak - # with the higher frequency on the left hand side and - # have to flip the bounds - is_anti_stokes = np.mean(region) <\ - f_rayleigh + for rayleigh_idx in range(rayleigh_peaks.shape[1]): + anti_stokes_limit = np.nanmean(rayleigh_peaks[:, rayleigh_idx]) + # We need to figure out, whether a peak is a + # Stokes or Anti-Stokes peak in order to decide + # to which Rayleigh peak we need to relate to. + tmp = region >= anti_stokes_limit + is_pure_region = (tmp == tmp[0]).all() + + is_anti_stokes_region = np.mean(region) >\ + anti_stokes_limit local_bound = [] for bound in bounds: + parsed_bound = [] + for limit in bound: + try: + parsed_bound.append(float(limit)) + except ValueError: + parsed_bound.append(np.NaN) + # We don't treat Inf as a value + with warnings.catch_warnings(): + warnings.filterwarnings( + action='ignore', + message='Mean of empty slice' + ) + is_anti_stokes_peak =\ + np.nanmean( + np.array(parsed_bound)[ + np.isfinite(parsed_bound)] + ) < 0 + + is_anti_stokes = is_anti_stokes_region if\ + is_pure_region else is_anti_stokes_peak + local_limit = [] for limit in bound: if limit.lower() == 'min': @@ -739,10 +764,11 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks): # a value in pixel depending on the time try: val = ((-1) ** is_anti_stokes)\ - * 1e9 * float(limit) + f_rayleigh + * 1e9 * abs(float(limit))\ + + rayleigh_peaks[ + int(is_anti_stokes), rayleigh_idx] except BaseException: val = np.Inf - # print(val) local_limit.append(val) # Check that the bounds are sorted ascendingly # (for anti-stokes, they might not). diff --git a/tests/test_evaluation_controller.py b/tests/test_evaluation_controller.py index 876e586..22bb453 100644 --- a/tests/test_evaluation_controller.py +++ b/tests/test_evaluation_controller.py @@ -635,6 +635,8 @@ def test_create_bounds(mocker): fit_bounds = EvaluationController().create_bounds(None, None) assert fit_bounds is None + # Test that we get correct bounds for regions + # that only contain one type of peaks (either Stokes or Anti-Stokes) evm.bounds = [['min', '5'], ['5.5', 'Inf'], ['-inf', 'max']] brillouin_regions = [(3.3e9, 7.0e9), (8.0e9, 12.1e9)] @@ -652,3 +654,76 @@ def test_create_bounds(mocker): [[10.3e9, 12.1e9], [-np.inf, 9.8e9], [8.0e9, np.inf]] ] ], fit_bounds, atol=1e6) + + # Test that we get the correct bounds for regions + # that contain both Stokes and Anti-Stokes peaks + evm.bounds = [['min', '5'], ['5.5', 'Inf'], ['-inf', 'max'], + ['min', '-5'], ['-5.5', 'Inf'], ['-inf', 'max']] + + brillouin_regions = [(2.3e9, 12.1e9), (3.3e9, 12.1e9)] + rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]]) + fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks) + + np.testing.assert_allclose([ + [ + [[2.3e9, 5.0e9], [5.5e9, np.inf], [-np.inf, 12.1e9], + [9.8e9, 12.1e9], [-np.inf, 9.3e9], [-np.inf, 12.1e9]], + [[2.3e9, 5.1e9], [5.6e9, np.inf], [-np.inf, 12.1e9], + [10.0e9, 12.1e9], [-np.inf, 9.5e9], [-np.inf, 12.1e9]], + [[2.3e9, 4.8e9], [5.3e9, np.inf], [-np.inf, 12.1e9], + [10.3e9, 12.1e9], [-np.inf, 9.8e9], [-np.inf, 12.1e9]] + ], + [ + [[3.3e9, 5.0e9], [5.5e9, np.inf], [-np.inf, 12.1e9], + [9.8e9, 12.1e9], [-np.inf, 9.3e9], [-np.inf, 12.1e9]], + [[3.3e9, 5.1e9], [5.6e9, np.inf], [-np.inf, 12.1e9], + [10.0e9, 12.1e9], [-np.inf, 9.5e9], [-np.inf, 12.1e9]], + [[3.3e9, 4.8e9], [5.3e9, np.inf], [-np.inf, 12.1e9], + [10.3e9, 12.1e9], [-np.inf, 9.8e9], [-np.inf, 12.1e9]] + ] + ], fit_bounds, atol=1e6) + + # Test real values for a two peak fit + evm.bounds = [['3.5', '5'], ['5.0', '7.5']] + + brillouin_regions = [(3.3e9, 7.5e9), (8.0e9, 12.1e9)] + rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]]) + fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks) + + np.testing.assert_allclose([ + [ + [[3.5e9, 5.0e9], [5.0e9, 7.5e9]], + [[3.6e9, 5.1e9], [5.1e9, 7.6e9]], + [[3.3e9, 4.8e9], [4.8e9, 7.3e9]] + ], [ + [[9.8e9, 11.3e9], [7.3e9, 9.8e9]], + [[10.0e9, 11.5e9], [7.5e9, 10.0e9]], + [[10.3e9, 11.8e9], [7.8e9, 10.3e9]] + ] + ], fit_bounds, atol=1e6) + + # Test real values for a four peak fit over the whole spectrum + evm.bounds = [['3.5', '5'], ['5.0', '7.5'], + ['-7.5', '-5.0'], ['-5', '-3.5']] + + brillouin_regions = [(2.3e9, 12.1e9), (3.3e9, 12.1e9)] + rayleigh_peaks = 1e9 * np.array([[0.0, 0.1, -0.2], [14.8, 15.0, 15.3]]) + fit_bounds = evc.create_bounds(brillouin_regions, rayleigh_peaks) + + np.testing.assert_allclose([ + [ + [[3.5e9, 5.0e9], [5.0e9, 7.5e9], + [7.3e9, 9.8e9], [9.8e9, 11.3e9]], + [[3.6e9, 5.1e9], [5.1e9, 7.6e9], + [7.5e9, 10.0e9], [10.0e9, 11.5e9]], + [[3.3e9, 4.8e9], [4.8e9, 7.3e9], + [7.8e9, 10.3e9], [10.3e9, 11.8e9]] + ], [ + [[3.5e9, 5.0e9], [5.0e9, 7.5e9], + [7.3e9, 9.8e9], [9.8e9, 11.3e9]], + [[3.6e9, 5.1e9], [5.1e9, 7.6e9], + [7.5e9, 10.0e9], [10.0e9, 11.5e9]], + [[3.3e9, 4.8e9], [4.8e9, 7.3e9], + [7.8e9, 10.3e9], [10.3e9, 11.8e9]] + ] + ], fit_bounds, atol=1e6) From 7315f1d3793bd58b07012a99440000e0e5cf3913 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Raimund=20Schl=C3=BC=C3=9Fler?= Date: Wed, 28 Sep 2022 13:38:12 +0200 Subject: [PATCH 4/4] Add test for quadruple lorentz fit w/o bounds --- tests/test_fits.py | 51 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/tests/test_fits.py b/tests/test_fits.py index 84f6d34..6267d12 100644 --- a/tests/test_fits.py +++ b/tests/test_fits.py @@ -4,7 +4,7 @@ from bmlab.fits import lorentz, fit_lorentz, fit_circle,\ fit_double_lorentz, calculate_exact_circle, fit_vipa, VIPA,\ - are_points_on_line + are_points_on_line, fit_quadruple_lorentz from bmlab.models.setup import AVAILABLE_SETUPS @@ -75,6 +75,55 @@ def test_fit_double_lorentz(): np.testing.assert_almost_equal(actual_offset, offset, decimal=3) +def test_fit_quadruple_lorentz(): + # Arrange + x = np.linspace(0, 60, 200) + + w0_0 = 10. + intensity_0 = 8. + fwhm_0 = 2. + + w0_1 = 20. + intensity_1 = 12. + fwhm_1 = 4. + + w0_2 = 30. + intensity_2 = 14. + fwhm_2 = 4. + + w0_3 = 40. + intensity_3 = 16. + fwhm_3 = 4. + + offset = 10. + + y_data = lorentz(x, w0_0, fwhm_0, intensity_0) + y_data += lorentz(x, w0_1, fwhm_1, intensity_1) + y_data += lorentz(x, w0_2, fwhm_2, intensity_2) + y_data += lorentz(x, w0_3, fwhm_3, intensity_3) + offset + + w0s, fwhms, intens, actual_offset\ + = fit_quadruple_lorentz(x, y_data) + + np.testing.assert_almost_equal(w0s[0], w0_0, decimal=3) + np.testing.assert_almost_equal(fwhms[0], fwhm_0, decimal=3) + np.testing.assert_almost_equal(intens[0], intensity_0, decimal=3) + + np.testing.assert_almost_equal(w0s[1], w0_1, decimal=3) + np.testing.assert_almost_equal(fwhms[1], fwhm_1, decimal=3) + np.testing.assert_almost_equal(intens[1], intensity_1, decimal=3) + + np.testing.assert_almost_equal(w0s[2], w0_2, decimal=3) + np.testing.assert_almost_equal(fwhms[2], fwhm_2, decimal=3) + np.testing.assert_almost_equal(intens[2], intensity_2, decimal=3) + + np.testing.assert_almost_equal(w0s[3], w0_3, decimal=3) + np.testing.assert_almost_equal(fwhms[3], fwhm_3, decimal=3) + np.testing.assert_almost_equal(intens[3], intensity_3, decimal=3) + + np.testing.assert_almost_equal(actual_offset, offset, decimal=3) + + def test_fit_double_lorentz_with_bounds(): # Arrange x = np.linspace(0, 60, 500)