# Licensed under a 3-clause BSD style license - see LICENSE.rst

"""

This module implements classes (called Fitters) which combine optimization

algorithms (typically from `scipy.optimize`) with statistic functions to perform

fitting. Fitters are implemented as callable classes. In addition to the data

to fit, the ``__call__`` method takes an instance of

`~astropy.modeling.core.FittableModel` as input, and returns a copy of the

model with its parameters determined by the optimizer.

Optimization algorithms, called "optimizers" are implemented in

`~astropy.modeling.optimizers` and statistic functions are in

`~astropy.modeling.statistic`. The goal is to provide an easy to extend

framework and allow users to easily create new fitters by combining statistics

with optimizers.

There are two exceptions to the above scheme.

`~astropy.modeling.fitting.LinearLSQFitter` uses Numpy's `~numpy.linalg.lstsq`

function. `~astropy.modeling.fitting.LevMarLSQFitter` uses

`~scipy.optimize.leastsq` which combines optimization and statistic in one

implementation.

"""

# pylint: disable=invalid-name

import abc

import inspect

import warnings

from functools import reduce, wraps

from importlib.metadata import entry_points

import numpy as np

from astropy.units import Quantity

from astropy.utils.decorators import deprecated_renamed_argument

from astropy.utils.exceptions import AstropyDeprecationWarning, AstropyUserWarning

from ._fitting_parallel import FitInfoArrayContainer, parallel_fit_dask

from .optimizers import DEFAULT_ACC, DEFAULT_EPS, DEFAULT_MAXITER, SLSQP, Simplex

from .spline import (

SplineExactKnotsFitter,

SplineInterpolateFitter,

SplineSmoothingFitter,

SplineSplrepFitter,

)

from .statistic import leastsquare

from .utils import _combine_equivalency_dict, poly_map_domain

__all__ = [

"DogBoxLSQFitter",

"FitInfoArrayContainer",

"Fitter",

"FittingWithOutlierRemoval",

"JointFitter",

"LMLSQFitter",

"LevMarLSQFitter",

"LinearLSQFitter",

"ModelLinearityError",

"ModelsError",

"SLSQPLSQFitter",

"SimplexLSQFitter",

"SplineExactKnotsFitter",

"SplineInterpolateFitter",

"SplineSmoothingFitter",

"SplineSplrepFitter",

"TRFLSQFitter",

"parallel_fit_dask",

]

# Statistic functions implemented in `astropy.modeling.statistic.py

STATISTICS = [leastsquare]

# Optimizers implemented in `astropy.modeling.optimizers.py

OPTIMIZERS = [Simplex, SLSQP]

class NonFiniteValueError(RuntimeError):

"""

Error raised when attempting to a non-finite value.

"""

class Covariance:

"""Class for covariance matrix calculated by fitter."""

def __init__(self, cov_matrix, param_names):

self.cov_matrix = cov_matrix

self.param_names = param_names

def pprint(self, max_lines, round_val):

# Print and label lower triangle of covariance matrix

# Print rows for params up to `max_lines`, round floats to 'round_val'

longest_name = max(len(x) for x in self.param_names)

ret_str = "parameter variances / covariances \n"

fstring = f"{'': <{longest_name}}| {{0}}\n"

for i, row in enumerate(self.cov_matrix):

if i <= max_lines - 1:

param = self.param_names[i]

ret_str += fstring.replace(" " * len(param), param, 1).format(

repr(np.round(row[: i + 1], round_val))[7:-2]

)

else:

ret_str += "..."

return ret_str.rstrip()

def __repr__(self):

return self.pprint(max_lines=10, round_val=3)

def __getitem__(self, params):

# index covariance matrix by parameter names or indices

if len(params) != 2:

raise ValueError("Covariance must be indexed by two values.")

if all(isinstance(item, str) for item in params):

i1, i2 = (

self.param_names.index(params[0]),

self.param_names.index(params[1]),

)

elif all(isinstance(item, int) for item in params):

