SAS Functions

Module functions

This module defines classes representing SAS functions.

Currently there are two classes : Piecewise and Continuous. The first allows the SAS function to be specified as a set of breakpoints in a piecewise linear form of the cumulative distribution. The second allows any distribution specified as a scipy.stats rv_continuous class to be used as a SAS function. However, it will be converted into a piecewise form when the model is run.

Calling an instance of either class supplied with values of ST (as an array, list, or number) evaluates the CDF, and returns corresponding cumulative probabilities. The Piecewise.inv() method takes cumulative probabilities and returns values of ST.

class mesas.sas.functions.Piecewise(ST=None, P=None, nsegment=1, ST_max=1.0, ST_min=0.0, auto='uniform')[source]

Bases: _SASFunctionBase

Piecewise-linear SAS function defined by ST breakpoints.

The CDF is represented as a piecewise-linear curve through a set of (ST, P) pairs. Linear interpolation is used to evaluate the CDF and its inverse.

You can define a piecewise sas function in several different ways. The default is a single segment from 0 to 1:

>>>sas_fun = Piecewise() >>>sas_fun([0.1, 0.3, 0.9]) array([0.1, 0.3, 0.9])

The range of the single segment can be modified by providing ST_min and ST_max keywords

>>>sas_fun = Piecewise(ST_min=100, ST_max=110) >>>sas_fun([99, 101, 103, 109, 111]) array([0. , 0.1, 0.3, 0.9, 1. ])

The (ST,P) pairs defining the CDF are given by the ST and P attributes

>>>sas_fun.ST array([100. , 103.00391164, 104.17022005, 110. ]) >>>sas_fun.P array([0. , 0.33333333, 0.66666667, 1. ])

The segments defining the function are stored in the segment_list attribute as a numpy array:

>>>seglst = sas_fun.segment_list array([100. , 3.00391164, 1.1663084 , 5.82977995])

The first item in the array is ST_min, and each following value is the interval along the ST axis to the end of each segment.

Once created, the values of segment_list or ST can be modified and the other will update too. Values of P can also be changed.

Parameters:

ST (array_like or None, optional) – Explicit ST breakpoints. Must be length > 1 with ST[0] >= 0.
P (array_like or None, optional) – Cumulative probabilities at each breakpoint. Must start at 0, end at 1, and be non-decreasing. Length must match ST.
nsegment (int, optional) – Number of segments when ST is not provided. Default is 1.
ST_max (float, optional) – Upper bound of the SAS function support. Default is 1.0.
ST_min (float, optional) – Lower bound of the SAS function support. Default is 0.0.
auto (str, optional) – How to auto-generate breakpoints: 'uniform' (default) for equally spaced, 'random' for random recursive splitting.

property argsS: ndarray

ST breakpoints (alias for ST).

Type:: np.ndarray

property argsP: ndarray

Probability breakpoints (alias for P).

Type:: np.ndarray

property ST: ndarray

Cumulative age-ranked storage breakpoints.

Type:: np.ndarray

property P: ndarray

Cumulative probabilities at each ST breakpoint.

Type:: np.ndarray

property parameter_list: ndarray

Segment-length parameterisation of the SAS function.

Type:: np.ndarray

subdivided_copy(segment, split_frac=0.5)[source]

Return a new Piecewise instance with one segment split in two.

Parameters:

segment (int) – Zero-based index of the segment to split.
split_frac (float, optional) – Fraction along the segment where the split occurs. Default 0.5.

Returns:

A new SAS function with nsegment + 1 segments.

Return type:

Piecewise

inv(P)[source]

Evaluate the inverse CDF of the SAS function.

Parameters:: P (array_like) – Cumulative probabilities (scalar, list, or array).
Returns:: Corresponding ST values, same shape as input.
Return type:: np.ndarray

get_jacobian(dCdSj, index=None, mode='segment', logtransform=True)[source]

Compute a limited Jacobian of predicted concentrations.

If \(\vec{y} = f(\vec{x})\) the Jacobian is a matrix where each i, j entry is

Here, the \(y_i\) values are the predicted concentrations at each timestep, and the \(x_j\) are either the lengths of each piecewise segment along with ST_min (default), or the values of ST at the endpoints of each segment (with mode=’endpoint’).

The Jacobian is limited as it represents the sensitivity of the model predictions of output concentrations at each timestep to variations in the sas function at that timestep, and neglects the cumulative effects that changing the sas function would have on the state variables (ST and MS).

Parameters:

dCdSj (np.ndarray) – Jacobian of concentrations w.r.t. ST breakpoints from the Fortran solver, shape (n_timesteps, nsegment + 1).
index (np.ndarray or None, optional) – Timestep indices to include. Defaults to all timesteps.
mode (str, optional) – 'endpoint' returns derivatives w.r.t. ST values; 'segment' (default) w.r.t. segment lengths.
logtransform (bool, optional) – If True and mode='segment', return derivatives w.r.t. log-transformed segment lengths. Default is True.

Returns:

The Jacobian matrix.

Return type:

np.ndarray

class mesas.sas.functions.Continuous(use, func, func_kwargs, P=None, nsegment=25, ST_max=1.797693134862315e+308)[source]

Bases: _SASFunctionBase

SAS function backed by a continuous probability distribution.

Wraps a scipy.stats continuous distribution (or one of the built-in distributions: gamma, beta, kumaraswamy) as a SAS function. When used with use='scipy.stats', the continuous CDF is also discretised into a piecewise-linear lookup table for use by the Fortran solver.

Parameters:

use (str) – Either 'builtin' for optimised built-in distributions or 'scipy.stats' for arbitrary rv_continuous distributions.
func (str or rv_continuous) – Distribution name ('gamma', 'beta', 'kumaraswamy') when use='builtin', or a scipy.stats.rv_continuous instance / name string when use='scipy.stats'.
func_kwargs (dict[str, Any]) – Keyword arguments passed to the distribution to create a frozen instance (e.g. {'a': 2, 'loc': 0, 'scale': 100}).
P (array_like or None, optional) – Probability breakpoints for the piecewise lookup table. Only used when use='scipy.stats'.
nsegment (int, optional) – Number of segments in the lookup table when P is not given. Default is 25.
ST_max (float, optional) – Upper storage bound for the lookup table. Default is float('inf').

func

The frozen distribution instance.

Type:: scipy.stats.rv_frozen

property argsS: list[float] | ndarray

Storage-axis distribution parameters.

Type:: list or np.ndarray

property argsP: ndarray

Probability-axis distribution parameters.

Type:: np.ndarray

property func: scipy.stats.rv_frozen

The frozen distribution instance.

Type:: scipy.stats.rv_frozen

property ST: ndarray

Cumulative age-ranked storage breakpoints.

Type:: np.ndarray

property P: ndarray

Cumulative probabilities at each ST breakpoint.

Type:: np.ndarray

inv(P)[source]

Evaluate the inverse CDF of the SAS function.

Parameters:: P (array_like) – Cumulative probabilities.
Returns:: Corresponding ST values.
Return type:: np.ndarray

get_jacobian(dCdSj, index=None, mode='segment', logtransform=True)[source]

Compute a limited Jacobian of predicted concentrations.

Parameters:

dCdSj (np.ndarray) – Jacobian from the Fortran solver.
index (np.ndarray or None, optional) – Timestep indices. Defaults to all.
mode (str, optional) – 'segment' or 'endpoint'.
logtransform (bool, optional) – Log-transform segment derivatives. Default True.

Returns:

The Jacobian matrix.

Return type:

np.ndarray

Raises:

NotImplementedError – Always raised; not yet implemented for continuous functions.