[roytsmart]

Adds ndoptimize.least_squares(), which fits a model to many independent datasets at once (every spectrum of a spectral cube, for example) with the Levenberg-Marquardt method, compiled and run in parallel with Numba.

Each fit can have its own independent variable, data, noise, initial guess, bounds and Gaussian priors, all broadcast against each other. NaN data points are left out of their fit. The noise of neighboring points may be correlated, as it is in resampled data. That makes the fit generalized least squares, with a banded Cholesky whitening, and gives correct parameter uncertainties. The result is a LeastSquaresResult with the parameters, their covariance, χ², the number of points used, the iteration count and a converged flag.

Algorithm

• Marquardt damping, A + λ diag(A), with Nielsen's gain-ratio update.
• Projected bounds: parameters at an active bound are held fixed for that step, as are parameters with no effect on the model.
• The Newton decrement gᵀA⁻¹g is the convergence test.
• The noise covariance and its Cholesky factor are stored as a band, so memory grows linearly with the number of points.
• The covariance is the inverse of the normal matrix. A parameter that nothing constrains gets infinite variance. If the rest of the normal matrix is singular, as in a fit with fewer usable points than parameters, the covariance is NaN.
• Invalid inputs raise ValueError: bounds with lower > upper, a non-positive prior_sigma, a non-finite prior_mean with a prior, a correlation that isn't finite and positive definite, or a guess that doesn't match the model.

Models

Models are Numba functions that write the model and its Jacobian into output arrays. The new decorator compiles them to one explicit signature and names their parameters:

@ndoptimize.model(parameters=("intercept", "slope"))
def line(x, parameters, f, jacobian):
    ...

It returns an ndoptimize.Model. least_squares() checks the guess against the model's parameter names, since the compiled function doesn't check its indices. Calling a Model directly checks the shapes of its arrays. A model is copied and pickled by reference, like a function, so a copy of an object that holds one keeps the same model. The kernel takes the model as a typed function pointer, so one cached kernel serves every model, including user-defined ones. A new process imports and runs its first fit in about 1 s once the cache is built. cache=False is for models defined outside a file, or ones that read globals computed at run time.

ndoptimize.models provides gaussian, gaussian_linear and two_gaussians_linear.

Evaluation

ndoptimize.evaluate(model, x, parameters) evaluates a model for many sets of parameters at once, in parallel. It broadcasts the batch axes of the points and the parameters like least_squares() does, so the best-fit model of every fit is one call. The named-arrays wrapper, na.ModelFitFunctionArray (sun-data/named-arrays#271), uses it to evaluate a fit at new inputs and for its predictions.

Package setup

• BSD-3-Clause license, with PEP 639 metadata, plus CITATION.cff.
• Tests and typing workflows copied from ndfilters, including the Codecov upload, plus a Codecov badge in the README.
• numpy and numba are the only runtime dependencies (astropy is dropped).
• A new README and documentation introduction. The least_squares docstring has an example that fits 1,000 noisy lines and plots the centroid errors against a unit normal.

The 2024 work-in-progress stubs were never committed and aren't included.

Testing

47 tests:

• Against scipy.optimize.least_squares for every combination of correlated noise, priors and bounds: the same minimum to 10⁻³σ and the same covariance.
• Calibration over 500 noise realizations. With the correlation given, the reported uncertainties match the scatter; without it, they're too small.
• Broadcasting (including arguments broadcast over only some batch axes), missing data, an unconstrained parameter, a degenerate fit, an exact linear fit, a parameter 10²⁷ times the data's scale, a NaN starting point, a poor initial guess, the iteration cap and invalid inputs.
• The models against finite differences, shape checks on direct calls, an uncached model defined with exec, and copying and pickling.
• evaluate against the model evaluated one set at a time, with broadcasting, at the best-fit parameters of noiseless data, and with invalid shapes.

The suite takes 3 s with Numba and 17 s with NUMBA_DISABLE_JIT=1, at 99% coverage. The misses are numerical-safety branches. pyright 1.1.411, ruff and black are clean, and the docs build without warnings.

On the IRIS Fe XII 2×4-binned raster (51,400 spectra), this kernel reproduces the existing numpy GLS fit (99.9% of parameters within 0.05σ) in about 0.1 s, about 600× faster. A synthetic benchmark of the same size (Gaussian plus line, priors, bounds and correlated noise) takes 74 ms on 48 threads.

🤖 Generated with Claude Code

https://claude.ai/code/session_01PZRXiFhGzzrF87vsFg73B9

[codecov]

Welcome to Codecov 🎉

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered ☂️