stacie.model module

Models to fit the low-frequency part of the spectrum.

class ExpPolyModel(degrees)[source]

Bases: SpectrumModel

Exponential function of a linear combination of simple monomials.

Parameters:

degrees (ndarray[tuple[int, ...], dtype[int]])

bounds()[source]

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)[source]

See SpectrumModel.compute().

Return type:

list[ndarray[tuple[int, ...], dtype[float]]]

Parameters:
degrees: ndarray[tuple[int, ...], dtype[int]]

The degree of the monomials.

derive_props(props)[source]

Add the autocorrelation integral (and other properties) derived from the parameters.

Parameters:

props (dict[str, ndarray[tuple[int, ...], dtype[float]]])

get_par_nonlinear()[source]

Return a boolean mask for the nonlinear parameters.

Return type:

ndarray[tuple[int, ...], dtype[bool]]

property name
property npar

Return the number of parameters.

property par_scales: ndarray[tuple[int, ...], dtype[float]]

Return the scales of the parameters and the cost function.

solve_linear(freqs, ndofs, amplitudes, weights, nonlinear_pars)[source]

Use linear linear regression to solve a subset of the parameters.

This is a specialized implementation that rewrites the problem in a different form to solve all parameters with a linear regression.

Return type:

ndarray[tuple[int, ...], dtype[float]]

Parameters:
class LorentzModel(*, ratio_weight=1.0, ratio_threshold=100.0)[source]

Bases: PadeModel

A model for the spectrum with a Lorentzian peak at zero frequency plus some white noise.

This is a special case of the PadeModel with numer_degrees = [0, 2] and denom_degrees = [2]. Furthermore, it will only accept parameters that correspond to a well-defined exponential correlation time.

For too small cutoffs (covering only the peak of the Lorentzian and not its decay), the estimates of the Lorentzian width, and consequently the exponential correlation time, become statistically unreliable. In this regime, the assumption of maximum a posteriori probability (MAP), on which STACIE relies to fit the model and estimate parameter uncertainties, also breaks down. Unreliable MAP estimates are inferred from the relative error of the exponential correlation time divided by the relative error of the autocorrelation integral. This implementation uses the ratio in two ways:

  1. When the ratio exceeds a predefined threshold (default value 100), the cutoff criterion is set to infinity.

  2. If the ratio remains below this threshold, the ratio times a weight (default value 1.0) is added to the cutoff criterion.

Note that this is an empirical penalty to mitigate MAP-related issues. Because the penalty is expressed as a ratio of relative errors, it is dimensionless and insensitive to the overall uncertainty of the spectrum.

The hyperparameters ratio_weight and ratio_threshold may be tuned to adjust the sensitivity of the heuristic, but it is recommended to keep their default values.

Parameters:
denom_degrees: ndarray[tuple[int, ...], dtype[int]]

The degrees of the monomials in the denominator.

Note that the leading term is always 1, and there is no need to include degree zero.

derive_props(props)[source]

Add the autocorrelation integral (and other properties) derived from the parameters.

The exponential correlation time is derived from the parameters, if the fitted model has a maximum at zero frequency. If not, the “criterion” is set to infinity and the “msg” is set accordingly, to discard the current fit from the average over the cutoff frequencies.

Parameters:

props (dict[str, ndarray[tuple[int, ...], dtype[float]]])

property name
numer_degrees: ndarray[tuple[int, ...], dtype[int]]

The degrees of the monomials in the numerator.

ratio_threshold: float

A threshold for the ratio of relative errors used to set the cutoff criterion to Inf.

ratio_weight: float

The penalty for the cutoff criterion is this weight times the ratio of relative errors.

class PadeModel(numer_degrees, denom_degrees)[source]

Bases: SpectrumModel

A rational function model for the spectrum, a.k.a. a Padé approximation.

Parameters:
bounds()[source]

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)[source]

See SpectrumModel.compute().

Return type:

list[ndarray[tuple[int, ...], dtype[float]]]

Parameters:
denom_degrees: ndarray[tuple[int, ...], dtype[int]]

The degrees of the monomials in the denominator.

Note that the leading term is always 1, and there is no need to include degree zero.

derive_props(props)[source]

Add the autocorrelation integral (and other properties) derived from the parameters.

Parameters:

props (dict[str, ndarray[tuple[int, ...], dtype[float]]])

get_par_nonlinear()[source]

Return a boolean mask for the nonlinear parameters.

Return type:

ndarray[tuple[int, ...], dtype[bool]]

property name
property npar

Return the number of parameters.

numer_degrees: ndarray[tuple[int, ...], dtype[int]]

The degrees of the monomials in the numerator.

property par_scales: ndarray[tuple[int, ...], dtype[float]]

Return the scales of the parameters and the cost function.

solve_linear(freqs, ndofs, amplitudes, weights, nonlinear_pars)[source]

Use linear linear regression to solve a subset of the parameters.

This is a specialized implementation that rewrites the problem in a different form to solve all parameters with a linear regression.

Return type:

ndarray[tuple[int, ...], dtype[float]]

Parameters:
class SpectrumModel[source]

