Inverse parameter fitting#

PyMieSim includes a small NumPy-only, derivative-free fitting layer. It adds no required optimization dependency. The forward model remains a normal PyMieSim Simulation or Experiment; the fitting layer only supplies parameters, collects predictions, and compares them with observations.

FitParameter defines a named initial value and mandatory bounds. MeasuredData holds measured values and optional positive standard uncertainties. The user supplies a callable receiving a mapping of parameter names to values and returning a prediction with the same shape as the observation.

The earlier names Parameter, Observation, and FitResult remain available as aliases for FitParameter, MeasuredData, and ParameterFitResult, respectively.

The default coordinate search is deterministic and intended for a small number of parameters. It is derivative-free and has no global-optimum guarantee; use physical bounds and repeat fits from several initial guesses for multimodal resonant problems.

Set show_progress=True to display the iteration, current objective, and normalized search step while a fit is running. Progress output is disabled by default, which keeps scripts and notebooks quiet.

ParameterFitResult contains fitted parameters, predictions, residuals, objective value, iteration and evaluation counts, and convergence status. Unit-bearing parameters are passed to the model with their original units, and uncertainties are used to weight residuals.

Quantity bounds are converted to the initial parameter’s units. Quantity predictions and uncertainties are converted to the observation’s units before computing residuals; incompatible dimensions raise a Pint DimensionalityError. Bare numbers are interpreted as magnitudes in the corresponding reference units. When the reference is a bare number, quantities must be dimensionless. Zero-valued initial quantities are supported.

ParameterFitResult.prediction and ParameterFitResult.observed are numerical arrays in the observation’s units. ParameterFitResult.residuals contains their numerical difference, divided by the uncertainty when one is supplied.

The fitted values are conditional on the forward model. A fit to a dense suspension using isolated-particle Mie calculations can compensate for missing multiple scattering or particle correlations. For correlation-aware scattering models, refer to PackLab’s scattering workflow.

Fit a noisy scattering spectrum#

The complete script below estimates diameter and refractive index from nine synthetic wavelength channels. It starts with a known sphere, adds Gaussian noise with a fixed seed, then fits a second sphere model. Wavelengths and diameter carry length units; Qsca and its uncertainties are dimensionless. No material catalog or external measurements are needed.

scattering_spectrum adapts single-particle simulations to an array of predictions with shape (9,).
MeasuredData pairs those observations with standard uncertainties of the same shape. One deliberately biased channel is assigned a larger uncertainty to demonstrate weighting.
FitParameter supplies physical bounds and initial guesses. The model receives diameter as a quantity and refractive index as a plain number.
The overlay compares predictions with measurements; the lower panel exposes systematic differences that a single objective value can hide.

import matplotlib.pyplot as plt
import numpy as np

from PyMieSim import Measure, FitParameter, MeasuredData, PlaneWave, PolarizationState, Simulation, Sphere, fit_parameters, ureg

wavelengths = np.linspace(450, 750, 9) * ureg.nanometer


def scattering_spectrum(diameter, refractive_index):
    """Return isolated-sphere scattering efficiencies at all wavelengths."""
    spectrum = np.empty(wavelengths.size)
    for index, wavelength in enumerate(wavelengths):
        source = PlaneWave(
            wavelength=wavelength,
            polarization=PolarizationState(angle=0 * ureg.degree),
            amplitude=1 * ureg.volt / ureg.meter,
        )
        simulation = Simulation(
            scatterer=Sphere(
                diameter=diameter,
                material=refractive_index,
                medium=1.0,
            ),
            source=source,
        )
        spectrum[index] = simulation.get(Measure.QSCA).magnitude
    return spectrum


observed = scattering_spectrum(
    diameter=310 * ureg.nanometer,
    refractive_index=1.58,
)
uncertainty = np.full_like(observed, 0.015)
uncertainty[5] = 0.08
rng = np.random.default_rng(seed=42)
observed = observed + rng.normal(loc=0.0, scale=uncertainty)
observed[5] += 0.08


def model(parameters):
    return scattering_spectrum(
        diameter=parameters["diameter"],
        refractive_index=parameters["refractive_index"],
    )


