stacie.model module

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

class DebyeModel(*, ratio_weight=1.0, ratio_threshold=100.0)

Bases: PadeModel

A Debye relaxation on top of white noise, with a minimum at zero frequency.

This is a special case of the PadeModel with numer_degrees = [0, 2] and denom_degrees = [2]:

\[C(f) = \frac{p_0 + p_2 f^2}{1 + q_2 f^2} = C_\text{DC} + C_\Delta \frac{(2 \pi f \tau_\text{exp})^2}{1 + (2 \pi f \tau_\text{exp})^2}\]

with \(C_\text{DC} = p_0\), \(C_\Delta = p_2/q_2 - p_0\) and \(\tau_\text{exp} = \sqrt{q_2} / (2 \pi)\). (These parameters are stored in pars_debye in the same order.) The zero-frequency (DC) level \(C_\text{DC}\) is the autocorrelation integral. Only parameters with \(q_2 > 0\), \(C_\Delta > 0\) and \(C_\text{DC} > 0\) are accepted.

For API compatibility with the Lorentz model, the relaxation time is labeled as \(\tau_\text{exp}\), although it is not the decay time of an exponential tail in the autocorrelation function. It is still fair to think of it as the “slowest” timescale.

This model is completely analogous to the Lorentz model, except that it features a minimum at frequency zero. This type of spectrum is common in systems whose degrees of freedom are temporarily confined by their environment: on short timescales they oscillate in place, while on longer timescales they escape and exhibit diffusive (random-walk) behavior. To keep the number of fitting parameters minimal, the model only captures the onset towards the resonance peak of the confined oscillation.

Selection of cutoffs. The selection of cutoffs mirrors that of the LorentzModel (which requires a maximum at zero frequency instead): for too small cutoffs, covering only the bottom of the dip, the relaxation time is poorly determined and the MAP approximation breaks down. When the relative error of \(\tau_\text{exp}\), divided by that of the autocorrelation integral, exceeds ratio_threshold, the cutoff criterion is set to infinity. Otherwise, this ratio times ratio_weight is added to the cutoff criterion.

Parameters:

• ratio_weight (float)

• ratio_threshold (float)

denom_degrees: NDArray[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)

Add the autocorrelation integral and the Debye parameters (if well-defined).

The exponential correlation time is derived from the parameters, if the fitted model has a minimum 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[float]])

property name

numer_degrees: NDArray[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 ExpPolyModel(degrees)

Bases: SpectrumModel

Exponential function of a linear combination of simple monomials.

Parameters:

degrees (NDArray[int])

bounds()

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)

See SpectrumModel.compute().

Return type:

list[NDArray[float]]

Parameters:

• freqs (NDArray[float])

• pars (NDArray[float])

• deriv (int)

degrees: NDArray[int]

The degree of the monomials.

derive_props(props)

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

Parameters:

props (dict[str, NDArray[float]])

get_par_nonlinear()

Return a boolean mask for the nonlinear parameters.

Return type:

NDArray[bool]

property name

property npar

Return the number of parameters.

property par_scales: NDArray[float]

Return the scales of the parameters and the cost function.

solve_linear(freqs, ndofs, amplitudes, weights, nonlinear_pars)

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[float]

Parameters:

• freqs (NDArray[float])

• ndofs (NDArray[float])

• amplitudes (NDArray[float])

• weights (NDArray[float])

• nonlinear_pars (NDArray[float])

class LorentzModel(*, ratio_weight=1.0, ratio_threshold=100.0)

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]:

\[C(f) = \frac{p_0 + p_2 f^2}{1 + q_2 f^2} = C_0 + \frac{2 C_1 \tau_\text{exp}}{1 + (2 \pi f \tau_\text{exp})^2}\]

with \(C_0 = p_2/q_2\), \(C_1 = \pi (p_0 - p_2/q_2) / \sqrt{q_2}\) and \(\tau_\text{exp} = \sqrt{q_2} / (2 \pi)\). (These parameters are stored in pars_lorentz in the same order.) Furthermore, it will only accept parameters that correspond to a well-defined exponential correlation time and a positive autocorrelation integral, i.e. \(q_2 > 0\), \(C_1 > 0\) and \(p_0 > 0\).

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:

• ratio_weight (float)

• ratio_threshold (float)

denom_degrees: NDArray[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)

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[float]])

