Skip to content

Commit

Permalink
MAINT: Deprecate calc_*_from_diags functions (#111)
Browse files Browse the repository at this point in the history
* Move main code for calculating state 
probabilities into `calc_state_probs`, add 
`dir_edges` parameter, and add 
`DeprecationWarning` to 
`calc_state_probs_from_diags`.

* Implement the same changes as 
above for `calc_net_cycle_flux` and 
`calc_net_cycle_flux_from_diags`

* Update tests to use the updated
versions of `calc_state_probs` and 
`calc_net_cycle_flux`

* Add checks for `DeprecationWarning`s
to `test_function_inputs` for
`calc_state_probs_from_diags` and
`calc_net_cycle_flux_from_diags`

* Update docstrings for all functions 
using directional edges
  • Loading branch information
nawtrey authored Aug 15, 2024
1 parent 014d721 commit 2a7cc47
Show file tree
Hide file tree
Showing 2 changed files with 157 additions and 110 deletions.
232 changes: 135 additions & 97 deletions kda/calculations.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
.. footbibliography::
"""

import warnings
import math
import numpy as np
import networkx as nx
Expand Down Expand Up @@ -109,8 +109,11 @@ def calc_sigma(G, dirpar_edges, key="name", output_strings=True):
----------
G : ``NetworkX.MultiDiGraph``
A kinetic diagram
dirpar_edges : list
List of all directional diagrams for the input diagram ``G``.
dirpar_edges : ndarray
Array of all directional diagram edges (made from 2-tuples)
for the input diagram ``G``. Created using
:meth:`~kda.diagrams.generate_directional_diagrams`
with ``return_edges=True``.
key : str
Attribute key used to retrieve edge data from ``G.edges``. The default
``NetworkX`` edge key is ``"weight"``, however the ``kda`` edge keys
Expand Down Expand Up @@ -486,7 +489,7 @@ def calc_thermo_force(G, cycle, order, key="name", output_strings=True):
return parsed_thermo_force_str


def calc_state_probs(G, key="name", output_strings=True):
def calc_state_probs(G, key="name", output_strings=True, dir_edges=None):
r"""Generates the state probability expressions using the diagram
method developed by King and Altman :footcite:`king_schematic_1956` and
Hill :footcite:`hill_studies_1966`.
Expand All @@ -506,14 +509,20 @@ def calc_state_probs(G, key="name", output_strings=True):
probabilities using numbers. If ``True``, this will assume the input
``'key'`` will return strings of variable names to join into the
analytic state multplicity and normalization function.
dir_edges : ndarray (optional)
Array of all directional diagram edges (made from 2-tuples)
for the input diagram ``G``. Created using
:meth:`~kda.diagrams.generate_directional_diagrams`
with ``return_edges=True``.
Returns
-------
state_probs : ndarray
state_probs : ndarray or list of ``SymPy`` expressions
Array of state probabilities for ``N`` states
of the form ``[p1, p2, p3, ..., pN]``.
state_probs_sympy : ``SymPy`` expression
List of symbolic state probability expressions.
of the form ``[p1, p2, p3, ..., pN]``
(for ``output_strings=False``), or a
list of symbolic state probability expressions
in the same order (for ``output_strings=True``).
Notes
-----
Expand All @@ -538,16 +547,62 @@ def calc_state_probs(G, key="name", output_strings=True):
all :math:`\Omega_i` s) for the kinetic diagram.
"""
dirpar_edges = diagrams.generate_directional_diagrams(G, return_edges=True)
state_probs = calc_state_probs_from_diags(
G, dirpar_edges=dirpar_edges, key=key, output_strings=output_strings,
)
edge_is_str = isinstance(G.edges[list(G.edges)[0]][key], str)
if output_strings != edge_is_str:
msg = f"""Inputs `key={key}` and `output_strings={output_strings}`
do not match. If symbolic outputs are requested the input `key`
should retrieve edge data from `G` that corresponds to symbolic
variable names for all edges."""
raise TypeError(msg)

if dir_edges is None:
# generate the directional diagram edges
dir_edges = diagrams.generate_directional_diagrams(G=G, return_edges=True)
# get the number of nodes/states
n_states = G.number_of_nodes()
# get the number of directional diagrams
n_dirpars = dir_edges.shape[0]
# get the number of partial diagrams
n_partials = int(n_dirpars / n_states)
if output_strings:
state_probs = expressions.construct_sympy_prob_funcs(state_mult_funcs=state_probs)
dirpar_rate_products = np.empty(shape=(n_dirpars,), dtype=object)
for i, edge_list in enumerate(dir_edges):
rate_product_vals = []
for edge in edge_list:
rate_product_vals.append(G.edges[edge][key])
dirpar_rate_products[i] = "*".join(rate_product_vals)

state_mults = np.empty(shape=(n_states,), dtype=object)
dirpar_rate_products = dirpar_rate_products.reshape(n_states, n_partials)
for i, arr in enumerate(dirpar_rate_products):
state_mults[i] = "+".join(arr)
state_probs = expressions.construct_sympy_prob_funcs(state_mult_funcs=state_mults)
else:
# retrieve the rate matrix from G
Kij = graph_utils.retrieve_rate_matrix(G)
# create array of ones for storing rate products
dirpar_rate_products = np.ones(n_dirpars, dtype=float)
# iterate over the directional diagrams
for i, edge_list in enumerate(dir_edges):
# for each edge list, retrieve an array of the ith and jth indices,
# retrieve the values associated with each (i, j) pair, and
# calculate the product of those values
Ki = edge_list[:, 0]
Kj = edge_list[:, 1]
dirpar_rate_products[i] = np.prod(Kij[Ki, Kj])

state_mults = dirpar_rate_products.reshape(n_states, n_partials).sum(axis=1)
state_probs = state_mults / math.fsum(dirpar_rate_products)
if any(elem < 0 for elem in state_probs):
raise ValueError(
"Calculated negative state probabilities, overflow or underflow occurred."
)

return state_probs


def calc_net_cycle_flux(G, cycle, order, key="name", output_strings=True):
def calc_net_cycle_flux(G, cycle, order, key="name",
output_strings=True, dir_edges=None):
r"""Generates the expression for the net cycle flux for
some ``cycle`` in kinetic diagram ``G``.
Expand All @@ -574,11 +629,18 @@ def calc_net_cycle_flux(G, cycle, order, key="name", output_strings=True):
using numbers. If ``True``, this will assume the input ``'key'``
will return strings of variable names to join into the analytic
cycle flux function.
dir_edges : ndarray (optional)
Array of all directional diagram edges (made from 2-tuples)
for the input diagram ``G``. Given as an option for performance reasons
(when calculating net cycle fluxes for multiple cycles it is best to
generate the directional diagram edges up front and provide them).
Created using :meth:`~kda.diagrams.generate_directional_diagrams`
with ``return_edges=True``.
Returns
-------
net_cycle_flux : float or ``SymPy`` expression
Net cycle flux for input cycle.
Net cycle flux for the input ``cycle``.
Notes
-----
Expand All @@ -596,27 +658,44 @@ def calc_net_cycle_flux(G, cycle, order, key="name", output_strings=True):
(i.e. positive) direction is counter-clockwise (CCW).
"""
# generate the directional diagram edges
dir_edges = diagrams.generate_directional_diagrams(G=G, return_edges=True)
net_cycle_flux = calc_net_cycle_flux_from_diags(
G=G, dirpar_edges=dir_edges, cycle=cycle,
order=order, key=key, output_strings=output_strings,
)
if dir_edges is None:
# generate the directional diagram edges
dir_edges = diagrams.generate_directional_diagrams(G=G, return_edges=True)
# generate the flux diagrams
flux_diags = diagrams.generate_flux_diagrams(G=G, cycle=cycle)
# construct the expressions for (Pi+ - Pi-), sigma, and sigma_k
# from the directional diagram edges
pi_diff = calc_pi_difference(
G=G, cycle=cycle, order=order, key=key, output_strings=output_strings)
sigma_K = calc_sigma_K(
G=G, cycle=cycle, flux_diags=flux_diags,
key=key, output_strings=output_strings)
sigma = calc_sigma(
G=G, dirpar_edges=dir_edges, key=key, output_strings=output_strings)
if output_strings:
net_cycle_flux = expressions.construct_sympy_net_cycle_flux_func(
pi_diff_str=pi_diff, sigma_K_str=sigma_K, sigma_str=sigma)
else:
net_cycle_flux = pi_diff * sigma_K / sigma
return net_cycle_flux


def calc_state_probs_from_diags(G, dirpar_edges, key="name", output_strings=True):
"""Generates the state probability expressions using the diagram
method developed by King and Altman :footcite:`king_schematic_1956` and
Hill :footcite:`hill_studies_1966`. If directional diagram edges are already
generated this offers better performance than :func:`calc_state_probs`.
generated this offers better performance than
:meth:`~kda.calculations.calc_state_probs`.
Parameters
----------
G : ``NetworkX.MultiDiGraph``
A kinetic diagram
dirpar_edges : array
Array of all directional diagram edges (made from 2-tuples).
Array of all directional diagram edges (made from 2-tuples)
for the input diagram ``G``. Created using
:meth:`~kda.diagrams.generate_directional_diagrams`
with ``return_edges=True``.
key : str
Attribute key used to retrieve edge data from ``G.edges``. The default
``NetworkX`` edge key is ``"weight"``, however the ``kda`` edge keys
Expand All @@ -631,78 +710,42 @@ def calc_state_probs_from_diags(G, dirpar_edges, key="name", output_strings=True
Returns
-------
state_probabilities : ndarray
Array of state probabilities for ``N`` states of
the form ``[p1, p2, p3, ..., pN]``.
state_mults : list of str
List of algebraic state multiplicity expressions.
state_probs : ndarray or list of ``SymPy`` expressions
Array of state probabilities for ``N`` states
of the form ``[p1, p2, p3, ..., pN]``
(for ``output_strings=False``), or a
list of symbolic state probability expressions
in the same order (for ``output_strings=True``).
"""
# get the number of nodes/states
n_states = G.number_of_nodes()
# get the number of directional diagrams
n_dirpars = dirpar_edges.shape[0]
# get the number of partial diagrams
n_partials = int(n_dirpars / n_states)
# retrieve the rate matrix from G
Kij = graph_utils.retrieve_rate_matrix(G)

edge_value = G.edges[list(G.edges)[0]][key]
if not output_strings:
if isinstance(edge_value, str):
raise TypeError(
"To enter variable strings set parameter output_strings=True."
)
# create array of ones for storing rate products
dirpar_rate_products = np.ones(n_dirpars, dtype=float)
# iterate over the directional diagrams
for i, edge_list in enumerate(dirpar_edges):
# for each edge list, retrieve an array of the ith and jth indices,
# retrieve the values associated with each (i, j) pair, and
# calculate the product of those values
Ki = edge_list[:, 0]
Kj = edge_list[:, 1]
dirpar_rate_products[i] = np.prod(Kij[Ki, Kj])

state_mults = dirpar_rate_products.reshape(n_states, n_partials).sum(axis=1)
state_probs = state_mults / math.fsum(dirpar_rate_products)
if any(elem < 0 for elem in state_probs):
raise ValueError(
"Calculated negative state probabilities, overflow or underflow occurred."
)
return state_probs
elif output_strings:
if not isinstance(edge_value, str):
raise TypeError(
"To enter variable values set parameter output_strings=False."
)
dirpar_rate_products = np.empty(shape=(n_dirpars,), dtype=object)
for i, edge_list in enumerate(dirpar_edges):
rate_product_vals = []
for edge in edge_list:
rate_product_vals.append(G.edges[edge][key])
dirpar_rate_products[i] = "*".join(rate_product_vals)

state_mults = np.empty(shape=(n_states,), dtype=object)
dirpar_rate_products = dirpar_rate_products.reshape(n_states, n_partials)
for i, arr in enumerate(dirpar_rate_products):
state_mults[i] = "+".join(arr)
return state_mults
msg = """`kda.calculations.calc_state_probs_from_diags` will be deprecated.
Use `kda.calculations.calc_state_probs` with parameter `dir_edges`."""
warnings.warn(msg, DeprecationWarning)
state_probs = calc_state_probs(
G=G, dir_edges=dirpar_edges, key=key, output_strings=output_strings,
)
if output_strings:
state_probs = expressions.construct_sympy_prob_funcs(state_mult_funcs=state_probs)
return state_probs


def calc_net_cycle_flux_from_diags(
G, dirpar_edges, cycle, order, key="name", output_strings=True
):
"""Generates the expression for the net cycle flux for some ``cycle``
in kinetic diagram ``G``. If directional diagram edges are already
generated this offers better performance than :func:`calc_net_cycle_flux`.
generated this offers better performance than
:meth:`~kda.calculations.calc_net_cycle_flux`.
Parameters
----------
G : ``NetworkX.MultiDiGraph``
A kinetic diagram
dirpar_edges : ndarray
Array of all directional diagram edges (made from 2-tuples).
Array of all directional diagram edges (made from 2-tuples)
for the input diagram ``G``. Created using
:meth:`~kda.diagrams.generate_directional_diagrams`
with ``return_edges=True``.
cycle : list of int
List of node indices for cycle of interest, index zero. Order of node
indices does not matter.
Expand All @@ -725,22 +768,17 @@ def calc_net_cycle_flux_from_diags(
Returns
-------
net_cycle_flux : float or ``SymPy`` expression
Net cycle flux for input cycle.
Net cycle flux for the input ``cycle``.
"""
flux_diags = diagrams.generate_flux_diagrams(G=G, cycle=cycle)
# construct the expressions for (Pi+ - Pi-), sigma, and sigma_k
# from the directional diagram edges
pi_diff = calc_pi_difference(
G=G, cycle=cycle, order=order, key=key, output_strings=output_strings)
sigma_K = calc_sigma_K(
G=G, cycle=cycle, flux_diags=flux_diags,
key=key, output_strings=output_strings)
sigma = calc_sigma(
G=G, dirpar_edges=dirpar_edges, key=key, output_strings=output_strings)
if output_strings:
net_cycle_flux = expressions.construct_sympy_net_cycle_flux_func(
pi_diff_str=pi_diff, sigma_K_str=sigma_K, sigma_str=sigma)
else:
net_cycle_flux = pi_diff * sigma_K / sigma
return net_cycle_flux
msg = """`kda.calculations.calc_net_cycle_flux_from_diags` will be deprecated.
Use `kda.calculations.calc_net_cycle_flux` with parameter `dir_edges`."""
warnings.warn(msg, DeprecationWarning)
return calc_net_cycle_flux(
G=G,
cycle=cycle,
order=order,
key=key,
output_strings=output_strings,
dir_edges=dirpar_edges,
)
Loading

0 comments on commit 2a7cc47

Please sign in to comment.