result = fit_parameters(
    model=model,
    observation=MeasuredData(values=observed, uncertainty=uncertainty),
    parameters=[
        FitParameter(
            name="diameter",
            initial=240 * ureg.nanometer,
            bounds=(150 * ureg.nanometer, 450 * ureg.nanometer),
        ),
        FitParameter(name="refractive_index", initial=1.45, bounds=(1.2, 2.0)),
    ],
)

print(result.summary())
print(result.parameters)
print("Standardized residual shape:", result.residuals.shape)  # (9,)

wavelength_values = wavelengths.to("nanometer").magnitude
figure, axes = plt.subplots(2, 1, figsize=(8, 6.5), sharex=True, height_ratios=(2, 1))
axes[0].errorbar(wavelength_values, observed, yerr=uncertainty, fmt="o", label="Synthetic observation")
axes[0].plot(wavelength_values, result.prediction, linewidth=2, label="Fitted isolated-sphere model")
axes[0].set(ylabel="Qsca", title="Multi-wavelength fitting constrains a spectral response")
axes[0].legend()
axes[1].axhline(0, color="black", linewidth=1)
axes[1].plot(wavelength_values, result.residuals, "o-")
axes[1].set(xlabel="Vacuum wavelength [nm]", ylabel="Standardized residual")
for axis in axes:
    axis.grid(alpha=0.2)
plt.show()

The minimized objective is the sum of squared standardized residuals:

\[r_i=\frac{f_i(\mathbf{p})-y_i}{\sigma_i},\qquad J(\mathbf{p})=\sum_i r_i^2.\]

Without uncertainties, residuals are numerical prediction-minus-observation values in the observation’s units. objective is the raw sum, not a mean, reduced chi-square, or parameter uncertainty estimate. Each sigma must be positive; correlations between measurement channels are not modeled.

Interpreting convergence and ambiguity#

For this script, expect a fitted spectrum near the noisy observations and parameters near the generating values (310 nm and 1.58), rather than exact recovery. There are nine predictions and nine residuals. A positive residual means the model lies above the observation; the deliberately biased channel can retain a negative residual even after convergence.

success=True means the coordinate search reached its step tolerance. It does not establish a global optimum, a unique physical solution, or a good model. Check message, objective, and evaluations as well as the plot. If the fit stops at its iteration limit, inspect the trajectory and bounds before increasing max_iterations.

Repeat the fit from several diameter and refractive-index starting points within the same bounds. Compare objectives and predictions: different parameter pairs with similar objectives indicate ambiguity. A parameter close to a bound can indicate an active physical constraint or an overly restrictive search range. Plot an objective grid across the two parameters to inspect valleys and competing minima; add informative wavelength channels when the data permit it. The returned result does not provide confidence intervals or a covariance matrix.

A visually good fit remains conditional on the isolated-particle forward model. Use residual structure to investigate calibration, material dispersion, or missing physics before interpreting fitted values as material properties. See the rendered diagnostic example for its plots.

Smaller examples#

Fit a scalar scale factor introduces the minimal API.
Fit a unit-aware diameter demonstrates quantity parameters.
Fit an isolated-sphere scattering measurement fits one sphere diameter.

API#

class FitParameter#

Bases: pybind11_object

A named bounded scalar parameter; values may carry Pint units.

property bounds#

property initial#

property name#

class MeasuredData#

Bases: pybind11_object

Measured values and optional positive standard uncertainties.

property name#

property uncertainty#

property values#

class ParameterFitResult#

Bases: pybind11_object

Fitted parameters, model predictions, residuals, and convergence diagnostics.

property evaluations#

property initial_parameters#

property iterations#

property message#

property objective#

property observed#

property parameter_names#

property parameters#

property prediction#

property residuals#

property success#

summary(self: PyMieSim.inverse.ParameterFitResult) → str#

fit_parameters(model: collections.abc.Callable, observation: PyMieSim.inverse.MeasuredData, parameters: collections.abc.Iterable, *, max_iterations: SupportsInt | SupportsIndex = 200, initial_step: SupportsFloat | SupportsIndex = 0.1, step_tolerance: SupportsFloat | SupportsIndex = 1e-06, show_progress: bool = False) → object#

Fit bounded parameters with deterministic coordinate search.

Set show_progress=True to display the current iteration, objective value, and normalized step size. Progress output is disabled by default.