stratage package

Submodules

stratage.geochron module

class stratage.geochron.Geochron(h, rv, dt, prob_threshold=1e-08, neighbors='all')[source]

Bases: object

Class for handling geochronologic constraints in stratigraphic sections.

model_weights(units)[source]

Calculate the weights for units and contacts based on their bounding by geochron constraints.

Parameters:

units (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

Returns:

Tuple[List[List[float]], List[List[float]]] –

  • alpha: 2D list (NxL) of weights for each of N units and L pairs of geochron constraints.

  • beta: 2D list (MxL) of weights for each of M contacts and L pairs of geochron constraints.

Raises:

AssertionError – If any rows in alpha or beta are entirely zero, indicating unconstrained model parameters.

Notes

  • Run after stratagemc.trim_units().

plot_constraints(ax=None, tol=3, scale=1, **kwargs)[source]

Plot the geochronologic constraints.

Parameters:

  • ax (matplotlib.axes.Axes) – The axes on which to plot the constraints. If None, a new figure is created.

  • tol (int) – The number of orders of magnitude below the maximum pdf value to plot.

  • scale (float) – The scale factor for the pdfs.

  • **kwargs – Additional keyword arguments passed to the matplotlib fill_between function.

Returns:

  • matplotlib.axes._axes.Axes – The axes on which the constraints are plotted.

  • list – The handles for the plotted constraints.

plot_increment_pdfs(ax=None, tol=3, scale=1)[source]

Plot the time increment pdfs.

Parameters:

  • ax (matplotlib.axes.Axes) – The axes on which to plot the increment pdfs. If None, a new figure is created.

  • tol (int) – The number of orders of magnitude below the maximum pdf value to plot.

  • scale (float) – The scale factor for the pdfs.

Returns:

matplotlib.axes.Axes – The axes on which the increment pdfs are plotted.

straddling(height)[source]

Returns the indices of pairs of constraints that straddle a given height. :param height: The height in section. :type height: float

Returns:

array-like – The indices (into self.lower_idx, self.upper_idx) that straddle the given height.

stratage.stratage module

class stratage.stratage.AgeModel(units, geochron, sed_rates_prior, hiatuses_prior, ls_kwargs=None)[source]

Bases: object

Age model object for Bayesian inference of sedimentation rates and hiatuses.

units

nx2 array of unit bottom and top heights for n units.

Type:

numpy.ndarray

geochron

Geochron object containing geochron constraints.

Type:

geochron.Geochron

sed_rates_prior

Prior distribution for sedimentation rates. Must be valid as dist argument to pymc.CustomDist(dist=dist). Signature is sed_rate_prior(size=size).

Type:

function

hiatuses_prior

Prior distribution for hiatuses. Must be valid as dist argument to pymc.CustomDist(dist=dist). Signature is hiatus_prior(size=size).

Type:

function

n_units

Number of units.

Type:

int

n_contacts

Number of contacts.

Type:

int

units_trim

Trimmed unit heights after adjusting for the top and bottom units.

Type:

numpy.ndarray

sed_rates_ls

Sedimentation rates from least squares model.

Type:

numpy.ndarray

hiatuses_ls

Hiatuses from least squares model.

Type:

numpy.ndarray

model

PyMC model object for Bayesian inference.

Type:

pymc.Model

vars_list

List of variables in the pymc.model.

Type:

list

static fit_absolute_age(ii, sed_rates_post, hiatuses_post, units_trim, geochron, h, method)[source]

Fit the age model for the ii-th sample of the posterior.

Static method to work with joblib.Parallel.

Parameters:

  • ii (int) – Index of the posterior sample.

  • sed_rates_post (ndarray) – Sedimentation rates for each unit, ii-th posterior sample.

  • hiatuses_post (ndarray) – Hiatuses between units, ii-th posterior sample.

  • units_trim (ndarray) – Trimmed unit heights after adjusting for the top and bottom units.

  • geochron (Geochron) – Geochron object containing geochron constraints.

  • h (arraylike) – Heights at which to evaluate the age model.

  • method (str) – Method for optimization in fit_floating_model(). ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood.

Returns:

ndarray – Ages at the given height(s).

sample(draws=1000, **kwargs)[source]

Sample posterior distribution of sedimentation rates and hiatus durations.

The output of this function must be transformed to age models.

Parameters:

draws (int, optional) – Number of posterior draws. Defaults to 1000.

Returns:

arviz.InferenceData – ArviZ InferenceData object containing the MCMC trace.

sample_prior(draws=100)[source]

Sample prior predictive distribution of age models.

Parameters:

draws (int, optional) – Number of prior predictive draws. Defaults to 100.

Returns:

list – List of prior predictive samples of times; each element is a nx2 array of unit bottom and top times for n units.

trace2ages(trace, h, n_posterior=None, n_jobs=1, method='random')[source]

Transform MCMC trace to age models.

Parameters:

  • trace (arviz.InferenceData) – ArviZ InferenceData object containing the MCMC trace.

  • h (arraylike) – Heights at which to evaluate the age model.

  • n_posterior (int, optional) – Number of posterior samples. Defaults to None.

  • n_jobs (int, optional) – Number of parallel jobs. Defaults to 1. If 1, no parallelization is used. Uses joblib.Parallel for parallelization.

  • method (str, optional) – Method for optimization in fit_floating_model(). ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood. Defaults to ‘random’.

Returns:

  • list – List of age models; each element is a nx2 array of unit bottom and top times for n units.

  • numpy.ndarray – Array of age models at heights h; each row is an age model. Shape is (n_posterior, len(h)). Only returned if h is not None.

stratage.stratage.DT_logp_l_gen(pdt)[source]

Generates a log-probability function for a numerical time increment distribution. Returned function interpolates a log probability and takes as inputs time increment(s) at which to evaluate as well as (required) time increment coordinates corresponding to the numerical pdf associated with the function. Sped up with numba.

Parameters:

pdt (array-like) – The probability density function values for the time increments.

Returns:

numba.core.registry.CPUDispatcher – A function that computes the log-probability of a given time increment. The returned function, DT_logp, takes the following parameters:

dt_query (float): The time increment for which the log-probability is to be computed. dt (array-like): The array of time increments corresponding to the probability density function values. Must be same length as pdt.

The DT_logp function computes the log-probability by interpolating the probability density function values and applying sigmoid functions to ensure the values are within the valid range.

stratage.stratage.age(times, units, heights)[source]

Interpolate time at given heights in strat. Times at contact heights are returned as the age of the top of the unit below the contact.

Parameters:

  • times (arraylike) – nx2 array of unit bottom and top times for n units.

  • units (arraylike) – nx2 array of unit bottom and top heights for n units.

  • heights (arraylike) – Height(s) at which to evaluate the age model.

Returns:

array – Age(s) at the given height(s).

stratage.stratage.age_depth(units, times)[source]

Convert units, times arrays to age-depth curve

Parameters:

  • units (numpy.ndarray) – An nx2 array of bottom and top elevations for each of n units.

  • times (numpy.ndarray) – An nx2 array of bottom and top times for each of n units.

Returns:

  • numpy.ndarray – Vector of ages for each unit.

  • numpy.ndarray – Vector of heights for each unit.

stratage.stratage.fit_floating_model(sed_rates, hiatuses, units, geochron, tol=1e-06, method='random')[source]

Fit a floating age model to geochronologic constraints

Two methods are implemented. One (method=’max’) maximizes the likelihood of alignment with the geochronologic constraints by using scipy.optimize.minimize_scalar to minimize the negative log likelihood of a time offset with respect to the distributions in absolute time of the geochronologic constraints.

The other (method=’random’) samples from the total alignment likelihood and selects a random offset that is consistent with the likelihood. This method is the default.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit.

  • hiatuses (arraylike) – Hiatuses between units. Must be len(sed_rates)-1.

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

  • geochron (geochron.Geochron) – Geochron object containing geochron constraints.

  • tol (float, optional) – Tolerance for bounds on optimization. Defaults to 1e-6.

  • method (str, optional) – Method for optimization. ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood. Defaults to ‘random’.

Returns:

ndarray – nx2 array of unit bottom and top ages for n units in absolute time

stratage.stratage.floating_age(sed_rates, hiatuses, units, heights)[source]

Floating ages at heights in strat. Assumes t=0 at base of section. Times at contact heights are returned as the age of the top of the unit below the contact.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit.

  • hiatuses (arraylike) – Hiatuses between units. Must be len(sed_rates)-1.

  • units (arraylike) – nx2 array of unit bottom and top heights for n units.

  • heights (arraylike | float) – Height(s) at which to evaluate the age model.

Returns:

array – Age(s) at the given height(s).

stratage.stratage.geochron_height_check(unit_heights, geochron_heights)[source]

Ensure that geochron heights are within units and not at contacts. Also ensure that geochron heights are within the section. Run before trimming section to the top and bottom of the geochron constraints.

Parameters:

  • unit_heights (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

  • geochron_heights (arraylike) – Heights of geochron constraints in the section

Raises:

  • AssertionError – If any of the geochron heights are below the section or above the section.

  • AssertionError – If any of the geochron heights are at contacts.

stratage.stratage.get_times(sed_rates, hiatuses, units)[source]

Floating times array for units.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit

  • hiatuses (arraylike) – Hiatuses between units

  • units (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

Returns:

times (numpy.ndarray) – Floating times array for units, same shape as units

stratage.stratage.loglike_gen(geochron, units)[source]

Returns log-likelihood object for use in PyMC.CustomDist as a likelihood for posterior sampling. See subfunction eps_logp for details on the log-likelihood computation.

Parameters:

  • geochron (geochron.Geochron) – Geochron object containing geochronologic constraints.

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

Returns:

LogLike – PyTensor Op class for likelihood evaluation.

stratage.stratage.model_ls(units, geochron, sed_rate_bounds=None, hiatus_bounds=None)[source]

Least squares age model fitting.

Parameters:

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

  • geochron (geochron.Geochron) – Geochron object containing geochron constraints.

  • sed_rate_bounds (list, optional) – Bounds on sedimentation rates. Defaults to None.

  • hiatus_bounds (list, optional) – Bounds on hiatuses. Defaults to None.

Returns:

  • sed_rates (numpy.ndarray) – Sedimentation rates for each unit.

  • hiatuses (numpy.ndarray) – Hiatuses between units.

stratage.stratage.randlike_gen(geochron, units)[source]

Generate function for generating random draws from likelihood. Permits generating prior_predictive samples from a model with the CustomDist likelihood.

Parameters:

  • geochron (geochron.Geochron) – Geochron object containing geochronologic constraints.

  • units (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

Returns:

function – Function for generating random draws from likelihood. The returned function, eps_rand, takes the following parameters:

params (array-like): Array of parameters; concatenation of sedimentation rates and hiatuses.
rng (numpy.random.Generator): Random number generator.
size (int or tuple of ints, optional): Output shape. Default is None. Last dimension must be the number of pairs of geochron constraints.

The function returns a numpy.ndarray of random draws.

stratage.stratage.sigmoid(t, scale=0.001)[source]

Sigmoid function for numerical stability.

Parameters:

  • t (float | arraylike) – Value to apply sigmoid to.

  • scale (float, optional) – Transition width scale of sigmoid. Defaults to 0.001.

Returns:

float | arraylike – Sigmoid applied to t.

stratage.stratage.trim_units(unit_heights, geochron_heights)[source]

Trim the top and bottom of the section to the top and bottom of the geochron constraints. Run after geochron_height_check.

Parameters:

  • unit_heights (numpy.ndarray) – Bottom and top heights of units in section. 2d array-like (nx2) of bottom and top heights for each of n units

  • geochron_heights (arraylike) – Heights of geochron constraints in section

Returns:

numpy.ndarray – Trimmed unit heights after adjusting for the top and bottom units

Raises:

None –

Module contents

class stratage.AgeModel(units, geochron, sed_rates_prior, hiatuses_prior, ls_kwargs=None)[source]

Bases: object

Age model object for Bayesian inference of sedimentation rates and hiatuses.

units

nx2 array of unit bottom and top heights for n units.

Type:

numpy.ndarray

geochron

Geochron object containing geochron constraints.

Type:

geochron.Geochron

sed_rates_prior

Prior distribution for sedimentation rates. Must be valid as dist argument to pymc.CustomDist(dist=dist). Signature is sed_rate_prior(size=size).

Type:

function

hiatuses_prior

Prior distribution for hiatuses. Must be valid as dist argument to pymc.CustomDist(dist=dist). Signature is hiatus_prior(size=size).

Type:

function

n_units

Number of units.

Type:

int

n_contacts

Number of contacts.

Type:

int

units_trim

Trimmed unit heights after adjusting for the top and bottom units.

Type:

numpy.ndarray

sed_rates_ls

Sedimentation rates from least squares model.

Type:

numpy.ndarray

hiatuses_ls

Hiatuses from least squares model.

Type:

numpy.ndarray

model

PyMC model object for Bayesian inference.

Type:

pymc.Model

vars_list

List of variables in the pymc.model.

Type:

list

static fit_absolute_age(ii, sed_rates_post, hiatuses_post, units_trim, geochron, h, method)[source]

Fit the age model for the ii-th sample of the posterior.

Static method to work with joblib.Parallel.

Parameters:

  • ii (int) – Index of the posterior sample.

  • sed_rates_post (ndarray) – Sedimentation rates for each unit, ii-th posterior sample.

  • hiatuses_post (ndarray) – Hiatuses between units, ii-th posterior sample.

  • units_trim (ndarray) – Trimmed unit heights after adjusting for the top and bottom units.

  • geochron (Geochron) – Geochron object containing geochron constraints.

  • h (arraylike) – Heights at which to evaluate the age model.

  • method (str) – Method for optimization in fit_floating_model(). ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood.

Returns:

ndarray – Ages at the given height(s).

sample(draws=1000, **kwargs)[source]

Sample posterior distribution of sedimentation rates and hiatus durations.

The output of this function must be transformed to age models.

Parameters:

draws (int, optional) – Number of posterior draws. Defaults to 1000.

Returns:

arviz.InferenceData – ArviZ InferenceData object containing the MCMC trace.

sample_prior(draws=100)[source]

Sample prior predictive distribution of age models.

Parameters:

draws (int, optional) – Number of prior predictive draws. Defaults to 100.

Returns:

list – List of prior predictive samples of times; each element is a nx2 array of unit bottom and top times for n units.

trace2ages(trace, h, n_posterior=None, n_jobs=1, method='random')[source]

Transform MCMC trace to age models.

Parameters:

  • trace (arviz.InferenceData) – ArviZ InferenceData object containing the MCMC trace.

  • h (arraylike) – Heights at which to evaluate the age model.

  • n_posterior (int, optional) – Number of posterior samples. Defaults to None.

  • n_jobs (int, optional) – Number of parallel jobs. Defaults to 1. If 1, no parallelization is used. Uses joblib.Parallel for parallelization.

  • method (str, optional) – Method for optimization in fit_floating_model(). ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood. Defaults to ‘random’.

Returns:

  • list – List of age models; each element is a nx2 array of unit bottom and top times for n units.

  • numpy.ndarray – Array of age models at heights h; each row is an age model. Shape is (n_posterior, len(h)). Only returned if h is not None.

stratage.DT_logp_l_gen(pdt)[source]

Generates a log-probability function for a numerical time increment distribution. Returned function interpolates a log probability and takes as inputs time increment(s) at which to evaluate as well as (required) time increment coordinates corresponding to the numerical pdf associated with the function. Sped up with numba.

Parameters:

pdt (array-like) – The probability density function values for the time increments.

Returns:

numba.core.registry.CPUDispatcher – A function that computes the log-probability of a given time increment. The returned function, DT_logp, takes the following parameters:

dt_query (float): The time increment for which the log-probability is to be computed. dt (array-like): The array of time increments corresponding to the probability density function values. Must be same length as pdt.

The DT_logp function computes the log-probability by interpolating the probability density function values and applying sigmoid functions to ensure the values are within the valid range.

stratage.age(times, units, heights)[source]

Interpolate time at given heights in strat. Times at contact heights are returned as the age of the top of the unit below the contact.

Parameters:

  • times (arraylike) – nx2 array of unit bottom and top times for n units.

  • units (arraylike) – nx2 array of unit bottom and top heights for n units.

  • heights (arraylike) – Height(s) at which to evaluate the age model.

Returns:

array – Age(s) at the given height(s).

stratage.age_depth(units, times)[source]

Convert units, times arrays to age-depth curve

Parameters:

  • units (numpy.ndarray) – An nx2 array of bottom and top elevations for each of n units.

  • times (numpy.ndarray) – An nx2 array of bottom and top times for each of n units.

Returns:

  • numpy.ndarray – Vector of ages for each unit.

  • numpy.ndarray – Vector of heights for each unit.

stratage.fit_floating_model(sed_rates, hiatuses, units, geochron, tol=1e-06, method='random')[source]

Fit a floating age model to geochronologic constraints

Two methods are implemented. One (method=’max’) maximizes the likelihood of alignment with the geochronologic constraints by using scipy.optimize.minimize_scalar to minimize the negative log likelihood of a time offset with respect to the distributions in absolute time of the geochronologic constraints.

The other (method=’random’) samples from the total alignment likelihood and selects a random offset that is consistent with the likelihood. This method is the default.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit.

  • hiatuses (arraylike) – Hiatuses between units. Must be len(sed_rates)-1.

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

  • geochron (geochron.Geochron) – Geochron object containing geochron constraints.

  • tol (float, optional) – Tolerance for bounds on optimization. Defaults to 1e-6.

  • method (str, optional) – Method for optimization. ‘max’ for maximum likelihood, ‘random’ for random sampling of alignment likelihood. Defaults to ‘random’.

Returns:

ndarray – nx2 array of unit bottom and top ages for n units in absolute time

stratage.floating_age(sed_rates, hiatuses, units, heights)[source]

Floating ages at heights in strat. Assumes t=0 at base of section. Times at contact heights are returned as the age of the top of the unit below the contact.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit.

  • hiatuses (arraylike) – Hiatuses between units. Must be len(sed_rates)-1.

  • units (arraylike) – nx2 array of unit bottom and top heights for n units.

  • heights (arraylike | float) – Height(s) at which to evaluate the age model.

Returns:

array – Age(s) at the given height(s).

stratage.geochron_height_check(unit_heights, geochron_heights)[source]

Ensure that geochron heights are within units and not at contacts. Also ensure that geochron heights are within the section. Run before trimming section to the top and bottom of the geochron constraints.

Parameters:

  • unit_heights (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

  • geochron_heights (arraylike) – Heights of geochron constraints in the section

Raises:

  • AssertionError – If any of the geochron heights are below the section or above the section.

  • AssertionError – If any of the geochron heights are at contacts.

stratage.get_times(sed_rates, hiatuses, units)[source]

Floating times array for units.

Parameters:

  • sed_rates (arraylike) – Sedimentation rates for each unit

  • hiatuses (arraylike) – Hiatuses between units

  • units (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

Returns:

times (numpy.ndarray) – Floating times array for units, same shape as units

stratage.loglike_gen(geochron, units)[source]

Returns log-likelihood object for use in PyMC.CustomDist as a likelihood for posterior sampling. See subfunction eps_logp for details on the log-likelihood computation.

Parameters:

  • geochron (geochron.Geochron) – Geochron object containing geochronologic constraints.

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

Returns:

LogLike – PyTensor Op class for likelihood evaluation.

stratage.model_ls(units, geochron, sed_rate_bounds=None, hiatus_bounds=None)[source]

Least squares age model fitting.

Parameters:

  • units (ndarray) – nx2 array of unit bottom and top heights for n units.

  • geochron (geochron.Geochron) – Geochron object containing geochron constraints.

  • sed_rate_bounds (list, optional) – Bounds on sedimentation rates. Defaults to None.

  • hiatus_bounds (list, optional) – Bounds on hiatuses. Defaults to None.

Returns:

  • sed_rates (numpy.ndarray) – Sedimentation rates for each unit.

  • hiatuses (numpy.ndarray) – Hiatuses between units.

stratage.randlike_gen(geochron, units)[source]

Generate function for generating random draws from likelihood. Permits generating prior_predictive samples from a model with the CustomDist likelihood.

Parameters:

  • geochron (geochron.Geochron) – Geochron object containing geochronologic constraints.

  • units (numpy.ndarray) – Bottom and top heights of units in section. 2d array (nx2) of bottom and top heights for each of n units

Returns:

function – Function for generating random draws from likelihood. The returned function, eps_rand, takes the following parameters:

params (array-like): Array of parameters; concatenation of sedimentation rates and hiatuses.
rng (numpy.random.Generator): Random number generator.
size (int or tuple of ints, optional): Output shape. Default is None. Last dimension must be the number of pairs of geochron constraints.

The function returns a numpy.ndarray of random draws.

stratage.sigmoid(t, scale=0.001)[source]

Sigmoid function for numerical stability.

Parameters:

  • t (float | arraylike) – Value to apply sigmoid to.

  • scale (float, optional) – Transition width scale of sigmoid. Defaults to 0.001.

Returns:

float | arraylike – Sigmoid applied to t.

stratage.trim_units(unit_heights, geochron_heights)[source]

Trim the top and bottom of the section to the top and bottom of the geochron constraints. Run after geochron_height_check.

Parameters:

  • unit_heights (numpy.ndarray) – Bottom and top heights of units in section. 2d array-like (nx2) of bottom and top heights for each of n units

  • geochron_heights (arraylike) – Heights of geochron constraints in section

Returns:

numpy.ndarray – Trimmed unit heights after adjusting for the top and bottom units

Raises:

None –