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[Any, ...], 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[Any, ...], dtype[float]]]

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

The degree of the monomials.

derive_props(pars, covar)[source]

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

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[float]]]

Parameters:
get_par_nonlinear()[source]

Return a boolean mask for the nonlinear parameters.

Return type:

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

property name
property npar

Return the number of parameters.

property par_scales: ndarray[tuple[Any, ...], 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[Any, ...], dtype[float]]

Parameters:
class LorentzModel(relerr_threshold=0.1)[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.

Parameters:

relerr_threshold (float)

denom_degrees: ndarray[tuple[Any, ...], 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(pars, covar)[source]

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

The exponential correlation time is derived from the parameters, if possible. If not, or if the variance of the estimate is too large, the “criterion” is set to infinity and the “msg” is set accordingly, to discard the current fit from the average over the cutoff frequencies.

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[float]]]

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

The degrees of the monomials in the numerator.

relerr_threshold: float

A threshold for the relative error on the exponential correlation time.

The relative error is defined as ratio of the standard deviation to the mean of the correlation time.

To accept the current cutoff frequency, the relative error must be below this threshold. This eliminates high-variance cases for which the maximum a posteriori estimate tends to be a poor approximation.

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[Any, ...], dtype[float]]]

Parameters:
denom_degrees: ndarray[tuple[Any, ...], 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(pars, covar)[source]

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

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[float]]]

Parameters:
get_par_nonlinear()[source]

Return a boolean mask for the nonlinear parameters.

Return type:

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

property name
property npar

Return the number of parameters.

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

The degrees of the monomials in the numerator.

property par_scales: ndarray[tuple[Any, ...], 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[Any, ...], 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[Any, ...], dtype[float]]) – The frequencies for which the model spectrum amplitudes are computed.

  • pars (ndarray[tuple[Any, ...], 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[Any, ...], 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[Any, ...], dtype[float]]

Parameters:
derive_props(pars, covar)[source]

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

Parameters:
Return type:

dict[str, ndarray[tuple[Any, ...], dtype[float]]]

Returns:

props – A dictionary with additional properties, whose calculation requires model-specific knowledge.

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[Any, ...], 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[Any, ...], dtype[float]]]

Parameters:
property npar

Return the number of parameters.

property par_scales: ndarray[tuple[Any, ...], 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[Any, ...], 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[Any, ...], 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[Any, ...], dtype[float]])

which_invalid(pars)[source]

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

Return type:

ndarray[tuple[Any, ...], 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[Any, ...], dtype[float]]) – The frequencies for which the model spectrum amplitudes are computed.

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

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

  • weights (ndarray[tuple[Any, ...], 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.