本文整理汇总了Python中msmbuilder.msm_analysis.check_transition函数的典型用法代码示例。如果您正苦于以下问题：Python check_transition函数的具体用法？Python check_transition怎么用？Python check_transition使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了check_transition函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: calculate_all_to_all_mfpt
def calculate_all_to_all_mfpt(tprob, populations=None):
"""
Calculate the all-states by all-state matrix of mean first passage
times.
This uses the fundamental matrix formalism, and should be much faster
than GetMFPT for calculating many MFPTs.
Parameters
----------
tprob : matrix
transition probability matrix
populations : array_like, float
optional argument, the populations of each state. If not supplied,
it will be computed from scratch
Returns
-------
MFPT : array, float
MFPT in time units of LagTime, square array for MFPT from i -> j
See Also
--------
GetMFPT : function
for calculating a subset of the MFPTs, with functionality for including
a set of sinks
"""
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
tprob = tprob.toarray()
logger.warning('calculate_all_to_all_mfpt does not support sparse linear algebra')
if populations is None:
eigens = msm_analysis.get_eigenvectors(tprob, 5)
if np.count_nonzero(np.imag(eigens[1][:,0])) != 0:
raise ValueError('First eigenvector has imaginary parts')
populations = np.real(eigens[1][:,0])
# ensure that tprob is a transition matrix
msm_analysis.check_transition(tprob)
num_states = len(populations)
if tprob.shape[0] != num_states:
raise ValueError("Shape of tprob and populations vector don't match")
eye = np.transpose( np.matrix(np.ones(num_states)) )
limiting_matrix = eye * populations
#z = scipy.linalg.inv(scipy.sparse.eye(num_states, num_states) - (tprob - limiting_matrix))
z = scipy.linalg.inv(np.eye(num_states) - (tprob - limiting_matrix))
# mfpt[i,j] = z[j,j] - z[i,j] / pi[j]
mfpt = -z
for j in range(num_states):
mfpt[:, j] += z[j, j]
mfpt[:, j] /= populations[j]
return mfpt示例2: calculate_mfpt
def calculate_mfpt(sinks, tprob, lag_time=1.):
"""
Gets the Mean First Passage Time (MFPT) for all states to a *set*
of sinks.
Parameters
----------
sinks : array, int
indices of the sink states
tprob : matrix
transition probability matrix
LagTime : float
the lag time used to create T (dictates units of the answer)
Returns
-------
MFPT : array, float
MFPT in time units of LagTime, for each state (in order of state index)
See Also
--------
calculate_all_to_all_mfpt : function
A more efficient way to calculate all the MFPTs in a network
"""
sinks = _ensure_iterable(sinks)
msm_analysis.check_transition(tprob)
n = tprob.shape[0]
if scipy.sparse.isspmatrix(tprob):
tprob = tprob.tolil()
for state in sinks:
tprob[state,:] = 0.0
tprob[state,state] = 2.0
if scipy.sparse.isspmatrix(tprob):
tprob = tprob - scipy.sparse.eye(n,n)
tprob = tprob.tocsr()
else:
tprob = tprob - np.eye(n)
RHS = -1 * np.ones(n)
for state in sinks:
RHS[state] = 0.0
if scipy.sparse.isspmatrix(tprob):
MFPT = lag_time * scipy.sparse.linalg.spsolve(tprob, RHS)
else:
MFPT = lag_time * np.linalg.solve(tprob, RHS)
return MFPT示例3: calculate_net_fluxes
def calculate_net_fluxes(sources, sinks, tprob, populations=None, committors=None):
"""
Computes the transition path theory net flux matrix.
Parameters
----------
sources : array_like, int
The set of unfolded/reactant states.
sinks : array_like, int
The set of folded/product states.
tprob : mm_matrix
The transition matrix.
Returns
------
net_fluxes : mm_matrix
The net flux matrix
Optional Parameters
-------------------
populations : nd_array, float
The equilibrium populations, if not provided is re-calculated
committors : nd_array, float
The committors associated with `sources`, `sinks`, and `tprob`.
If not provided, is calculated from scratch. If provided, `sources`
and `sinks` are ignored.
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
dense = False
else:
dense = True
n = tprob.shape[0]
flux = calculate_fluxes(sources, sinks, tprob, populations, committors)
ind = flux.nonzero()
if dense:
net_flux = np.zeros((n, n))
else:
net_flux = scipy.sparse.lil_matrix((n, n))
for k in range(len(ind[0])):
i, j = ind[0][k], ind[1][k]
forward = flux[i, j]
reverse = flux[j, i]
net_flux[i, j] = max(0, forward - reverse)
return net_flux示例4: calculate_ensemble_mfpt
def calculate_ensemble_mfpt(sources, sinks, tprob, lag_time):
"""
Calculates the average 'Folding Time' of an MSM defined by T and a LagTime.
The Folding Time is the average of the MFPTs (to F) of all the states in U.
Note here 'Folding Time' is defined as the avg MFPT of {U}, to {F}.
Consider this carefully. This is probably NOT the experimental folding time!
Parameters
----------
sources : array, int
indices of the source states
sinks : array, int
indices of the sink states
tprob : matrix
transition probability matrix
lag_time : float
the lag time used to create T (dictates units of the answer)
Returns
-------
avg : float
the average of the MFPTs
std : float
the standard deviation of the MFPTs
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
X = calculate_mfpt(sinks, tprob, lag_time)
times = np.zeros(len(sources))
for i in range(len(sources)):
times[i] = X[sources[i]]
return np.average(times), np.std(times)示例5: calculate_ensemble_mfpt
def calculate_ensemble_mfpt(sources, sinks, tprob, lag_time):
"""
Calculates the average 'Folding Time' of an MSM defined by T and a LagTime.
The Folding Time is the average of the MFPTs (to F) of all the states in U.
Note here 'Folding Time' is defined as the avg MFPT of {U}, to {F}.
Consider this carefully. This is probably NOT the experimental folding time!
Parameters
----------
sources : array, int
indices of the source states
sinks : array, int
indices of the sink states
tprob : matrix
transition probability matrix
lag_time : float
the lag time used to create T (dictates units of the answer)
Returns
-------
avg : float
the average of the MFPTs
std : float
the standard deviation of the MFPTs
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
X = calculate_mfpt(sinks, tprob, lag_time)
times = np.zeros(len(sources))
for i in range(len(sources)):
times[i] = X[ sources[i] ]
return np.average(times), np.std(times)示例6: find_top_paths
def find_top_paths(sources, sinks, tprob, num_paths=10, node_wipe=False, net_flux=None):
r"""
Calls the Dijkstra algorithm to find the top 'NumPaths'.
Does this recursively by first finding the top flux path, then cutting that
path and relaxing to find the second top path. Continues until NumPaths
have been found.
Parameters
----------
sources : array_like, int
The indices of the source states
sinks : array_like, int
Indices of sink states
num_paths : int
The number of paths to find
Returns
-------
Paths : list of lists
The nodes transversed in each path
Bottlenecks : list of tuples
The nodes between which exists the path bottleneck
Fluxes : list of floats
The flux through each path
Optional Parameters
-------------------
node_wipe : bool
If true, removes the bottleneck-generating node from the graph, instead
of just the bottleneck (not recommended, a debugging functionality)
net_flux : sparse matrix
Matrix of the net flux from `sources` to `sinks`, see function `net_flux`.
If not provided, is calculated from scratch. If provided, `tprob` is
ignored.
To Do
-----
-- Add periodic flow check
References
----------
.. [1] Dijkstra, E. W. (1959). "A note on two problems in connexion with
graphs". Numerische Mathematik 1: 269–271. doi:10.1007/BF01386390.
"""
# first, do some checking on the input, esp. `sources` and `sinks`
# we want to make sure all objects are iterable and the sets are disjoint
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
# check to see if we get net_flux for free, otherwise calculate it
if not net_flux:
net_flux = calculate_net_fluxes(sources, sinks, tprob)
# initialize objects
paths = []
fluxes = []
bottlenecks = []
if scipy.sparse.issparse(net_flux):
net_flux = net_flux.tolil()
# run the initial Dijkstra pass
pi, b = Dijkstra(sources, sinks, net_flux)
logger.info("Path Num | Path | Bottleneck | Flux")
i = 1
done = False
while not done:
# First find the highest flux pathway
(path, (b1, b2), flux) = _backtrack(sinks, b, pi, net_flux)
# Add each result to a Paths, Bottlenecks, Fluxes list
if flux == 0:
logger.info("Only %d possible pathways found. Stopping backtrack.", i)
break
paths.append(path)
bottlenecks.append((b1, b2))
fluxes.append(flux)
logger.info("%s | %s | %s | %s ", i, path, (b1, b2), flux)
# Cut the bottleneck, start relaxing from B side of the cut
if node_wipe:
net_flux[:, b2] = 0
logger.info("Wiped node: %s", b2)
else:
net_flux[b1, b2] = 0
G = scipy.sparse.find(net_flux)
Q = [b2]
b, pi, net_flux = _back_relax(b2, b, pi, net_flux)
# Then relax the graph and repeat
# But only if we still need to
if i != num_paths - 1:
while len(Q) > 0:
w = Q.pop()
#.........这里部分代码省略.........
示例7: calculate_all_to_all_mfpt
def calculate_all_to_all_mfpt(tprob, populations=None):
"""
Calculate the all-states by all-state matrix of mean first passage
times.
This uses the fundamental matrix formalism, and should be much faster
than GetMFPT for calculating many MFPTs.
Parameters
----------
tprob : matrix
transition probability matrix
populations : array_like, float
optional argument, the populations of each state. If not supplied,
it will be computed from scratch
Returns
-------
MFPT : array, float
MFPT in time units of LagTime, square array for MFPT from i -> j
See Also
--------
GetMFPT : function
for calculating a subset of the MFPTs, with functionality for including
a set of sinks
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
tprob = tprob.toarray()
logger.warning('calculate_all_to_all_mfpt does not support sparse linear algebra')
if populations is None:
eigens = msm_analysis.get_eigenvectors(tprob, 1)
if np.count_nonzero(np.imag(eigens[1][:, 0])) != 0:
raise ValueError('First eigenvector has imaginary parts')
populations = np.real(eigens[1][:, 0])
# ensure that tprob is a transition matrix
msm_analysis.check_transition(tprob)
num_states = len(populations)
if tprob.shape[0] != num_states:
raise ValueError("Shape of tprob and populations vector don't match")
eye = np.transpose(np.matrix(np.ones(num_states)))
limiting_matrix = eye * populations
#z = scipy.linalg.inv(scipy.sparse.eye(num_states, num_states) - (tprob - limiting_matrix))
z = scipy.linalg.inv(np.eye(num_states) - (tprob - limiting_matrix))
# mfpt[i,j] = z[j,j] - z[i,j] / pi[j]
mfpt = -z
for j in range(num_states):
mfpt[:, j] += z[j, j]
mfpt[:, j] /= populations[j]
return mfpt示例8: calculate_mfpt
def calculate_mfpt(sinks, tprob, lag_time=1.):
"""
Gets the Mean First Passage Time (MFPT) for all states to a *set*
of sinks.
Parameters
----------
sinks : array, int
indices of the sink states
tprob : matrix
transition probability matrix
LagTime : float
the lag time used to create T (dictates units of the answer)
Returns
-------
MFPT : array, float
MFPT in time units of LagTime, for each state (in order of state index)
See Also
--------
calculate_all_to_all_mfpt : function
A more efficient way to calculate all the MFPTs in a network
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
sinks = _ensure_iterable(sinks)
msm_analysis.check_transition(tprob)
n = tprob.shape[0]
if scipy.sparse.isspmatrix(tprob):
tprob = tprob.tolil()
for state in sinks:
tprob[state, :] = 0.0
tprob[state, state] = 2.0
if scipy.sparse.isspmatrix(tprob):
tprob = tprob - scipy.sparse.eye(n, n)
tprob = tprob.tocsr()
else:
tprob = tprob - np.eye(n)
RHS = -1 * np.ones(n)
for state in sinks:
RHS[state] = 0.0
if scipy.sparse.isspmatrix(tprob):
MFPT = lag_time * scipy.sparse.linalg.spsolve(tprob, RHS)
else:
MFPT = lag_time * np.linalg.solve(tprob, RHS)
return MFPT示例9: calculate_avg_TP_time
def calculate_avg_TP_time(sources, sinks, tprob, lag_time):
"""
Calculates the Average Transition Path Time for MSM with: T, LagTime.
The TPTime is the average of the MFPTs (to F) of all the states
immediately adjacent to U, with the U states effectively deleted.
Note here 'TP Time' is defined as the avg MFPT of all adjacent states to {U},
to {F}, ignoring {U}.
Consider this carefully.
Parameters
----------
sources : array, int
indices of the unfolded states
sinks : array, int
indices of the folded states
tprob : matrix
transition probability matrix
lag_time : float
the lag time used to create T (dictates units of the answer)
Returns
-------
avg : float
the average of the MFPTs
std : float
the standard deviation of the MFPTs
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
n = tprob.shape[0]
if scipy.sparse.issparse(tprob):
T = tprob.tolil()
P = scipy.sparse.lil_matrix((n, n))
else:
P = np.zeros((n, n))
for u in sources:
for i in range(n):
if i not in sources:
P[u, i] = T[u, i]
for u in sources:
T[u, :] = np.zeros(n)
T[:, u] = 0
for i in sources:
N = T[i, :].sum()
T[i, :] = T[i, :] / N
X = calculate_mfpt(sinks, tprob, lag_time)
TP = P * X.T
TPtimes = []
for time in TP:
if time != 0:
TPtimes.append(time)
return np.average(TPtimes), np.std(TPtimes)示例10: calculate_net_fluxes
def calculate_net_fluxes(sources, sinks, tprob, populations=None, committors=None):
"""
Computes the transition path theory net flux matrix.
Parameters
----------
sources : array_like, int
The set of unfolded/reactant states.
sinks : array_like, int
The set of folded/product states.
tprob : mm_matrix
The transition matrix.
Returns
------
net_fluxes : mm_matrix
The net flux matrix
Optional Parameters
-------------------
populations : nd_array, float
The equilibrium populations, if not provided is re-calculated
committors : nd_array, float
The committors associated with `sources`, `sinks`, and `tprob`.
If not provided, is calculated from scratch. If provided, `sources`
and `sinks` are ignored.
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
dense = False
else:
dense = True
n = tprob.shape[0]
flux = calculate_fluxes(sources, sinks, tprob, populations, committors)
ind = flux.nonzero()
if dense:
net_flux = np.zeros((n, n))
else:
net_flux = scipy.sparse.lil_matrix((n, n))
for k in range(len(ind[0])):
i, j = ind[0][k], ind[1][k]
forward = flux[i, j]
reverse = flux[j, i]
net_flux[i, j] = max(0, forward - reverse)
return net_flux示例11: calculate_fluxes
def calculate_fluxes(sources, sinks, tprob, populations=None, committors=None):
"""
Compute the transition path theory flux matrix.
Parameters
----------
sources : array_like, int
The set of unfolded/reactant states.
sinks : array_like, int
The set of folded/product states.
tprob : mm_matrix
The transition matrix.
Returns
------
fluxes : mm_matrix
The flux matrix
Optional Parameters
-------------------
populations : nd_array, float
The equilibrium populations, if not provided is re-calculated
committors : nd_array, float
The committors associated with `sources`, `sinks`, and `tprob`.
If not provided, is calculated from scratch. If provided, `sources`
and `sinks` are ignored.
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
dense = False
else:
dense = True
# check if we got the populations
if populations is None:
eigens = msm_analysis.get_eigenvectors(tprob, 5)
if np.count_nonzero(np.imag(eigens[1][:,0])) != 0:
raise ValueError('First eigenvector has imaginary components')
populations = np.real(eigens[1][:,0])
# check if we got the committors
if committors is None:
committors = calculate_committors(sources, sinks, tprob)
# perform the flux computation
Indx, Indy = tprob.nonzero()
n = tprob.shape[0]
if dense:
X = np.zeros((n, n))
Y = np.zeros((n, n))
X[(np.arange(n), np.arange(n))] = populations * (1.0 - committors)
Y[(np.arange(n), np.arange(n))] = committors
else:
X = scipy.sparse.lil_matrix((n,n))
Y = scipy.sparse.lil_matrix((n,n))
X.setdiag( populations * (1.0 - committors))
Y.setdiag(committors)
if dense:
fluxes = np.dot(np.dot(X, tprob), Y)
fluxes[(np.arange(n), np.arange(n))] = np.zeros(n)
else:
fluxes = np.dot(np.dot(X.tocsr(), tprob.tocsr()), Y.tocsr())
fluxes = fluxes.tolil()
fluxes.setdiag(np.zeros(n))
return fluxes示例12: calculate_hub_score
def calculate_hub_score(tprob, waypoint):
"""
Calculate the hub score for the states `waypoint`.
The "hub score" is a measure of how well traveled a certain state or
set of states is in a network. Specifically, it is the fraction of
times that a walker visits a state en route from some state A to another
state B, averaged over all combinations of A and B.
Parameters
----------
tprob : matrix
The transition probability matrix
waypoints : int
The indices of the intermediate state(s)
Returns
-------
Hc : float
The hub score for the state composed of `waypoints`
See Also
--------
calculate_fraction_visits : function
Calculate the fraction of times a state is visited on pathways going
from a set of "sources" to a set of "sinks".
calculate_all_hub_scores : function
A more efficient way to compute the hub score for every state in a
network.
Notes
-----
Employs dense linear algebra,
memory use scales as N^2
cycle use scales as N^5
References
----------
..[1] Dickson & Brooks (2012), J. Chem. Theory Comput.,
Article ASAP DOI: 10.1021/ct300537s
"""
msm_analysis.check_transition(tprob)
# typecheck
if type(waypoint) != int:
if hasattr(waypoint, '__len__'):
if len(waypoint) == 1:
waypoint = waypoint[0]
else:
raise ValueError('Must pass waypoints as int or list/array of ints')
else:
raise ValueError('Must pass waypoints as int or list/array of ints')
# find out which states to include in A, B (i.e. everything but C)
N = tprob.shape[0]
states_to_include = list(range(N))
states_to_include.remove(waypoint)
# calculate the hub score
Hc = 0.0
for s1 in states_to_include:
for s2 in states_to_include:
if (s1 != s2) and (s1 != waypoint) and (s2 != waypoint):
Hc += calculate_fraction_visits(tprob, waypoint,
s1, s2, return_cond_Q=False)
Hc /= ((N - 1) * (N - 2))
return Hc示例13: calculate_committors
def calculate_committors(sources, sinks, tprob):
"""
Get the forward committors of the reaction sources -> sinks.
Parameters
----------
sources : array_like, int
The set of unfolded/reactant states.
sinks : array_like, int
The set of folded/product states.
tprob : mm_matrix
The transition matrix.
Returns
-------
Q : array_like
The forward committors for the reaction U -> F.
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
dense = False
tprob = tprob.tolil()
else:
dense = True
# construct the committor problem
n = tprob.shape[0]
if dense:
T = np.eye(n) - tprob
else:
T = scipy.sparse.eye(n, n, 0, format='lil') - tprob
T = T.tolil()
for a in sources:
T[a,:] = 0.0 #np.zeros(n)
T[:,a] = 0.0
T[a,a] = 1.0
for b in sinks:
T[b,:] = 0.0 # np.zeros(n)
T[:,b] = 0.0
T[b,b] = 1.0
IdB = np.zeros(n)
IdB[sinks] = 1.0
if dense:
RHS = np.dot(tprob, IdB)
else:
RHS = tprob * IdB
RHS[sources] = 0.0
RHS[sinks] = 1.0
# solve for the committors
if dense == False:
Q = scipy.sparse.linalg.spsolve(T.tocsr(), RHS)
else:
Q = np.linalg.solve(T, RHS)
assert np.all( Q <= 1.0 )
assert np.all( Q >= 0.0 )
return Q示例14: calculate_avg_TP_time
def calculate_avg_TP_time(sources, sinks, tprob, lag_time):
"""
Calculates the Average Transition Path Time for MSM with: T, LagTime.
The TPTime is the average of the MFPTs (to F) of all the states
immediately adjacent to U, with the U states effectively deleted.
Note here 'TP Time' is defined as the avg MFPT of all adjacent states to {U},
to {F}, ignoring {U}.
Consider this carefully.
Parameters
----------
sources : array, int
indices of the unfolded states
sinks : array, int
indices of the folded states
tprob : matrix
transition probability matrix
lag_time : float
the lag time used to create T (dictates units of the answer)
Returns
-------
avg : float
the average of the MFPTs
std : float
the standard deviation of the MFPTs
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
n = tprob.shape[0]
if scipy.sparse.issparse(tprob):
T = tprob.tolil()
P = scipy.sparse.lil_matrix((n, n))
else:
p = np.zeros((n, n))
for u in sources:
for i in range(n):
if i not in sources:
P[u, i] = T[u, i]
for u in sources:
T[u, :] = np.zeros(n)
T[:, u] = 0
for i in sources:
N = T[i, :].sum()
T[i,:] = T[i, :]/N
X = calculate_mfpt(sinks, tprob, lag_time)
TP = P * X.T
TPtimes = []
for time in TP:
if time != 0: TPtimes.append(time)
return np.average(TPtimes), np.std(TPtimes)示例15: calculate_committors
def calculate_committors(sources, sinks, tprob):
"""
Get the forward committors of the reaction sources -> sinks.
Parameters
----------
sources : array_like, int
The set of unfolded/reactant states.
sinks : array_like, int
The set of folded/product states.
tprob : mm_matrix
The transition matrix.
Returns
-------
Q : array_like
The forward committors for the reaction U -> F.
References
----------
.. [1] Metzner, P., Schutte, C. & Vanden-Eijnden, E. Transition path theory
for Markov jump processes. Multiscale Model. Simul. 7, 1192–1219
(2009).
.. [2] Berezhkovskii, A., Hummer, G. & Szabo, A. Reactive flux and folding
pathways in network models of coarse-grained protein dynamics. J.
Chem. Phys. 130, 205102 (2009).
"""
sources, sinks = _check_sources_sinks(sources, sinks)
msm_analysis.check_transition(tprob)
if scipy.sparse.issparse(tprob):
dense = False
tprob = tprob.tolil()
else:
dense = True
# construct the committor problem
n = tprob.shape[0]
if dense:
T = np.eye(n) - tprob
else:
T = scipy.sparse.eye(n, n, 0, format='lil') - tprob
T = T.tolil()
for a in sources:
T[a, :] = 0.0 # np.zeros(n)
T[:, a] = 0.0
T[a, a] = 1.0
for b in sinks:
T[b, :] = 0.0 # np.zeros(n)
T[:, b] = 0.0
T[b, b] = 1.0
IdB = np.zeros(n)
IdB[sinks] = 1.0
if dense:
RHS = np.dot(tprob, IdB)
else:
RHS = tprob.dot(IdB)
# This should be the same as below
#RHS = tprob * IdB
RHS[sources] = 0.0
RHS[sinks] = 1.0
# solve for the committors
if dense == False:
Q = scipy.sparse.linalg.spsolve(T.tocsr(), RHS)
else:
Q = np.linalg.solve(T, RHS)
epsilon = 0.001
assert np.all(Q <= 1.0 + epsilon)
assert np.all(Q >= 0.0 - epsilon)
return Q