demestats Core Functions API Reference

Model Constraint

demestats.event_tree.EventTree.variables()

List out all of the demestats variables corresponding to a demes.Graph.

Returns:

All associated variables given a demes graph

Return type:

Sequence[Variable]

Notes

To understand what a demestats and demes path are, please refer to the Notation documentation.

To use this function you must first create an EventTree object.

Example:

et = EventTree(demo.to_demes())
print(et.variables)
demestats.event_tree.EventTree.variable_for(self, path: tuple[str | int, ...]) → Set[tuple[str | int, ...]]

Get the variable corresponding to a demestats path

Parameters:

path (Path) – a tuple representing a demestats path

Returns:

The associated demestats variable given a path

Return type:

Variable

Notes

To understand what a demestats path is, please refer to the Notation documentation.

To use this function you must create an EventTree.

Example:

et = EventTree(demo.to_demes())

# EventTree.variable_for, example taken from momi3 Tutorial
parameters = [
    ('demes', 0, 'epochs', 0, 'end_size'), # The ancestral population size
    ('migrations', 0, 'rate'), # Rate of migration from P0 to P1
    ('demes', 0, 'epochs', 0, 'end_time') # Time of divergence
]

momi3_parameters = [et.variable_for(param) for param in parameters]
demestats.constr.constraints_for(et: EventTree, *vars_: Set[tuple[str | int, ...]]) → ConstraintSet[source]

Return a list of constraints for the given variables.

Parameters:

  • et (EventTree) – The event tree to extract constraints from

  • var (Variable) – The variables for which to extract constraints, which must exist in et.variables.

Returns:

A dictionary containing the equality and inequality constraints, as well as a mapping of columns of the constraint matrices to the variables.

Return type:

ConstraintSet

Notes

Example:

et = EventTree(demo.to_demes())
variables = et.variables
final_constraints = constraints_for(et, *variables)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

demestats.constr.print_constraints(constraint_dict, variable_list)[source]

Print out all equality and inequality constraints in an interpretable way.

Parameters:

  • constraint_dict (dict) – A dictionary of equality and inequality constraints

  • variable (list) – List of associated variable names

Notes

The constraint_dict is obtained from the output of constraints_for function.

Example:

et = EventTree(demo.to_demes())

# See example in the tutorial
parameters = [
    ('demes', 0, 'epochs', 0, 'end_size'), # The ancestral population size
    ('migrations', 0, 'rate'), # Rate of migration from P0 to P1
    ('demes', 0, 'epochs', 0, 'end_time') # Time of divergence
]

variable_list = [et.variable_for(param) for param in parameters]
constraint = constraints_for(et, *variable_list)
print_constraints(constraint, variable_list)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.constr.constraints_for

demestats.fit.util.alternative_constraint_rep(A, b)[source]

Returns an alternative representation of inequality constraints with a lower and upper bound. Depending on the numerical optimizer one would like to use, sometimes it’s more preferable to express inequality constraints explicitly with a lower and upper bound. The input for the function are inequality constraints of the form Ax <= b which is exactly what the demestats.constr.constraints_for function returns.

Parameters:

  • A (array_like) – Coefficients for inequalities of the form Ax <= b

  • b (array_list) – Values for inequalities of the form Ax <= b

Notes

Example:

# See tutorial for a detailed example
parameters = [
    ('demes', 0, 'epochs', 0, 'end_size'), # The ancestral population size
    ('migrations', 0, 'rate'), # Rate of migration from P0 to P1
    ('demes', 0, 'epochs', 0, 'end_time') # Time of divergence
]

momi3_parameters = [et.variable_for(param) for param in parameters]
constraint = constraints_for(et, *momi3_parameters)
G, h = constraint["ineq"]
A_alt, ub_alt, lb_alt = alternative_constraint_rep(G, h)
print(A_alt)
print("lower bound: ", lb_alt)
print("upper bound: ", ub_alt)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.constr.constraint_for

demestats.fit.util.modify_constraints_for_equality(constraint, indices_for_equality)[source]

Returns a modified version of the input constraint where all parameters associated with indicies in indicies_for_equality will now have an equality constraint.

