pyemma.msm.ui.EstimatedMSM¶
-
class
pyemma.msm.ui.
EstimatedMSM
(dtrajs, lag, reversible=True, sparse=False, connectivity='largest', estimate=True, dt='1 step', **kwargs)¶ Estimates a Markov model from discrete trajectories.
Parameters: - dtrajs (list containing ndarrays(dtype=int) or ndarray(n, dtype=int)) – discrete trajectories, stored as integer ndarrays (arbitrary size) or a single ndarray for only one trajectory.
- lag (int) – lagtime for the MSM estimation in multiples of trajectory steps
- reversible (bool, optional, default = True) – If true compute reversible MSM, else non-reversible MSM
- sparse (bool, optional, default = False) – If true compute count matrix, transition matrix and all derived quantities using sparse matrix algebra. In this case python sparse matrices will be returned by the corresponding functions instead of numpy arrays. This behavior is suggested for very large numbers of states (e.g. > 4000) because it is likely to be much more efficient.
- connectivity (str, optional, default = 'largest') –
Connectivity mode. Three methods are intended (currently only ‘largest’ is implemented) ‘largest’ : The active set is the largest reversibly connected set. All estimation will be done on this
subset and all quantities (transition matrix, stationary distribution, etc) are only defined on this subset and are correspondingly smaller than the full set of states- ‘all’ : The active set is the full set of states. Estimation will be conducted on each reversibly connected
- set separately. That means the transition matrix will decompose into disconnected submatrices, the stationary vector is only defined within subsets, etc. Currently not implemented.
- ‘none’ : The active set is the full set of states. Estimation will be conducted on the full set of states
- without ensuring connectivity. This only permits nonreversible estimation. Currently not implemented.
- estimate (bool, optional, default=True) – If true estimate the MSM when creating the MSM object.
- dt (str, optional, default='1 step') –
Description of the physical time corresponding to the lag. May be used by analysis algorithms such as plotting tools to pretty-print the axes. By default ‘1 step’, i.e. there is no physical time unit. Specify by a number, whitespace and unit. Permitted units are (* is an arbitrary string):
‘fs’, ‘femtosecond*’‘ps’, ‘picosecond*’‘ns’, ‘nanosecond*’‘us’, ‘microsecond*’‘ms’, ‘millisecond*’‘s’, ‘second*’ - **kwargs –
- = 1000000 (maxiter) – Optional parameter with reversible = True. maximum number of iterations before the transition matrix estimation method exits
- = 1e-8 (maxerr) – Optional parameter with reversible = True. convergence tolerance for transition matrix estimation. This specifies the maximum change of the Euclidean norm of relative stationary probabilities (\(x_i = \sum_k x_{ik}\)). The relative stationary probability changes \(e_i = (x_i^{(1)} - x_i^{(2)})/(x_i^{(1)} + x_i^{(2)})\) are used in order to track changes in small probabilities. The Euclidean norm of the change vector, \(|e_i|_2\), is compared to maxerr.
Notes
You can postpone the estimation of the MSM using estimate=False and initiate the estimation procedure by manually calling the MSM.estimate() method.
-
__init__
(dtrajs, lag, reversible=True, sparse=False, connectivity='largest', estimate=True, dt='1 step', **kwargs)¶
Methods
__init__
(dtrajs, lag[, reversible, sparse, ...])committor_backward
(A, B)Backward committor from set A to set B committor_forward
(A, B)Forward committor (also known as p_fold or splitting probability) from set A to set B correlation
(a[, b, maxtime, k, ncv])Time-correlation for equilibrium experiment. eigenvalues
([k, ncv])Compute the transition matrix eigenvalues eigenvectors_left
([k, ncv])Compute the left transition matrix eigenvectors eigenvectors_right
([k, ncv])Compute the right transition matrix eigenvectors estimate
()Runs msm estimation. expectation
(a)Equilibrium expectation value of a given observable. fingerprint_correlation
(a[, b, k, ncv])Dynamical fingerprint for equilibrium time-correlation experiment. fingerprint_relaxation
(p0, a[, k, ncv])Dynamical fingerprint for perturbation/relaxation experiment. generate_traj
(N[, start, stop, stride])Generates a synthetic discrete trajectory of length N and simulation time stride * lag time * N mfpt
(A, B)Mean first passage times from set A to set B, in units of the input trajectory time step pcca
(m)Runs PCCA++ [1]_ in order to compute a fuzzy metastable decomposition of MSM states relaxation
(p0, a[, maxtime, k, ncv])Simulates a perturbation-relaxation experiment. sample_by_distributions
(distributions, nsample)Generates samples according to given probability distributions sample_by_state
(nsample[, subset, replace])Generates samples of the connected states. timescales
([k, ncv])The relaxation timescales corresponding to the eigenvalues trajectory_weights
()Uses the MSM to assign a probability weight to each trajectory frame. Attributes
active_count_fraction
The fraction of counts in the active set. active_set
The active set of states on which all computations and estimations will be done active_state_fraction
The fraction of states in the active set. active_state_indexes
Ensures that the connected states are indexed and returns the indices computed
Returns whether this msm has been estimated yet connected_sets
The reversible connected sets of states, sorted by size (descending) count_matrix_active
The count matrix on the active set given the connectivity mode used. count_matrix_full
The count matrix on full set of discrete states, irrespective as to whether they are connected or not. discrete_trajectories_active
A list of integer arrays with the discrete trajectories mapped to the connectivity mode used. discrete_trajectories_full
A list of integer arrays with the original (unmapped) discrete trajectories effective_count_matrix
Statistically uncorrelated transition counts within the active set of states is_reversible
Returns whether the MSM is reversible is_sparse
Returns whether the MSM is sparse lagtime
The lag time at which the Markov model was estimated largest_connected_set
The largest reversible connected set of states metastable_assignments
Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_ metastable_distributions
Computes the probability distributions of active set states within each metastable set using the PCCA++ method [1]_ using Bayesian inversion as described in [2]_. metastable_memberships
Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_. metastable_sets
Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_ nstates
The active set of states on which all computations and estimations will be done stationary_distribution
The stationary distribution, estimated on the active set. timestep
Returns the physical time corresponding to one step of the transition matrix as string, e.g. transition_matrix
The transition matrix, estimated on the active set. -
active_count_fraction
¶ The fraction of counts in the active set.
-
active_set
¶ The active set of states on which all computations and estimations will be done
-
active_state_fraction
¶ The fraction of states in the active set.
-
active_state_indexes
¶ Ensures that the connected states are indexed and returns the indices
-
committor_backward
(A, B)¶ Backward committor from set A to set B
Parameters: - A (int or int array) – set of starting states
- B (int or int array) – set of target states
-
committor_forward
(A, B)¶ Forward committor (also known as p_fold or splitting probability) from set A to set B
Parameters: - A (int or int array) – set of starting states
- B (int or int array) – set of target states
-
computed
¶ Returns whether this msm has been estimated yet
-
connected_sets
¶ The reversible connected sets of states, sorted by size (descending)
-
correlation
(a, b=None, maxtime=None, k=None, ncv=None)¶ Time-correlation for equilibrium experiment.
In order to simulate a time-correlation experiment (e.g. fluorescence correlation spectroscopy [1]_, dynamical neutron scattering [2]_, ...), first compute the mean values of your experimental observable \(a\) by MSM state:
where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\). If a cross-correlation function is wanted, also apply the above computation to a second experimental observable \(b\).
Then the precise (i.e. without statistical error) autocorrelation function of \(f(x_t)\) given the Markov model is computed by correlation(a), and the precise cross-correlation function is computed by correlation(a,b). This is done by evaluating the equation
where \(acf\) stands for autocorrelation function and \(ccf\) stands for cross-correlation function, \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.
Note that instead of using this method you could generate a long synthetic trajectory from the MSM using
generate_traj()
and then estimating the time-correlation of your observable(s) directly from this trajectory. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely long synthetic trajectory the two results will agree exactly. The correlation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.Parameters: - a ((M,) ndarray) – Observable, represented as vector on state space
- maxtime (int or float) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\) By default (None), the maxtime will be set equal to the 3 times the slowest relaxation time of the MSM, because after this time the signal is constant.
- b ((M,) ndarray (optional)) – Second observable, for cross-correlations
- k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power.
- ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: - times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the correlation has been computed
- correlations (ndarray (N)) – Correlation values at given times
Examples
This example computes the autocorrelation function of a simple observable on a three-state Markov model and plots the result using matplotlib:
>>> import numpy as np >>> import pyemma.msm as msm >>> >>> P = np.array([[0.99, 0.01, 0], [0.01, 0.9, 0.09], [0, 0.1, 0.9]]) >>> a = np.array([0.0, 0.5, 1.0]) >>> M = msm.markov_model(P) >>> times, acf = M.correlation(a) >>> >>> import matplotlib.pylab as plt >>> plt.plot(times, acf)
References
[1] Noe, F., S. Doose, I. Daidone, M. Loellmann, J. D. Chodera, M. Sauer and J. C. Smith. 2011 Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. Proc. Natl. Acad. Sci. USA 108, 4822-4827. [2] Lindner, B., Z. Yi, J.-H. Prinz, J. C. Smith and F. Noe. 2013. Dynamic Neutron Scattering from Conformational Dynamics I: Theory and Markov models. J. Chem. Phys. 139, 175101.
-
count_matrix_active
¶ The count matrix on the active set given the connectivity mode used.
For example, for connectivity=’largest’, the count matrix is given only on the largest reversibly connected set. Attention: This count matrix has been obtained by sliding a window of length tau across the data. It contains a factor of tau more counts than are statistically uncorrelated. It’s fine to use this matrix for maximum likelihood estimated, but it will give far too small errors if you use it for uncertainty calculations. In order to do uncertainty calculations, use the effective count matrix, see:
effective_count_matrix()
See also
effective_count_matrix
- For a count matrix with effective (statistically uncorrelated) counts.
-
count_matrix_full
¶ The count matrix on full set of discrete states, irrespective as to whether they are connected or not. Attention: This count matrix has been obtained by sliding a window of length tau across the data. It contains a factor of tau more counts than are statistically uncorrelated. It’s fine to use this matrix for maximum likelihood estimated, but it will give far too small errors if you use it for uncertainty calculations. In order to do uncertainty calculations, use the effective count matrix, see: :attribute:`effective_count_matrix` (only implemented on the active set), or divide this count matrix by tau.
See also
effective_count_matrix
- For a active-set count matrix with effective (statistically uncorrelated) counts.
-
discrete_trajectories_active
¶ A list of integer arrays with the discrete trajectories mapped to the connectivity mode used. For example, for connectivity=’largest’, the indexes will be given within the connected set. Frames that are not in the connected set will be -1.
-
discrete_trajectories_full
¶ A list of integer arrays with the original (unmapped) discrete trajectories
-
effective_count_matrix
¶ Statistically uncorrelated transition counts within the active set of states
You can use this count matrix for any kind of estimation, in particular it is mean to give reasonable error bars in uncertainty measurements (error perturbation or Gibbs sampling of the posterior).
The effective count matrix is obtained by dividing the sliding-window count matrix by the lag time. This can be shown to provide a likelihood that is the geometrical average over shifted subsamples of the trajectory, \((s_1,\:s_{tau+1},\:...),\:(s_2,\:t_{tau+2},\:...),\) etc. This geometrical average converges to the correct likelihood in the statistical limit _[1].
[1] Trendelkamp-Schroer B, H Wu, F Paul and F Noe. 2015: Reversible Markov models of molecular kinetics: Estimation and uncertainty. in preparation.
-
eigenvalues
(k=None, ncv=None)¶ Compute the transition matrix eigenvalues
Parameters: - k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1
- ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: ts – transition matrix eigenvalues \(\lambda_i, i = 1,...,k\)., sorted by descending norm.
Return type: ndarray(m)
-
eigenvectors_left
(k=None, ncv=None)¶ Compute the left transition matrix eigenvectors
Parameters: - k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1
- ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: L – left eigenvectors in a row matrix. l_ij is the j’th component of the i’th left eigenvector
Return type: ndarray(k,n)
-
eigenvectors_right
(k=None, ncv=None)¶ Compute the right transition matrix eigenvectors
Parameters: - k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1
- ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: R – right eigenvectors in a column matrix. r_ij is the i’th component of the j’th right eigenvector
Return type: ndarray(n,k)
-
estimate
()¶ Runs msm estimation.
Only need to call this method if the msm was initialized with compute=False - otherwise it will have been called at time of initialization.
-
expectation
(a)¶ Equilibrium expectation value of a given observable.
Parameters: a ((M,) ndarray) – Observable vector Returns: val – Equilibrium expectation value fo the given observable Return type: float Notes
The equilibrium expectation value of an observable a is defined as follows
\[\mathbb{E}_{\mu}[a] = \sum_i \mu_i a_i\]\(\mu=(\mu_i)\) is the stationary vector of the transition matrix \(T\).
-
fingerprint_correlation
(a, b=None, k=None, ncv=None)¶ Dynamical fingerprint for equilibrium time-correlation experiment.
Parameters: - a ((M,) ndarray) – Observable, represented as vector on state space
- b ((M,) ndarray (optional)) – Second observable, for cross-correlations
- k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power
- ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: - timescales ((N,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
- amplitudes ((N,) ndarray) – Amplitudes for the correlation experiment
References
[1] Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.
-
fingerprint_relaxation
(p0, a, k=None, ncv=None)¶ Dynamical fingerprint for perturbation/relaxation experiment.
Parameters: - p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
- a ((n,) ndarray) – Observable, represented as vector on state space
- lag (int or int array) – List of lag time or lag times (in units of the transition matrix lag time tau) at which to compute correlation
- k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation
- ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: - timescales ((N,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
- amplitudes ((N,) ndarray) – Amplitudes for the relaxation experiment
References
[1] Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.
-
generate_traj
(N, start=None, stop=None, stride=1)¶ Generates a synthetic discrete trajectory of length N and simulation time stride * lag time * N
This information can be used in order to generate a synthetic molecular dynamics trajectory - see
pyemma.coordinates.save_traj()
Note that the time different between two samples is the Markov model lag time :math:` au`. When comparing quantities computing from this synthetic trajectory and from the input trajectories, the time points of this trajectory must be scaled by the lag time in order to have them on the same time scale.
Parameters: - N (int) – Number of time steps in the output trajectory. The total simulation time is stride * lag time * N
- start (int, optional, default = None) – starting state. If not given, will sample from the stationary distribution of P
- stop (int or int-array-like, optional, default = None) – stopping set. If given, the trajectory will be stopped before N steps once a state of the stop set is reached
- stride (int, optional, default = 1) – Multiple of lag time used as a time step. By default, the time step is equal to the lag time
Returns: indexes – trajectory and time indexes of the simulated trajectory. Each row consist of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory. Note that the time different between two samples is the Markov model lag time :math:` au`.
Return type: ndarray( (N, 2) )
See also
pyemma.coordinates.save_traj()
- in order to save this synthetic trajectory as a trajectory file with molecular structures
-
is_reversible
¶ Returns whether the MSM is reversible
-
is_sparse
¶ Returns whether the MSM is sparse
-
lagtime
¶ The lag time at which the Markov model was estimated
-
largest_connected_set
¶ The largest reversible connected set of states
-
metastable_assignments
¶ Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_
pcca()
needs to be called first before this attribute is available.This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!
Returns: Return type: For each active set state, the metastable state it is located in. See also
pcca
- to compute the metastable decomposition
References
[1] Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
metastable_distributions
¶ Computes the probability distributions of active set states within each metastable set using the PCCA++ method [1]_ using Bayesian inversion as described in [2]_.
pcca()
needs to be called first before this attribute is available.Returns: p_out – A matrix containing the probability distribution of each active set state, given that we are in one of the m metastable sets. i.e. p(state | metastable). The row sums of p_out are 1. Return type: ndarray((m,n)) See also
pcca
- to compute the metastable decomposition
References
[1] Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179 [2] F. Noe, H. Wu, J.-H. Prinz and N. Plattner: Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114 (2013)
-
metastable_memberships
¶ Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_.
pcca()
needs to be called first before this attribute is available.Returns: M – A matrix containing the probability or membership of each state to be assigned to each metastable set. i.e. p(metastable | state). The row sums of M are 1. Return type: ndarray((n,m)) See also
pcca
- to compute the metastable decomposition
References
[1] Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
metastable_sets
¶ Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_
pcca()
needs to be called first before this attribute is available.This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!
Returns: Return type: A list of length equal to metastable states. Each element is an array with microstate indexes contained in it See also
pcca
- to compute the metastable decomposition
References
[1] Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
mfpt
(A, B)¶ Mean first passage times from set A to set B, in units of the input trajectory time step
Parameters: - A (int or int array) – set of starting states
- B (int or int array) – set of target states
-
nstates
¶ The active set of states on which all computations and estimations will be done
-
pcca
(m)¶ Runs PCCA++ [1]_ in order to compute a fuzzy metastable decomposition of MSM states
After calling this method you can access
metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
Parameters: m (int) – Number of metastable sets Returns: pcca_obj – An object containing all PCCA quantities. However, you can also ingore this return value and instead retrieve the quantities of your interest with the following MSM functions: metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
.Return type: PCCA
References
[1] Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
relaxation
(p0, a, maxtime=None, k=None, ncv=None)¶ Simulates a perturbation-relaxation experiment.
In perturbation-relaxation experiments such as temperature-jump, pH-jump, pressure jump or rapid mixing experiments, an ensemble of molecules is initially prepared in an off-equilibrium distribution and the expectation value of some experimental observable is then followed over time as the ensemble relaxes towards equilibrium.
In order to simulate such an experiment, first determine the distribution of states at which the experiment is started, \(p_0\) and compute the mean values of your experimental observable \(a\) by MSM state:
where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\).
Then the precise (i.e. without statistical error) time-dependent expectation value of \(f(x_t)\) given the Markov model is computed by relaxation(p0, a). This is done by evaluating the equation
where \(E\) stands for the expectation value that relaxes to its equilibrium value that is identical to expectation(a), \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.
Note that instead of using this method you could generate many synthetic trajectory from the MSM using
generate_traj()
that with starting points drawn from the initial distribution and then estimating the time-dependent expectation value by an ensemble average. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely many trajectories the two results will agree exactly. The relaxation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.Parameters: - p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
- a ((n,) ndarray) – Observable, represented as vector on state space
- maxtime (int or float, optional, default = None) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\). By default (None), the maxtime will be set equal to the 3 times the slowest relaxation time of the MSM, because after this time the signal is constant.
- k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation
- ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: - times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the relaxation has been computed
- res (ndarray) – Array of expectation value at given times
-
sample_by_distributions
(distributions, nsample)¶ Generates samples according to given probability distributions
Parameters: - distributions (list or array of ndarray ( (n) )) – m distributions over states. Each distribution must be of length n and must sum up to 1.0
- nsample (int) – Number of samples per distribution. If replace = False, the number of returned samples per state could be smaller if less than nsample indexes are available for a state.
Returns: indexes – List of the sampled indices by distribution. Each element is an index array with a number of rows equal to nsample, with rows consisting of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory.
Return type: length m list of ndarray( (nsample, 2) )
-
sample_by_state
(nsample, subset=None, replace=True)¶ Generates samples of the connected states.
For each state in the active set of states, generates nsample samples with trajectory/time indexes. This information can be used in order to generate a trajectory of length nsample * nconnected using
pyemma.coordinates.save_traj()
or nconnected trajectories of length nsample each usingpyemma.coordinates.save_traj()
Parameters: - N (int) – Number of time steps in the output trajectory. The total simulation time is stride * lag time * N
- nsample (int) – Number of samples per state. If replace = False, the number of returned samples per state could be smaller if less than nsample indexes are available for a state.
- subset (ndarray((n)), optional, default = None) – array of states to be indexed. By default all states in the connected set will be used
- replace (boolean, optional) – Whether the sample is with or without replacement
- start (int, optional, default = None) – starting state. If not given, will sample from the stationary distribution of P
Returns: indexes – list of trajectory/time index arrays with an array for each state. Within each index array, each row consist of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory.
Return type: list of ndarray( (N, 2) )
See also
pyemma.coordinates.save_traj()
- in order to save the sampled frames sequentially in a trajectory file with molecular structures
pyemma.coordinates.save_trajs()
- in order to save the sampled frames in nconnected trajectory files with molecular structures
-
stationary_distribution
¶ The stationary distribution, estimated on the active set.
For example, for connectivity=’largest’ it will be the transition matrix amongst the largest set of reversibly connected states
-
timescales
(k=None, ncv=None)¶ The relaxation timescales corresponding to the eigenvalues
Parameters: - k (int) – number of timescales to be computed. As a result, k+1 eigenvalues will be computed
- ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns: ts – relaxation timescales in units of the input trajectory time step, defined by \(-tau / ln | \lambda_i |, i = 2,...,k+1\).
Return type: ndarray(m)
-
timestep
¶ Returns the physical time corresponding to one step of the transition matrix as string, e.g. ‘10 ps’
-
trajectory_weights
()¶ Uses the MSM to assign a probability weight to each trajectory frame.
This is a powerful function for the calculation of arbitrary observables in the trajectories one has started the analysis with. The stationary probability of the MSM will be used to reweigh all states. Returns a list of weight arrays, one for each trajectory, and with a number of elements equal to trajectory frames. Given \(N\) trajectories of lengths \(T_1\) to \(T_N\), this function returns corresponding weights: .. math:
(w_{1,1}, ..., w_{1,T_1}), (w_{N,1}, ..., w_{N,T_N})
that are normalized to one: .. math:
\sum_{i=1}^N \sum_{t=1}^{T_i} w_{i,t} = 1
Suppose you are interested in computing the expectation value of a function \(a(x)\), where \(x\) are your input configurations. Use this function to compute the weights of all input configurations and obtain the estimated expectation by: .. math:
\langle a
- angle = sum_{i=1}^N sum_{t=1}^{T_i} w_{i,t} a(x_{i,t})
Or if you are interested in computing the time-lagged correlation between functions \(a(x)\) and \(b(x)\) you could do: .. math:
\langle a(t) b(t+ au)
angle_t = sum_{i=1}^N sum_{t=1}^{T_i} w_{i,t} a(x_{i,t}) a(x_{i,t+ au})
The normalized trajectory weights. Given \(N\) trajectories of lengths \(T_1\) to \(T_N\), returns the corresponding weights: .. math:
(w_{1,1}, ..., w_{1,T_1}), (w_{N,1}, ..., w_{N,T_N})
-
transition_matrix
¶ The transition matrix, estimated on the active set. For example, for connectivity=’largest’ it will be the transition matrix amongst the largest set of reversibly connected states