Fitting without rebuilding¶
Build the model once with the chained methods. During a fit, never rebuild it: use the flat parameter vector, and bind the energy grid so unchanged peaks are not recomputed.
The parameter vector¶
m = (Model()
.gaussian(933.8, 0.49, 4.1e-5)
.voigt(938.3, 0.7, 0.2, 5e-5)
.step(939.4, 6.1, 7.2e-5)
.continuum(3.6e-4, -3.7e-7))
m.param_names # ['gaussian[0].e0', 'gaussian[0].width', 'gaussian[0].amp', 'voigt[0].e0', ..., 'scale', 'shift']
theta = m.params # flat float64 array
Layout: each component's parameters in the order you added it, then continuum.bkg, continuum.slope,
scale, shift. A pseudo-Voigt is one component, so eta is a fittable parameter.
| operation | effect |
|---|---|
m(E, params=theta) |
evaluate at theta; the model is unchanged |
m.set_params(theta) |
update the model in place (validated: bad values raise ValueError, model unchanged) |
Bind a grid¶
ev = m.bind(E) # fix the energy grid once
n = ev(theta) # set parameters and evaluate
n = ev(theta, out=buf) # same, into a preallocated complex128 array
The evaluator caches each component's contribution. A call recomputes only the components whose parameters changed since the previous call:
scale,continuum.bkgandcontinuum.slopecost no peak work at all.- Changing one peak costs one peak, not all of them. This is exactly the access pattern of a finite-difference Jacobian, which perturbs one parameter at a time.
shiftmoves every component, so it recomputes all of them.
Use it with scipy¶
import numpy as np
from scipy.optimize import least_squares
ev = model.bind(E)
buf = np.empty(E.size, complex)
theta = model.params.copy()
free = np.array([n.startswith("gaussian") or n == "scale" for n in model.param_names])
def residual(p):
theta[free] = p
return ev(theta, out=buf).imag - beta # fit beta; delta comes along for free
sol = least_squares(residual, theta[free], x_scale="jac")
The complete version, with bounds and DFT-seeded starts, is in the ZnPc example.
Freezing a parameter is just leaving it out of free.
How much does it save?¶
Forward finite-difference Jacobian, one parameter changed per column, Gaussians plus two steps, a delta table and a scale, on 4001 energies (Apple Silicon, single core):
| peaks | rebuild the model per column | bound evaluator (cached) | speedup |
|---|---|---|---|
| 20 | 22.7 ms | 2.6 ms | 9× |
| 50 | 106 ms | 5.8 ms | 18× |
| 100 | 498 ms | 12.0 ms | 41× |
A full optimizer step where every parameter changes costs the same either way: the time is the peaks
themselves (about 0.6 ms for 50 peaks on 4001 energies). Reproduce with
python -m benchmarks.no_rebuild.
Evaluation speed¶
- Energies sorted ascending take a faster path: far from each peak the tails use asymptotic series, and the delta table is walked with a moving pointer.
m(E, parallel=True)uses all cores and returns identical results. It already helps at a few thousand energies (the Cu model on 4001 points: about 1.6× faster) and approaches linear scaling on large grids (400,000 points: 12.8 ms down to 3.2 ms on this machine).- The Rust crate exposes the same model, parameter vector and evaluator; see the Rust API.
Analytic Jacobians
Each peak's derivatives with respect to position, width and height have closed forms, which would remove the per-column cost entirely. Not implemented yet; the evaluator is the supported path today.