From 292f6a0619c807c77a71434529f5a839a9b9fec3 Mon Sep 17 00:00:00 2001 From: JoepVanlier Date: Wed, 25 Oct 2023 11:11:34 +0200 Subject: [PATCH] calibration: generalize param vector format --- .../power_spectrum_calibration.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/lumicks/pylake/force_calibration/power_spectrum_calibration.py b/lumicks/pylake/force_calibration/power_spectrum_calibration.py index f5bc4d04d..ce30338e9 100644 --- a/lumicks/pylake/force_calibration/power_spectrum_calibration.py +++ b/lumicks/pylake/force_calibration/power_spectrum_calibration.py @@ -177,12 +177,11 @@ def calculate_power_spectrum( def lorentzian_loss(p, model: ScaledModel, power_spectrum: PowerSpectrum): - expectation = model(power_spectrum.frequency, p[0], p[1], *p[2:]) + expectation = model(power_spectrum.frequency, *p) gamma = expectation / power_spectrum.num_points_per_block**0.5 # fmt: off return np.sum(np.log( - 1 + 0.5 * ((power_spectrum.power - model(power_spectrum.frequency, p[0], p[1], *p[2:])) - / gamma) ** 2 + 1 + 0.5 * ((power_spectrum.power - model(power_spectrum.frequency, *p)) / gamma) ** 2 )) # fmt: on @@ -296,20 +295,18 @@ def fit_power_spectrum( # The actual curve fitting process is driven by a set of fit parameters that are of order unity. # This increases the robustness of the fit (see ref. 3). The `ScaledModel` model class takes # care of this parameter rescaling. - scaled_model = ScaledModel( - lambda f, fc, D, *filter_pars: 1 / model(f, fc, D, *filter_pars), initial_params + scaled_model = ScaledModel(lambda f, *params: 1 / model(f, *params), initial_params) + lower_bounds = scaled_model.normalize_params([0.0, 0.0, *model._filter.lower_bounds()]) + upper_bounds = scaled_model.normalize_params( + [np.inf, np.inf, *model._filter.upper_bounds(power_spectrum.sample_rate)] ) - + # What we *actually* have to minimize, is the chi^2 expression in Eq. 39 of ref. 1. We're # "hacking" `scipy.optimize.curve_fit` for this purpose, and passing in "y" values of "1/ps.power" # instead of just "ps.power", and a model function "1/model(...)" instead of "model(...)". This # effectively transforms the curve fitter's objective function # "np.sum( ((f(xdata, *popt) - ydata) / sigma)**2 )" into the expression in Eq. 39 of ref. 1. sigma = (1 / power_spectrum.power) / math.sqrt(power_spectrum.num_points_per_block) - lower_bounds = scaled_model.normalize_params([0.0, 0.0, *model._filter.lower_bounds()]) - upper_bounds = scaled_model.normalize_params( - [np.inf, np.inf, *model._filter.upper_bounds(power_spectrum.sample_rate)] - ) (solution_params_rescaled, pcov) = scipy.optimize.curve_fit( scaled_model, power_spectrum.frequency,