property name

numer_degrees: NDArray[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)

Bases: SpectrumModel

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

Parameters:

• numer_degrees (NDArray[int])

• denom_degrees (NDArray[int])

bounds()

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)

See SpectrumModel.compute().

Return type:

list[NDArray[float]]

Parameters:

• freqs (NDArray[float])

• pars (NDArray[float])

• deriv (int)

denom_degrees: NDArray[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)

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

Parameters:

props (dict[str, NDArray[float]])

get_par_nonlinear()

Return a boolean mask for the nonlinear parameters.

Return type:

NDArray[bool]

property name

property npar

Return the number of parameters.

numer_degrees: NDArray[int]

The degrees of the monomials in the numerator.

property par_scales: NDArray[float]

Return the scales of the parameters and the cost function.

solve_linear(freqs, ndofs, amplitudes, weights, nonlinear_pars)

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[float]

Parameters:

• freqs (NDArray[float])

• ndofs (NDArray[float])

• amplitudes (NDArray[float])

• weights (NDArray[float])

• nonlinear_pars (NDArray[float])

class SpectrumModel

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()

Return parameter bounds for the optimizer.

Return type:

list[tuple[float, float]]

compute(freqs, pars, *, deriv=0)

Compute the amplitudes of the spectrum model.

Parameters:

• freqs (NDArray[float]) – The frequencies for which the model spectrum amplitudes are computed.

• pars (NDArray[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[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)

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[float]

Parameters:

• timestep (float)

• freqs (NDArray[float])

• amplitudes (NDArray[float])

derive_props(props)

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

Parameters:

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

get_par_nonlinear()

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[bool]

property name

neglog_prior(pars, *, deriv=0)

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

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

Return type:

list[NDArray[float]]

Parameters:

• pars (NDArray[float])

• deriv (int)

property npar

Return the number of parameters.

property par_scales: NDArray[float]

Return the scales of the parameters and the cost function.

sample_nonlinear_pars(rng, budget)

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[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)

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:

• freqs (NDArray[float]) – The frequencies for which the model spectrum amplitudes are computed.

• amplitudes (NDArray[float]) – The amplitudes of the spectrum.

• ndofs (NDArray[float]) – The number of degrees of freedom at each frequency.

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

• nonlinear_pars (NDArray[float]) – The values of the nonlinear parameters for which the basis functions are computed.

Return type:

NDArray[float]

Returns:

• linear_pars – The solved linear parameters.

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

valid(pars)

Return True when the parameters are within the feasible region.

Return type:

bool

Parameters:

pars (NDArray[float])

which_invalid(pars)

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

Return type:

NDArray[bool]

convert_pade022_debye(pars, covar)

Convert parameters and covariance from Pade(0,2;2) to the Debye model.

Parameters:

• pars (NDArray[float]) – The parameters of the Pade(0,2;2) model.

• covar (NDArray[float]) – The covariance matrix of the Pade(0,2;2) model.

Return type:

tuple[NDArray[float], NDArray[float]]

Returns:

• pars_debye – The DC level p0, the step height C_delta = p2/q2 - p0 and the relaxation time sqrt(q2)/(2 pi).

• pars_debye_covar – The covariance matrix of the Debye parameters (linear error propagation).

convert_pade022_lorentz(pars, covar)

Convert parameters and covariance from Pade(0,2;2) to Lorentz model.

Parameters:

• pars (NDArray[float]) – The parameters of the Pade(0,2;2) model.

• covar (NDArray[float]) – The covariance matrix of the Pade(0,2;2) model.

Return type:

tuple[NDArray[float], NDArray[float]]

Returns:

• pars_lorentz – The white noise level C0 = p2/q2, the amplitude C1 = pi (p0 - p2/q2) / sqrt(q2) and the exponential correlation time sqrt(q2)/(2 pi).

• pars_lorentz_covar – The covariance matrix of the Lorentz parameters (linear error propagation).

guess(freqs, ndofs, amplitudes, weights, model, rng, nonlinear_budget)

Guess initial values of the parameters for a model.

Parameters:

• freqs (NDArray[float]) – The frequencies for which the model spectrum amplitudes are computed.

• ndofs (NDArray[float]) – The number of degrees of freedom at each frequency.

• amplitudes (NDArray[float]) – The amplitudes of the spectrum.

• weights (NDArray[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.