Parameters:

  • constraint (dict) – A dictionary of equality and inequality constraints

  • indices_for_equality (list) – List of tuples, where each tuple are the indices of parameters you want to impose an equality constraint

Returns:

A modified constraint with the new equality constraints

Return type:

dict

Notes

Example:

# See tutorial for a detailed example
parameters = [
    ('demes', 0, 'epochs', 0, 'end_size'), # The ancestral population size
    ('migrations', 0, 'rate'), # Rate of migration from P0 to P1
    ('demes', 0, 'epochs', 0, 'end_time') # Time of divergence
]

momi3_parameters = [et.variable_for(param) for param in parameters]
constraint = constraints_for(et, *momi3_parameters)

# new_constraint will have the 2nd and 3rd variable be constrained to be equal
new_constraint = modify_constraints_for_equality(constraint, [(1, 2)])
print(new_constraint)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.constr.constraint_for

momi3 (SFS) Statistics

demestats.sfs.ExpectedSFS(demo: Graph, num_samples: dict[str, Int[Array, ' '] | Int[ndarray, ' '] | number | int], prune: Mapping[str, int] | Sequence[tuple[str, int, tuple[str | int, ...] | int | None]] | None = None) → None[source]

Build an ExpectedSFS object that can be later used to compute the full expected site frequency spectrum or the projected site frequency spectrum.

Parameters:

  • demo (demes.Graph) – A demes graph

  • num_samples (dict) – A dictionary specifying how many haploids per population to use to compute the expected SFS. The name of the populations must match the exact names used in demo.

  • prune (mapping or sequence, optional) – Optional manual downsampling events. Provide either a mapping {deme_name: m} to downsample directly above leaves, or a sequence of (deme_name, m[, at]) tuples where at is a node id or a demes time path to insert the downsample event above that node.

Returns:

an ExpectedSFS object used to compute the expected site frequency spectrum

Return type:

ExpectedSFS

Notes

From a user perspective, understanding the underlying structure of an ExpectedSFS object is not necessary. The only functions that a user would use is ExpectedSFS.__call__ which computes the full expected site frequency spectrum and ExpectedSFS.tensor_prod which computes the projected site frequency spectrum.

Example:

ESFS = ExpectedSFS(demo.to_demes(), num_samples=afs_samples)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.sfs.ExpectedSFS.__call__, demestats.sfs.ExpectedSFS.tensor_prod

demestats.sfs.ExpectedSFS.__call__(self, params: dict[Set[tuple[str | int, ...]], Shaped[Array, ''] | Shaped[ndarray, ''] | bool | number | bool | int | float | complex] = {}) → Float[Array, '*T']

Computes the full expected site frequency spectrum under a given set of model parameters and values

Parameters:

params (dict) – A dictionary of model parameters and their value

Returns:

An array of float values that represent the full expected site frequency spectrum

Return type:

Float[Array]

Notes

You must first construct an ExpectedSFS object. See the ExpectedSFS API.

Example:

ESFS = ExpectedSFS(demo.to_demes(), num_samples=afs_samples)
params = {param_key: val}
esfs = ESFS(params)

Please refer to the tutorial for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.sfs.ExpectedSFS

demestats.loglik.sfs_loglik.prepare_projection(afs, afs_samples, sequence_length, num_projections, seed)[source]

Creates the specified number of random projection vectors and appropriate inputs for the Einstein summation for tensor operations that are used in ExpectedSFS.tensor_prod.

Parameters:

  • afs (array_like) – Observed allele frequency spectrum

  • esfs (array_like) – Expected allele frequency spectrum. Must be the same shape as afs

  • sequence_length (int, optional) – Total number of sites in the sequence. Required if theta is given

  • num_projections (int) – Number of desired random projection vectors

  • seed (int) – Seed for reproducibility

Returns:

  • dict – dictionary of random projection vectors

  • str – string containing axes names separated by commas for Einstein summation

  • list – list containing a dictionary specifying number of haploids per population and the afs. This list is a necessary input for the Einstein summation

Notes