i1, i2 = params

else:

raise TypeError(

"Covariance can be indexed by two parameter names or integer indices."

)

return self.cov_matrix[i1][i2]

class StandardDeviations:

"""Class for fitting uncertainties."""

def __init__(self, cov_matrix, param_names):

self.param_names = param_names

self.stds = self._calc_stds(cov_matrix)

def _calc_stds(self, cov_matrix):

# sometimes scipy lstsq returns a non-sensical negative vals in the

# diagonals of the cov_x it computes.

stds = [np.sqrt(x) if x > 0 else None for x in np.diag(cov_matrix)]

return stds

def pprint(self, max_lines, round_val):

longest_name = max(len(x) for x in self.param_names)

ret_str = "standard deviations\n"

for i, std in enumerate(self.stds):

if i <= max_lines - 1:

param = self.param_names[i]

ret_str += (

f"{param}{' ' * (longest_name - len(param))}| "

f"{np.round(std, round_val)}\n"

)

else:

ret_str += "..."

return ret_str.rstrip()

def __repr__(self):

return self.pprint(max_lines=10, round_val=3)

def __getitem__(self, param):

if isinstance(param, str):

i = self.param_names.index(param)

elif isinstance(param, int):

i = param

else:

raise TypeError(

"Standard deviation can be indexed by parameter name or integer."

)

return self.stds[i]

class ModelsError(Exception):

"""Base class for model exceptions."""

class ModelLinearityError(ModelsError):

"""Raised when a non-linear model is passed to a linear fitter."""

class UnsupportedConstraintError(ModelsError, ValueError):

"""

Raised when a fitter does not support a type of constraint.

"""

def fitter_unit_support(func):

"""

This is a decorator that can be used to add support for dealing with

quantities to any __call__ method on a fitter which may not support

quantities itself. This is done by temporarily removing units from all

parameters then adding them back once the fitting has completed.

"""

@wraps(func)

def wrapper(self, model, x, y, z=None, **kwargs):

equivalencies = kwargs.pop("equivalencies", None)

data_has_units = (

isinstance(x, Quantity)

or isinstance(y, Quantity)

or isinstance(z, Quantity)

)

model_has_units = model._has_units

if data_has_units or model_has_units:

if model._supports_unit_fitting:

# We now combine any instance-level input equivalencies with user

# specified ones at call-time.

input_units_equivalencies = _combine_equivalency_dict(

model.inputs, equivalencies, model.input_units_equivalencies

)

# If input_units is defined, we transform the input data into those

# expected by the model. We hard-code the input names 'x', and 'y'

# here since FittableModel instances have input names ('x',) or

# ('x', 'y')

if model.input_units is not None:

if isinstance(x, Quantity):

x = x.to(

model.input_units[model.inputs[0]],

equivalencies=input_units_equivalencies[model.inputs[0]],

)

if isinstance(y, Quantity) and z is not None:

y = y.to(

model.input_units[model.inputs[1]],

equivalencies=input_units_equivalencies[model.inputs[1]],

)

# Create a dictionary mapping the real model inputs and outputs

# names to the data. This remapping of names must be done here, after

# the input data is converted to the correct units.

rename_data = {model.inputs[0]: x}

if z is not None:

rename_data[model.outputs[0]] = z

rename_data[model.inputs[1]] = y

else:

rename_data[model.outputs[0]] = y

rename_data["z"] = None

# We now strip away the units from the parameters, taking care to

# first convert any parameters to the units that correspond to the

# input units (to make sure that initial guesses on the parameters)

# are in the right unit system

model = model.without_units_for_data(**rename_data)

if isinstance(model, tuple):

rename_data["_left_kwargs"] = model[1]

rename_data["_right_kwargs"] = model[2]

model = model[0]

# We strip away the units from the input itself

add_back_units = False

if isinstance(x, Quantity):

add_back_units = True

xdata = x.value

else:

xdata = np.asarray(x)

if isinstance(y, Quantity):