Bases: object

Abstract base class for spectrum models.

Subclasses must override all methods that raise NotImplementedError.

The first parameter must have a property that is used when constructing an initial guess: When the first parameter increases, the model should increase everywhere, and must allow for an arbitrary increase of the spectrum at all points. This is used to repair initial guesses that result in a partially negative spectrum.

bounds()[source]

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)[source]

Compute the amplitudes of the spectrum model.

Parameters:

  • freqs (ndarray[tuple[int, ...], dtype[float]]) – The frequencies for which the model spectrum amplitudes are computed.

  • pars (ndarray[tuple[int, ...], dtype[float]]) – The parameter vector. For vectorized calculations, the last axis corresponds to the parameter index.

  • deriv (int) – The maximum order of derivatives to compute: 0, 1 or 2.

Return type:

list[ndarray[tuple[int, ...], dtype[float]]]

Returns:

results – A results list, index corresponds to order of derivative. The shape of the arrays in the results list is as follows:

  • For deriv=0, the shape is (*vec_shape, len(freqs)).

  • For deriv=1, the shape is (*vec_shape, len(pars), len(freqs)).

  • For deriv=2, the shape is (*vec_shape, len(pars), len(pars), len(freqs))

If some derivatives are independent of the parameters, broadcasting rules may be used to reduce the memory footprint. This means that vec_shape may be replaced by a tuple of ones with the same length.

configure_scales(timestep, freqs, amplitudes)[source]

Store essential scale information in the scales attribute.

Other methods may access this information, so this method should be called before performing any computations.

Return type:

ndarray[tuple[int, ...], dtype[float]]

Parameters:
derive_props(props)[source]

Add the autocorrelation integral (and other properties) derived from the parameters.

Parameters:

props (dict[str, ndarray[tuple[int, ...], dtype[float]]]) – The properties dictionary, including the parameters and their uncertainties. Subclasses may add additional properties to this dictionary.

get_par_nonlinear()[source]

Return a boolean mask for the nonlinear parameters.

The returned parameters cannot be solved with the solve_linear method. Models are free to decide which parameters can be solved with linear regression. For example, some non-linear parameters may be solved with a linear regression after rewriting the regression problem in a different form.

Return type:

ndarray[tuple[int, ...], dtype[bool]]

property name
neglog_prior(pars, *, deriv=0)[source]

Minus logarithm of the prior probability density function, if any.

Subclasses may implement (a very weak) prior, if any.

Return type:

list[ndarray[tuple[int, ...], dtype[float]]]

Parameters:
property npar

Return the number of parameters.

property par_scales: ndarray[tuple[int, ...], dtype[float]]

Return the scales of the parameters and the cost function.

sample_nonlinear_pars(rng, budget)[source]

Return samples of the nonlinear parameters.

Parameters:

  • rng (Generator) – The random number generator.

  • budget (int) – The number of samples to generate.

  • freqs – The frequencies for which the model spectrum amplitudes are computed.

  • par_scales – The scales of the parameters and the cost function.

Return type:

ndarray[tuple[int, ...], dtype[float]]

Returns:

samples – The samples of the nonlinear parameters, array with shape (budget, num_nonlinear), where num_nonlinear is the number of nonlinear parameters.

scales: dict[str, float]

A dictionary with essential scale information for the parameters and the cost function.

solve_linear(freqs, ndofs, amplitudes, weights, nonlinear_pars)[source]

Use linear linear regression to solve a subset of the parameters.

The default implementation in the base class assumes that the linear parameters are genuinely linear without rewriting the regression problem in a different form.

Parameters:
Return type:

ndarray[tuple[int, ...], dtype[float]]

Returns:

  • linear_pars – The solved linear parameters.

  • amplitudes_model – The model amplitudes computed with the solved parameters.

valid(pars)[source]

Return True when the parameters are within the feasible region.

Return type:

bool

Parameters:

pars (ndarray[tuple[int, ...], dtype[float]])

which_invalid(pars)[source]

Return a boolean mask for the parameters outside the feasible region.

Return type:

ndarray[tuple[int, ...], dtype[bool]]

guess(freqs, ndofs, amplitudes, weights, model, rng, nonlinear_budget)[source]

Guess initial values of the parameters for a model.

Parameters:

  • freqs (ndarray[tuple[int, ...], dtype[float]]) – The frequencies for which the model spectrum amplitudes are computed.

  • ndofs (ndarray[tuple[int, ...], dtype[float]]) – The number of degrees of freedom at each frequency.

  • amplitudes (ndarray[tuple[int, ...], dtype[float]]) – The amplitudes of the spectrum.

  • weights (ndarray[tuple[int, ...], dtype[float]]) – Fitting weights, in range [0, 1], to use for each grid point.

  • model (SpectrumModel) – The model for which the parameters are guessed.

  • rng (Generator) – The random number generator.

  • nonlinear_budget (int) – The number of samples of the nonlinear parameters is computed as nonlinear_budget ** num_nonlinear, where num_nonlinear is the number of nonlinear parameters.

Returns:

pars – An initial guess of the parameters.