proj_dict contains the random projection vectors that define the low-dimensional subspace for approximating the full expected SFS, einsum_str is a string specifying the Einstein summation for tensor operations, and input_arrays are preprocessed arrays that serve as inputs to the jax.numpy.einsum call, optimized for JAX’s just-in-time compilation

Example:

proj_dict, einsum_str, input_arrays = prepare_projection(afs, afs_samples, sequence_length, num_projections, seed)

Please refer to Random Projection section for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.loglik.sfs_loglik.prepare_projection

demestats.sfs.ExpectedSFS.tensor_prod(self, X: PyTree[jaxtyping.Shaped[Array, 'B ?D'] | jaxtyping.Shaped[ndarray, 'B ?D'], 'T'], params: dict[Set[tuple[str | int, ...]], Shaped[Array, ''] | Shaped[ndarray, ''] | bool | number | bool | int | float | complex] = {}) → Float[Array, 'B']

Computes the projected expected site frequency spectrum under a given set random projection vectors and model parameters. A tensor product operation between the random projections and the expected site frequency spectrum is applied to obtain the projected SFS. To obtain the appropriate projection vectors, one can use the function demestats.loglik.sfs_loglik.prepare_projection.

Parameters:

  • X (dict) – A dictionary of random projection vectors

  • params (dict) – A dictionary of model parameters and their value

Returns:

An array of float values that represent the projected expected site frequency spectrum

Return type:

Float[Array]

Notes

You must first construct an ExpectedSFS object. See the ExpectedSFS API.

Example:

proj_dict, einsum_str, input_arrays = prepare_projection(afs, afs_samples, sequence_length, num_projections, seed)
esfs_obj = ExpectedSFS(demo.to_demes(), num_samples=afs_samples)
lowdim_esfs = esfs_obj.tensor_prod(proj_dict, paths)

Please refer to Random Projection section for a specific example, the above provided codes are just outlines of how to call on the functions.

See also

demestats.loglik.sfs_loglik.prepare_projection

ICR Statistics

CCR Statistics

demestats.ccr.exact.CCRCurve(demo: Graph, k: int) → None[source]

Build an CCRCurve object that can be later used to evaluate the cross-coalescence rate (CCR) through time for a demographic model.

Parameters:

  • demo (demes.Graph) – A demes graph describing the demographic history.

  • k (int) – The number of sampled lineages used to define the CCR curve.

Returns:

An CCRCurve object that can be called directly on a time grid, sampling configuration, and parameter value.

Return type:

CCRCurve

Notes

From a user perspective, understanding the underlying structure of an CCRCurve object is not necessary. The only function that a user would use is __call__, which evaluates the CCR on a grid of time points given a sampling configuration.

One can choose between two computational backends depending on the size of the problem, see “CCR: Exact vs Mean-Field” in the CCR tutorial.

::

# pairwise coalescence uses k = 10 ccr = CCRCurve(demo.to_demes(), k=10)

See also

demestats.ccr.CCRMeanFieldCurve

demestats.ccr.mf.CCRMeanFieldCurve(demo: Graph, k: int) → None[source]

Build an CCRMeanFieldCurve object that can be later used to evaluate the mean-field approximation of the cross-coalescence rate (CCR) through time for a demographic model.

Parameters:

  • demo (demes.Graph) – A demes graph describing the demographic history.

  • k (int) – The number of sampled lineages used to define the CCR curve.

Returns:

An CCRMeanFieldCurve object that can be called directly on a time grid, sampling configuration, and parameter value to construct the mean-field approximation of the CCR curve.

Return type:

CCRMeanFieldCurve

Notes

From a user perspective, understanding the underlying structure of an CCRMeanFieldCurve object is not necessary. The only function that a user would use is __call__, which evaluates the CCR on a grid of time points given a sampling configuration.

One can choose between two computational backends depending on the size of the problem, see “CCR: Exact vs Mean-Field” in the CCR tutorial.

::

# pairwise coalescence uses k = 10 ccr_mf = CCRMeanFieldCurve(demo.to_demes(), k=10)

See also

demestats.ccr.CCRMeanFieldCurve.__call__, demestats.ccr.CCRCurve