add_back_units = True

ydata = y.value

else:

ydata = np.asarray(y)

if z is not None:

if isinstance(z, Quantity):

add_back_units = True

zdata = z.value

else:

zdata = np.asarray(z)

# We run the fitting

if z is None:

model_new = func(self, model, xdata, ydata, **kwargs)

else:

model_new = func(self, model, xdata, ydata, zdata, **kwargs)

# And finally we add back units to the parameters

if add_back_units:

model_new = model_new.with_units_from_data(**rename_data)

return model_new

else:

raise NotImplementedError(

"This model does not support being fit to data with units."

)

else:

return func(self, model, x, y, z=z, **kwargs)

return wrapper

class Fitter:

"""

Base class for all fitters.

Parameters

----------

optimizer : callable

A callable implementing an optimization algorithm

statistic : callable

Statistic function

"""

_subclass_registry = set()

def __init_subclass__(cls) -> None:

if not (inspect.isabstract(cls) or cls.__name__.startswith("_")):

Fitter._subclass_registry.add(cls)

supported_constraints = []

def __init__(self, optimizer, statistic):

if optimizer is None:

raise ValueError("Expected an optimizer.")

if statistic is None:

raise ValueError("Expected a statistic function.")

if isinstance(optimizer, type):

# a callable class

self._opt_method = optimizer()

elif inspect.isfunction(optimizer):

self._opt_method = optimizer

else:

raise ValueError("Expected optimizer to be a callable class or a function.")

if isinstance(statistic, type):

self._stat_method = statistic()

else:

self._stat_method = statistic

def objective_function(self, fps, *args):

"""

Function to minimize.

Parameters

----------

fps : list

parameters returned by the fitter

args : list

[model, [other_args], [input coordinates]]

other_args may include weights or any other quantities specific for

a statistic

Notes

-----

The list of arguments (args) is set in the `__call__` method.

Fitters may overwrite this method, e.g. when statistic functions

require other arguments.

"""

model = args[0]

meas = args[-1]

fitter_to_model_params(model, fps)

res = self._stat_method(meas, model, *args[1:-1])

return res

@staticmethod

def _add_fitting_uncertainties(*args):

"""

When available, calculate and sets the parameter covariance matrix

(model.cov_matrix) and standard deviations (model.stds).

"""

return None

@abc.abstractmethod

def __call__(self):

"""

This method performs the actual fitting and modifies the parameter list

of a model.

Fitter subclasses should implement this method.

"""

raise NotImplementedError("Subclasses should implement this method.")

class LinearLSQFitter(Fitter):

"""

A class performing a linear least square fitting.

Uses `numpy.linalg.lstsq` to do the fitting.

Given a model and data, fits the model to the data and changes the

model's parameters. Keeps a dictionary of auxiliary fitting information.

Notes

-----

Note that currently LinearLSQFitter does not support compound models.

"""

supported_constraints = ["fixed"]

supports_masked_input = True

def __init__(self, calc_uncertainties=False):

self.fit_info = {

"residuals": None,

"rank": None,

"singular_values": None,

"params": None,

}

self._calc_uncertainties = calc_uncertainties

@staticmethod

def _is_invertible(m):

"""Check if inverse of matrix can be obtained."""

if m.shape[0] != m.shape[1]:

return False

if np.linalg.matrix_rank(m) < m.shape[0]:

return False

return True

def _add_fitting_uncertainties(self, model, a, n_coeff, x, y, z=None, resids=None): # ruff: ignore[PLR0917]

"""

Calculate and parameter covariance matrix and standard deviations

and set `cov_matrix` and `stds` attributes.

"""

x_dot_x_prime = np.dot(a.T, a)

masked = False or hasattr(y, "mask")

# check if invertible. if not, can't calc covariance.

if not self._is_invertible(x_dot_x_prime):

return model

inv_x_dot_x_prime = np.linalg.inv(x_dot_x_prime)

if z is None: # 1D models

if len(model) == 1: # single model

