v2.16.0

This release adds na.ModelFitFunctionArray, which fits a model compiled with ndoptimize to every set of outputs of a function array at once. It also fixes Cartesian2dMatrixArray.power() for matrices with units and for equal or nearly equal eigenvalues.

Fitting a model to a function array

na.ModelFitFunctionArray is a function array whose outputs are fit by a nonlinear least-squares model, with the Levenberg-Marquardt method of ndoptimize.least_squares() (#271). It follows na.PolynomialFitFunctionArray, with a model compiled by ndoptimize.model in place of the polynomial, so one array can fit a Gaussian to every spectrum of a spectral cube. Its abstract base, na.AbstractModelFunctionArray, is the counterpart of na.AbstractPolynomialFunctionArray.

fit = na.ModelFitFunctionArray(
    inputs=velocity,
    outputs=intensity,
    model=ndoptimize.models.gaussian,
    guess=na.CartesianNdVectorArray(dict(
        amplitude=5 * u.DN,
        center=0 * u.km / u.s,
        width=15 * u.km / u.s,
    )),
    variance=np.square(sigma),
    axis_model="velocity",
)
fit.parameters.center          # the best-fit centers, in km / s
fit.covariance.center.width    # in km2 / s2
fit(velocity_fine)             # the best-fit model on a finer grid

• Batches: the points of each fit vary along axis_model, and every other axis is a batch axis, along which each set of outputs gets a fit of its own. A guess which varies along an axis of its own starts several fits from different points.
• Parameters: the guess, the bounds, the priors and the results are vectors whose components are named after ndoptimize.Model.parameters, so a partial bound or prior is a smaller vector.
• Units: the model sees the inputs and outputs in their own units, and each parameter in the unit of its guess, so the guess must be in the units which the model expects for them. The variance, bounds and priors are converted, and the results get their units back, with the product of the units of two parameters for each element of the covariance.
• Vector inputs: the components of vector inputs are the inputs of the model, in order, or the ones components_model selects, like the "wavelength" of a SpectralPositionalVectorArray. More than one points axis is flattened into one, for a 2D Gaussian fit to image cutouts.
• Bins: along an axis where the inputs are the edges of bins, each input becomes the lower and upper edges of each bin, for a model of the average over each bin, like ndoptimize.models.gaussian_binned. The inputs are checked against the ones the model declares.
• Masks and noise: where_model and NaN outputs leave points out of their fit, and correlation describes noise which is correlated between neighboring points.
• Uncertainties: any field may be uncertain, including a single component of a vector. Then the nominal values and the samples are fit separately, and every result is uncertain.
• A function array like any other: every field is carried through indexing, stacking and arithmetic. The fit is computed when one of its results is first used, and again after a field, or elements of the array, are assigned to.

named-arrays now depends on ndoptimize~=0.2.0, which is imported only when a fit is computed.

Behavior change: numpy.copyto on a function array assigns its inputs and outputs again after writing them in place, as assigning to its elements does, so that a subclass which caches a result of its fields, like ModelFitFunctionArray, computes it again.

Powers of 2×2 matrices

Cartesian2dMatrixArray.power() is fixed in two ways:

• Units (#272): it took the square root of the discriminant with np.emath.sqrt, which astropy doesn't implement, so it warned and dropped the units of the root. A matrix with dimensions, such as one in mm, raised a UnitConversionError, and one in a dimensionless unit with a scale, such as mm / m, came out wrong. One in an unscaled dimensionless unit, like the transfer matrix of a periodic multilayer of optika, came out right, but warned on every call. The root is now taken of the value, and the square root of the unit is applied to it.
• Equal and nearly equal eigenvalues (#273): with equal eigenvalues, such as a multiple of the identity, a Jordan block or a periodic multilayer at the edge of its band, every entry came out NaN, and with nearly equal ones, digits were lost. The power is now the polynomial which interpolates it at the eigenvalues, in the form of Newton, which holds where they are equal too. It agrees with numpy.linalg.matrix_power to within 5e-16.

Full Changelog: v2.15.0...v2.16.0