# Source code for catkit.gen.route

from . import utils
import numpy as np
import itertools

[docs]def get_reaction_routes(nu, sigma, empty_routes=True, independent_only=False):
"""Returns an array of reaction routes. Returns all full
reaction routes by default.

Parameters
----------
nu : ndarray (n, m)
The stoichiometric matrix of n species and m mechanisms.
sigma : ndarray (m, j)
A linearly independent set of reaction routes.
empty_routes : bool
Return the empty routes along with the full routes.
independent_only : bool
Return only a linearly independent set of full reaction routes.
Can take less time.

Returns
-------
FR : ndarray (m, k)
Enumerated full reaction routes.
ER : ndarray (m, l)
Enumerated empty reaction routes.
"""
p = sigma.shape
m = sigma.shape

isFR = ~np.all(np.dot(sigma, nu) == 0, axis=1)
rsample = np.arange(1, p)

n = len(np.where(isFR))
liRR = np.zeros(sigma.shape, dtype=int)
liRR[:n] = sigma[isFR]

ER = [rr for rr in sigma[~isFR]]
FR = [rr for rr in sigma[isFR]]
route_combinations = itertools.combinations(rsample, r=m - 1)
for c in route_combinations:
S = np.repeat(sigma[:, c][None], sigma.shape, axis=0)
eye = np.eye(sigma.shape)[:, :, None]
R = np.concatenate([S, eye], axis=2)

# Does not convert to correct integers without round
values = np.round(np.linalg.det(R)).astype(int)

# Screen trivial solutions
if np.all(values == 0):
continue

# Normalize first index as positive and reduce by gcd
_gcd = utils.list_gcd(values)
values = (values / _gcd).astype(int)
route = (sigma * values[:, None]).sum(axis=0)
route *= np.sign(route[np.nonzero(route)])

OR = np.dot(route, nu)
if not np.all(OR == 0):
match = False
for rr in FR:
match = np.all(rr == route)
if match:
break

if not match:
FR += [route]

if independent_only:
liRR[n] = route
if np.linalg.matrix_rank(liRR) == n + 1:
n += 1

if n == m:
return liRR
elif empty_routes:
match = False
for rr in ER:
match = np.all(rr == route)
if match:
break

if not match:
ER += [route]

if empty_routes:
return np.array(FR), np.array(ER)

return np.array(FR)

[docs]def get_heppel_sellers(nu, terminal):
"""Returns an array of linearly independent reaction routes
as described by Heppel-Sellers reaction route enumeration.

Parameters
----------
nu : ndarray (n, m)
The stoichiometric matrix of n species and m mechanisms.
terminal : ndarray (j,)
Indices of the m species to be considered as terminal

Returns
-------
sigma : ndarray (m, k)
Linearly independent set of Heppel-Sellers reaction routes.
"""
inter = ~np.in1d(np.arange(nu.shape), terminal)

# Setting up some commonly used parameters
n = nu.shape
m = np.linalg.matrix_rank(nu)
inter = np.where(inter)
selection = inter == inter
nuT = nu.T[inter]

# Reduce the intermetiates to those which are linearly independent
for i, s in enumerate(inter[~selection]):
selection[i + 1] = True

A = nuT[selection]
rank = np.linalg.matrix_rank(A)
if rank == A.shape:
if rank == m - 1:
break
continue

selection[i + 1] = False

inter = inter[selection]

# Now choose m - 1 linearly independent elementary steps
p = nu.shape
reactions = np.arange(p)
selection = reactions == 0

for n in reactions[~selection]:
selection[n] = True
R = nu[selection][:, inter]

rank = np.linalg.matrix_rank(R)
if rank == R.shape:
if rank == m - 1:
break
continue

selection[n] = False

# Find a linearly independent set of RR
sigma = np.zeros((p - m + 1, p), dtype=int)

for i, n in enumerate(reactions[~selection]):
selection[n] = True

_nu = nu[selection][:, inter][None]
S = np.repeat(_nu, len(reactions[selection]), axis=0)
eye = np.eye(len(reactions[selection]))[:, :, None]
R = np.concatenate([S, eye], axis=2)

# Does not convert to correct integers without round
values = np.round(np.linalg.det(R)).astype(int)

# Screen trivial solutions
if np.all(values == 0):
selection[n] = False
continue

# Normalize first index as positive and reduce by gcd
_gcd = utils.list_gcd(values)
values *= np.sign(values[np.nonzero(values)])
values = (values / _gcd).astype(int)

sigma[i, selection] = values
selection[n] = False

return sigma

[docs]def get_response_reactions(epsilon, selection=None, species=False):
"""Returns an array of possible response reaction routes for a given
chemical formula array.

Parameters
----------
epsilon : ndarray (n, m)
Chemical formula array of n elements by m molecular species.
selection : ndarray (j,)
Indices of the m species to be considered as terminal
species : bool
Return the indices of the chemical species used.

Returns
-------
RER : ndarray (m, k)
Possible response reaction routes.
index : ndarray (j, k)
Indices of the k terminal species use to produce the
l response reactions.
"""
s = np.linalg.matrix_rank(epsilon)
RER, index = [], []
if not selection:
selection = np.arange(epsilon.shape)

possible_routes = itertools.combinations(selection, r=s + 1)
for sel in possible_routes:
values = np.zeros(epsilon.shape, dtype=int)

sigma = np.repeat(epsilon[[sel]][None], s + 1, axis=0)
eye = np.eye(s + 1)[:, :, None]
R = np.concatenate([sigma, eye], axis=2)

# Does not convert to correct integers without round
values[[sel]] = np.round(np.linalg.det(R))

# Screen trivial solutions
if np.all(values == 0):
continue

# Normalize first index as positive and reduce by gcd
_gcd = utils.list_gcd(values)
values *= np.sign(values[np.nonzero(values)])
values = (values / _gcd).astype(int)

# Screen the stoichiometric matches
match = False
for v in RER:
match = np.all(v == values)
if match:
break

if not match:
RER += [values]
index += [list(sel)]

RER = np.array(RER)
if species:
index = np.array(index)
return RER, index