mask = None

if masked:

mask = y.mask

xx = np.ma.array(x, mask=mask)

RSS = [(1 / (xx.count() - n_coeff)) * resids]

if len(model) > 1: # model sets

RSS = [] # collect sum residuals squared for each model in set

for j in range(len(model)):

mask = None

if masked:

mask = y.mask[..., j].ravel()

xx = np.ma.array(x, mask=mask)

eval_y = model(xx, model_set_axis=False)

eval_y = np.rollaxis(eval_y, model.model_set_axis)[j]

RSS.append(

(1 / (xx.count() - n_coeff)) * np.sum((y[..., j] - eval_y) ** 2)

)

else: # 2D model

if len(model) == 1:

mask = None

if masked:

warnings.warn(

"Calculation of fitting uncertainties "

"for 2D models with masked values not "

"currently supported.\n",

AstropyUserWarning,

)

return

xx, _ = np.ma.array(x, mask=mask), np.ma.array(y, mask=mask)

# len(xx) instead of xx.count. this will break if values are masked?

RSS = [(1 / (len(xx) - n_coeff)) * resids]

else:

RSS = []

for j in range(len(model)):

eval_z = model(x, y, model_set_axis=False)

mask = None # need to figure out how to deal w/ masking here.

if model.model_set_axis == 1:

# model_set_axis passed when evaluating only refers to input shapes

# so output must be reshaped for model_set_axis=1.

eval_z = np.rollaxis(eval_z, 1)

eval_z = eval_z[j]

RSS.append(

[(1 / (len(x) - n_coeff)) * np.sum((z[j] - eval_z) ** 2)]

)

covs = [inv_x_dot_x_prime * r for r in RSS]

free_param_names = [

x

for x in model.fixed

if (model.fixed[x] is False) and (model.tied[x] is False)

]

if len(covs) == 1:

model.cov_matrix = Covariance(covs[0], model.param_names)

model.stds = StandardDeviations(covs[0], free_param_names)

else:

model.cov_matrix = [Covariance(cov, model.param_names) for cov in covs]

model.stds = [StandardDeviations(cov, free_param_names) for cov in covs]

@staticmethod

def _deriv_with_constraints(model, param_indices, x=None, y=None):

if y is None:

d = np.array(model.fit_deriv(x, *model.parameters))

else:

d = np.array(model.fit_deriv(x, y, *model.parameters))

if model.col_fit_deriv:

return d[param_indices]

else:

return d[..., param_indices]

def _map_domain_window(self, model, x, y=None):

"""

Maps domain into window for a polynomial model which has these

attributes.

"""

if y is None:

if hasattr(model, "domain") and model.domain is None:

model.domain = [x.min(), x.max()]

if hasattr(model, "window") and model.window is None:

model.window = [-1, 1]

return poly_map_domain(x, model.domain, model.window)

else:

if hasattr(model, "x_domain") and model.x_domain is None:

model.x_domain = [x.min(), x.max()]

if hasattr(model, "y_domain") and model.y_domain is None:

model.y_domain = [y.min(), y.max()]

if hasattr(model, "x_window") and model.x_window is None:

model.x_window = [-1.0, 1.0]

if hasattr(model, "y_window") and model.y_window is None:

model.y_window = [-1.0, 1.0]

xnew = poly_map_domain(x, model.x_domain, model.x_window)

ynew = poly_map_domain(y, model.y_domain, model.y_window)

return xnew, ynew

@fitter_unit_support

def __call__( # ruff: ignore[PLR0917]

self,

model,

x,

y,

z=None,

weights=None,

rcond=None,

*,

inplace=False,

):

