aare/python/tests/test_Minuit2_scurve_1d.ipynb

In [ ]:

import time
import numpy as np
np.set_printoptions(suppress=True, precision=6)

import matplotlib.pyplot as plt
from scipy.special import erf
from aare import RisingScurve, FallingScurve
from pprint import pprint

Model curves

In [ ]:

def scurve(x, p): # rising Scurve
    p0, p1, p2, p3, p4, p5 = p
    return (p0 + p1*x) + 0.5 *(1 + erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))

def scurve2(x, p): #falling Scurve
    p0, p1, p2, p3, p4, p5 = p
    return (p0 + p1*x) + 0.5 *(1 - erf((x-p2) / (np.sqrt(2)*p3))) * (p4 + p5*(x-p2))

Generate data (1D)

In [ ]:

x = np.linspace(0,120, 121)

rng = np.random.default_rng(42)
0
p_true_rising = np.array([100.0, 0.25, 60.0, 6.0, 120.0,  1.0])
p_true_falling = np.array([100.0, 0.25, 60.0, 6.0, 120.0, -1.0])

y_true_rising = scurve(x, p_true_rising)
y_true_falling = scurve2(x, p_true_falling)

noise_sigma = 4

y_rising = y_true_rising + rng.normal(0, noise_sigma, size=x.shape)
# y_err_r = np.full_like(x, noise_sigma)

y_falling = y_true_falling + rng.normal(0, noise_sigma, size=x.shape)
# y_err_f = np.full_like(x, noise_sigma)

Plot synthetic data

In [ ]:

fig, ax = plt.subplots(1, 2, figsize=(12,4), sharex=True)

ax[0].plot(x, y_true_rising, label="true")
ax[0].scatter(x, y_rising, s=12, label="noisy")
ax[0].set_title("Rising S-curve")
ax[0].legend()

ax[1].plot(x, y_true_falling, label="true")
ax[1].scatter(x, y_falling, s=12, label="noisy")
ax[1].set_title("Falling S-curve")
ax[1].legend()

plt.show()

Fit the Rising S-curve

In [ ]:

model_r = RisingScurve()
print("Paramter list of the Scurve:")
# print(model_r.GetParNames())
print(model_r.par_names)
print(model_r.n_par, '\n')

model_r.FixParameter('p0', 100)  # Fixed p0 = 100 (optimizer does not touch the value)
model_r.SetParameter('A', 100) # Start value A = 100

print("== Tuned fit settings ==")
model_r.compute_errors = True
model_r.max_calls = 500
model_r.tolerance = 0.01
print(f"max_calls : {model_r.max_calls}")
print(f"tolerance : {model_r.tolerance}")
print(f"compute_erros: {model_r.compute_errors}")
print("\n")
print("== Results ==")
res_r_munuit_obj = model_r.fit(x, y_rising, np.sqrt(y_rising))

print("True rising params:         ", p_true_rising)
print("minuit_grad (obj api ) rising result:")
pprint(res_r_munuit_obj, sort_dicts=False)

Fit the Falling S-curve

In [ ]:

model_f = FallingScurve()
model_f.SetParameter(2, 20) # set start p2 = 20
model_f.SetParameter(4, 90) # set start p4 = 90
model_f.SetParLimits(4, 80, 130)
model_f.max_calls = 150 # this is limits the number of calls done by the minimizer (increasing the `max_calls` may help if the fit fails) 
res_f_munuit_obj = model_f.fit(x, y_falling)

print("True falling params:  ", p_true_falling)
print("minuit_grad (obj api ) falling result")
pprint(res_f_munuit_obj, sort_dicts=False)

In [ ]:

fig, ax = plt.subplots(1, 2, figsize=(12,4))

# rising
ax[0].plot(x, y_true_rising, label="true")
ax[0].scatter(x, y_rising, s=12, label="data")
ax[0].plot(x, model_r(x, res_r_munuit_obj['par']), linewidth=1, color="green", label="minuit")
ax[0].set_title("Rising S-curve fits")
ax[0].legend()

# falling
ax[1].plot(x, y_true_falling, label="true")
ax[1].scatter(x, y_falling, s=12, label="data")
ax[1].plot(x, model_f(x, res_f_munuit_obj['par']), linewidth=1, color="green", label="minuit")
ax[1].set_title("Falling S-curve fits")
ax[1].legend()

plt.show()

Quick error check

In [ ]:

print("Rising absolute parameter error:", np.abs(res_r_munuit_obj['par'] - p_true_rising))
print("Falling absolute parameter error:", np.abs(res_f_munuit_obj['par'] - p_true_falling))

Benchmark

In [ ]:

def bench(fn, n_repeats=200):
    for _ in range(3):
        fn()

    t0 = time.perf_counter()
    for _ in range(n_repeats):
        res = fn()
    t1 = time.perf_counter()
    return res, (t1 - t0) / n_repeats

model_rising = RisingScurve()
model_falling = FallingScurve()

_, t_rising = bench(lambda: model_rising.fit(x, y_rising))
_, t_falling = bench(lambda: model_falling.fit(x, y_falling))

print(f"Rising Minuit2:  {1e3*t_rising:.3f} ms")
print(f"Falling Minuit2: {1e3*t_falling:.3f} ms")

In [ ]: