Skip to content

Commit

Permalink
calibration: generalize param vector format
Browse files Browse the repository at this point in the history
  • Loading branch information
JoepVanlier committed Oct 25, 2023
1 parent 9912f93 commit 292f6a0
Showing 1 changed file with 7 additions and 10 deletions.
17 changes: 7 additions & 10 deletions lumicks/pylake/force_calibration/power_spectrum_calibration.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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,
Expand Down

0 comments on commit 292f6a0

Please sign in to comment.