"""

Fit data to this model.

Parameters

----------

model : `~astropy.modeling.FittableModel`

model to fit to x, y, z

x : array

Input coordinates

y : array-like

Input coordinates

z : array-like, optional

Input coordinates.

If the dependent (``y`` or ``z``) coordinate values are provided

as a `numpy.ma.MaskedArray`, any masked points are ignored when

fitting. Note that model set fitting is significantly slower when

there are masked points (not just an empty mask), as the matrix

equation has to be solved for each model separately when their

coordinate grids differ.

weights : array, optional

Weights for fitting.

For data with Gaussian uncertainties, the weights should be

1/sigma.

rcond : float, optional

Cut-off ratio for small singular values of ``a``.

Singular values are set to zero if they are smaller than ``rcond``

times the largest singular value of ``a``.

equivalencies : list or None, optional, keyword-only

List of *additional* equivalencies that are should be applied in

case x, y and/or z have units. Default is None.

inplace : bool, optional

If `False` (the default), a copy of the model with the fitted

parameters set will be returned. If `True`, the returned model will

be the same instance as the model passed in, and the parameter

values will be changed inplace.

Returns

-------

fitted_model : `~astropy.modeling.FittableModel`

If ``inplace`` is `False` (the default), this is a copy of the

input model with parameters set by the fitter. If ``inplace`` is

`True`, this is the same model as the input model, with parameters

updated to be those set by the fitter.

"""

if not model.fittable:

raise ValueError("Model must be a subclass of FittableModel")

if not model.linear:

raise ModelLinearityError(

"Model is not linear in parameters, "

"linear fit methods should not be used."

)

if hasattr(model, "submodel_names"):

raise ValueError("Model must be simple, not compound")

_validate_constraints(self.supported_constraints, model)

model_copy = model if inplace else model.copy()

model_copy.sync_constraints = False

_, fit_param_indices, _ = model_to_fit_params(model_copy)

if model_copy.n_inputs == 2 and z is None:

raise ValueError("Expected x, y and z for a 2 dimensional model.")

farg = _convert_input(

x,

y,

z,

n_models=len(model_copy),

model_set_axis=model_copy.model_set_axis,

verify_dims=_verify_dims_in_fitting(model_copy),

)

n_fixed = sum(model_copy.fixed.values())

# This is also done by _convert_inputs, but we need it here to allow

# checking the array dimensionality before that gets called:

if weights is not None:

weights = np.asarray(weights, dtype=float)

if n_fixed:

# The list of fixed params is the complement of those being fitted:

fixparam_indices = [

idx

for idx in range(len(model_copy.param_names))

if idx not in fit_param_indices

]

# Construct matrix of user-fixed parameters that can be dotted with

# the corresponding fit_deriv() terms, to evaluate corrections to

# the dependent variable in order to fit only the remaining terms:

fixparams = np.asarray(

[

getattr(model_copy, model_copy.param_names[idx]).value

for idx in fixparam_indices

]

)

if len(farg) == 2:

x, y = farg

if weights is not None:

# If we have separate weights for each model, apply the same

# conversion as for the data, otherwise check common weights

# as if for a single model:

_, weights = _convert_input(

x,

weights,

n_models=len(model_copy) if weights.ndim == y.ndim else 1,

model_set_axis=model_copy.model_set_axis,

)

# map domain into window

if hasattr(model_copy, "domain"):

x = self._map_domain_window(model_copy, x)

if n_fixed:

lhs = np.asarray(

self._deriv_with_constraints(model_copy, fit_param_indices, x=x)

)

fixderivs = self._deriv_with_constraints(

model_copy, fixparam_indices, x=x

)

else:

lhs = np.asarray(model_copy.fit_deriv(x, *model_copy.parameters))

sum_of_implicit_terms = model_copy.sum_of_implicit_terms(x)

rhs = y

else:

x, y, z = farg

if weights is not None:

# If we have separate weights for each model, apply the same

# conversion as for the data, otherwise check common weights

# as if for a single model:

_, _, weights = _convert_input(

x,

y,

weights,

n_models=len(model_copy) if weights.ndim == z.ndim else 1,

model_set_axis=model_copy.model_set_axis,

)

# map domain into window

if hasattr(model_copy, "x_domain"):

x, y = self._map_domain_window(model_copy, x, y)

