stratage package

Submodules

stratage.geochron module

class stratage.geochron.Geochron(h, rv, nt=500, ndt=100, prob_threshold=0.0001, neighbors='all', n_jobs=1)

Bases: object

Class for handling geochronologic constraints in stratigraphic sections.

Parameters:

• h (array-like) – List of stratigraphic heights at which temporal constraints are measured. The coordinate increases up section.

• rv (array-like) – List of random variables representing the temporal constraints at each height. The members of this list must be objects similar to stats.dist objects with the following methods: pdf, ppf (inverse of cdf).

• nt (int, optional) – Number of points for computing numerical pdfs from constraints. Defaults to 500.

• ndt (int, optional) – Number of points for computing numerical time increment pdfs. Defaults to 100. This value should be smaller than nt to avoid noisiness in the time increment pdfs.

• prob_threshold (float, optional) – Probability threshold for determining the temporal range over which to grid. Defaults to 1e-4.

• neighbors (str or int, optional) – Number of neighbors for pairing constraints. Defaults to ‘all’. Valid options are ‘all’ or integers up to n_constraints-1. All pairs are considered in the ‘all’ method. For integer values n, only n neighboring pairs are considered.

• n_jobs (int, optional) – Number of jobs to run in parallel for computing time increment pdfs. Defaults to 1

h

Sorted list of stratigraphic heights.

Type:

list

rv

Sorted list of random variables representing the temporal constraints.

Type:

list

n_constraints

Number of constraints.

Type:

int

n_pairs

Number of pairs of constraints.

Type:

int

nt

Number of points for computing numerical pdfs from constraints.

Type:

int

ndt

Number of points for computing numerical time increment pdfs.

Type:

int

lower_idx

List of indices into h and rv for the lower constraint in each pair.

Type:

list

upper_idx list

List of indices into h and rv for the upper constraint in each pair.

model_weights(units)

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:

• alpha (numpy.ndarray) – 2D array (NxL) of weights for each of N units and L pairs of geochron constraints.

• beta (numpy.ndarray) – 2D array (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, plotstyle=None, fillstyle=None)

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.

• plotstyle (dict, optional) – Additional keyword arguments passed to the matplotlib plot function for visualizing constraint pdfs, defaults to None. If None, the default plot style is used.

• fillstyle (dict, optional) – Additional keyword arguments passed to the matplotlib fill_between function for visualizing constraint pdfs, defaults to None. If None, the default fill style is used.

Returns:

• ax (matplotlib.axes.Axes) – The axes on which the constraints are plotted.

• h (list) – The handles for the plotted constraints. Each element of the list is a list of handles for the fill and line plots for a single constraint.

plot_increment_pdfs(ax=None, tol=3, scale=1)

Plot the time increment pdfs.

Parameters:

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

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

• scale (float, optional) – The scale factor for the pdfs. Defaults to 1.

Returns:

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

straddling(height)

Returns the indices of pairs of constraints that straddle a given height.

Parameters:

height (float) – The height in section.

Returns:

idx (numpy.ndarray) – 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)

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)

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)

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)

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

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)

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)

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)

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

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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

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)

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)

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)

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

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)

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)

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)

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)

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)

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)

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)

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)

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 –