if n_fixed:

lhs = np.asarray(

self._deriv_with_constraints(

model_copy, fit_param_indices, x=x, y=y

)

)

fixderivs = self._deriv_with_constraints(

model_copy, fixparam_indices, x=x, y=y

)

else:

lhs = np.asanyarray(model_copy.fit_deriv(x, y, *model_copy.parameters))

sum_of_implicit_terms = model_copy.sum_of_implicit_terms(x, y)

if len(model_copy) > 1:

# Just to be explicit (rather than baking in False == 0):

model_axis = model_copy.model_set_axis or 0

if z.ndim > 2:

# For higher-dimensional z, flatten all the axes except the

# dimension along which models are stacked and transpose so

# the model axis is *last* (I think this resolves Erik's

# pending generalization from 80a6f25a):

rhs = np.rollaxis(z, model_axis, z.ndim)

rhs = rhs.reshape(-1, rhs.shape[-1])

else:

# This "else" seems to handle the corner case where the

# user has already flattened x/y before attempting a 2D fit

# but z has a second axis for the model set. NB. This is

# ~5-10x faster than using rollaxis.

rhs = z.T if model_axis == 0 else z

if weights is not None:

# Same for weights

if weights.ndim > 2:

# Separate 2D weights for each model:

weights = np.rollaxis(weights, model_axis, weights.ndim)

weights = weights.reshape(-1, weights.shape[-1])

elif weights.ndim == z.ndim:

# Separate, flattened weights for each model:

weights = weights.T if model_axis == 0 else weights

else:

# Common weights for all the models:

weights = weights.ravel()

else:

rhs = z.ravel()

if weights is not None:

weights = weights.ravel()

# If the derivative is defined along rows (as with non-linear models)

if model_copy.col_fit_deriv:

lhs = np.asarray(lhs).T

# Some models (eg. Polynomial1D) don't flatten multi-dimensional inputs

# when constructing their Vandermonde matrix, which can lead to obscure

# failures below. Ultimately, np.linalg.lstsq can't handle >2D matrices,

# so just raise a slightly more informative error when this happens:

if np.asanyarray(lhs).ndim > 2:

raise ValueError(

f"{type(model_copy).__name__} gives unsupported >2D "

"derivative matrix for this x/y"

)

# Subtract any terms fixed by the user from (a copy of) the RHS, in

# order to fit the remaining terms correctly:

if n_fixed:

if model_copy.col_fit_deriv:

fixderivs = np.asarray(fixderivs).T # as for lhs above

rhs = rhs - fixderivs.dot(fixparams) # evaluate user-fixed terms

# Subtract any terms implicit in the model from the RHS, which, like

# user-fixed terms, affect the dependent variable but are not fitted:

if sum_of_implicit_terms is not None:

# If we have a model set, the extra axis must be added to

# sum_of_implicit_terms as its innermost dimension, to match the

# dimensionality of rhs after _convert_input "rolls" it as needed

# by np.linalg.lstsq. The vector then gets broadcast to the right

# number of sets (columns). This assumes all the models share the

# same input coordinates, as is currently the case.

if len(model_copy) > 1:

sum_of_implicit_terms = sum_of_implicit_terms[..., np.newaxis]

rhs = rhs - sum_of_implicit_terms

if weights is not None:

if rhs.ndim == 2:

if weights.shape == rhs.shape:

# separate weights for multiple models case: broadcast

# lhs to have more dimension (for each model)

lhs = lhs[..., np.newaxis] * weights[:, np.newaxis]

rhs = rhs * weights

else:

lhs *= weights[:, np.newaxis]

# Don't modify in-place in case rhs was the original

# dependent variable array

rhs = rhs * weights[:, np.newaxis]

else:

lhs *= weights[:, np.newaxis]

rhs = rhs * weights

scl = (lhs * lhs).sum(0)

lhs /= scl

masked = np.any(np.ma.getmask(rhs))

if weights is not None and not masked and np.any(np.isnan(lhs)):

raise ValueError(

"Found NaNs in the coefficient matrix, which "

"should not happen and would crash the lapack "

"routine. Maybe check that weights are not null."

)

a = None # need for calculating covariance

if (masked and len(model_copy) > 1) or (

weights is not None and weights.ndim > 1

):

# Separate masks or weights for multiple models case: Numpy's

# lstsq supports multiple dimensions only for rhs, so we need to

# loop manually on the models. This may be fixed in the future

# with https://github.com/numpy/numpy/pull/15777.

# Initialize empty array of coefficients and populate it one model

# at a time. The shape matches the number of coefficients from the

# Vandermonde matrix and the number of models from the RHS:

lacoef = np.zeros(lhs.shape[1:2] + rhs.shape[-1:], dtype=rhs.dtype)

# Arrange the lhs as a stack of 2D matrices that we can iterate

# over to get the correctly-orientated lhs for each model:

if lhs.ndim > 2:

lhs_stack = np.rollaxis(lhs, -1, 0)

else:

lhs_stack = np.broadcast_to(lhs, rhs.shape[-1:] + lhs.shape)

# Loop over the models and solve for each one. By this point, the

# model set axis is the second of two. Transpose rather than using,

# say, np.moveaxis(array, -1, 0), since it's slightly faster and

# lstsq can't handle >2D arrays anyway. This could perhaps be

# optimized by collecting together models with identical masks

# (eg. those with no rejected points) into one operation, though it

# will still be relatively slow when calling lstsq repeatedly.

for model_lhs, model_rhs, model_lacoef in zip(lhs_stack, rhs.T, lacoef.T):

# Cull masked points on both sides of the matrix equation:

good = ~model_rhs.mask if masked else slice(None)

model_lhs = model_lhs[good]

model_rhs = model_rhs[good][..., np.newaxis]

a = model_lhs

# Solve for this model:

t_coef, resids, rank, sval = np.linalg.lstsq(

model_lhs, model_rhs, rcond

)

model_lacoef[:] = t_coef.T

else:

# If we're fitting one or more models over a common set of points,

# we only have to solve a single matrix equation, which is an order

# of magnitude faster than calling lstsq() once per model below:

good = ~rhs.mask if masked else slice(None) # latter is a no-op

a = lhs[good]

# Solve for one or more models:

lacoef, resids, rank, sval = np.linalg.lstsq(lhs[good], rhs[good], rcond)

self.fit_info["residuals"] = resids

self.fit_info["rank"] = rank

self.fit_info["singular_values"] = sval

lacoef /= scl[:, np.newaxis] if scl.ndim < rhs.ndim else scl

self.fit_info["params"] = lacoef

fitter_to_model_params(model_copy, lacoef.ravel())

# TODO: Only Polynomial models currently have an _order attribute;

# maybe change this to read isinstance(model, PolynomialBase)

if (

hasattr(model_copy, "_order")

and len(model_copy) == 1

and rank < (model_copy._order - n_fixed)

):

warnings.warn("The fit may be poorly conditioned\n", AstropyUserWarning)

# calculate and set covariance matrix and standard devs. on model

if self._calc_uncertainties:

if len(y) > len(lacoef):

self._add_fitting_uncertainties(

model_copy, a * scl, len(lacoef), x, y, z, resids

)

model_copy.sync_constraints = True

return model_copy

class FittingWithOutlierRemoval:

"""

This class combines an outlier removal technique with a fitting procedure.

Basically, given a maximum number of iterations ``niter``, outliers are

removed and fitting is performed for each iteration, until no new outliers

are found or ``niter`` is reached.

Parameters

----------

fitter : `Fitter`

An instance of any Astropy fitter, i.e., LinearLSQFitter,

LevMarLSQFitter, SLSQPLSQFitter, SimplexLSQFitter, JointFitter. For

model set fitting, this must understand masked input data (as

indicated by the fitter class attribute ``supports_masked_input``).

outlier_func : callable

A function for outlier removal.

If this accepts an ``axis`` parameter like the `numpy` functions, the

appropriate value will be supplied automatically when fitting model

sets (unless overridden in ``outlier_kwargs``), to find outliers for

each model separately; otherwise, the same filtering must be performed

in a loop over models, which is almost an order of magnitude slower.

niter : int, optional

Maximum number of iterations.

outlier_kwargs : dict, optional

Keyword arguments for outlier_func.

Attributes

----------

fit_info : dict

The ``fit_info`` (if any) from the last iteration of the wrapped

``fitter`` during the most recent fit. An entry is also added with the

keyword ``niter`` that records the actual number of fitting iterations

performed (as opposed to the user-specified maximum).

"""

def __init__(self, fitter, outlier_func, niter=3, **outlier_kwargs):

self.fitter = fitter

self.outlier_func = outlier_func

self.niter = niter

self.outlier_kwargs = outlier_kwargs

self.fit_info = {"niter": None}

def __str__(self):

return (

f"Fitter: {self.fitter.__class__.__name__}\n"

f"Outlier function: {self.outlier_func.__name__}\n"

f"Num. of iterations: {self.niter}\n"

f"Outlier func. args.: {self.outlier_kwargs}"

)

def __repr__(self):

return (

f"{self.__class__.__name__}(fitter: {self.fitter.__class__.__name__}, "

f"outlier_func: {self.outlier_func.__name__},"

f" niter: {self.niter}, outlier_kwargs: {self.outlier_kwargs})"

)

def __call__(self, model, x, y, z=None, weights=None, *, inplace=False, **kwargs):

"""

Parameters

----------

model : `~astropy.modeling.FittableModel`

An analytic model which will be fit to the provided data.

This also contains the initial guess for an optimization

algorithm.

x : array-like

Input coordinates.

y : array-like

Data measurements (1D case) or input coordinates (2D case).

z : array-like, optional

Data measurements (2D case).

weights : array-like, optional

Weights to be passed to the fitter.

kwargs : dict, optional

Keyword arguments to be passed to the fitter.

inplace : bool, optional

If `False` (the default), a copy of the model with the fitted

parameters set will be returned. If `True`, the returned model will

be the same instance as the model passed in, and the parameter

values will be changed inplace.

Returns

-------

fitted_model : `~astropy.modeling.FittableModel`

If ``inplace`` is `False` (the default), this is a copy of the

input model with parameters set by the fitter. If ``inplace`` is

`True`, this is the same model as the input model, with parameters

updated to be those set by the fitter.

mask : `numpy.ndarray`

Boolean mask array, identifying which points were used in the final

fitting iteration (False) and which were found to be outliers or

were masked in the input (True).

"""

# For single models, the data get filtered here at each iteration and

# then passed to the fitter, which is the historical behavior and

# works even for fitters that don't understand masked arrays. For model

# sets, the fitter must be able to filter masked data internally,

# because fitters require a single set of x/y coordinates whereas the

# eliminated points can vary between models. To avoid this limitation,

# we could fall back to looping over individual model fits, but it

# would likely be fiddly and involve even more overhead (and the

# non-linear fitters don't work with model sets anyway, as of writing).

if len(model) == 1:

model_set_axis = None

else:

if (

not hasattr(self.fitter, "supports_masked_input")

or self.fitter.supports_masked_input is not True

):

raise ValueError(

f"{type(self.fitter).__name__} cannot fit model sets with masked "

"values"

)

# Fitters use their input model's model_set_axis to determine how

# their input data are stacked:

model_set_axis = model.model_set_axis

# Construct input coordinate tuples for fitters & models that are

# appropriate for the dimensionality being fitted:

if z is None:

coords = (x,)

data = y

else:

coords = x, y

data = z

# For model sets, construct a numpy-standard "axis" tuple for the

# outlier function, to treat each model separately (if supported):

if model_set_axis is not None:

if model_set_axis < 0:

model_set_axis += data.ndim