pyemma.msm.HMSM

class pyemma.msm.HMSM(P, pobs, pi=None, dt_model='1 step')

Hidden Markov model on discrete states.

Parameters

hmm (DiscreteHMM) – Hidden Markov Model

__init__(P, pobs, pi=None, dt_model='1 step')
Parameters
  • Pcoarse (ndarray (m,m)) – coarse-grained or hidden transition matrix

  • Pobs (ndarray (m,n)) – observation probability matrix from hidden to observable discrete states

  • dt_model (str, optional, default='1 step') – time step of the model

Methods

__init__(P, pobs[, pi, dt_model])

param Pcoarse

coarse-grained or hidden transition matrix

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])

Compute the transition matrix eigenvalues

eigenvectors_left([k])

Compute the left transition matrix eigenvectors

eigenvectors_right([k])

Compute the right transition matrix eigenvectors

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.

get_model_params([deep])

Get parameters for this model.

load(file_name[, model_name])

Loads a previously saved PyEMMA object from disk.

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]_ to compute a metastable decomposition of MSM states

propagate(p0, k)

Propagates the initial distribution p0 k times

relaxation(p0, a[, maxtime, k, ncv])

Simulates a perturbation-relaxation experiment.

save(file_name[, model_name, overwrite, …])

saves the current state of this object to given file and name.

set_model_params([P, pobs, pi, reversible, …])

param P

coarse-grained or hidden transition matrix

simulate(N[, start, stop, dt])

Generates a realization of the Hidden Markov Model

submodel([states, obs])

Returns a HMM with restricted state space

timescales([k])

The relaxation timescales corresponding to the eigenvalues

transition_matrix_obs([k])

Computes the transition matrix between observed states

update_model_params(**params)

Update given model parameter if they are set to specific values

Attributes

P

The transition matrix on the active set.

dt_model

Description of the physical time corresponding to the lag.

eigenvectors_left_obs

eigenvectors_right_obs

is_reversible

Returns whether the MSM is reversible

is_sparse

Returns whether the MSM is sparse

lifetimes

Lifetimes of states of the hidden transition matrix

metastable_assignments

Computes the assignment to metastable sets for observable states

metastable_distributions

Returns the output probability distributions. Identical to

metastable_memberships

Computes the memberships of observable states to metastable sets by

metastable_sets

Computes the metastable sets of observable states within each

n_metastable

Number of states chosen for PCCA++ computation.

neig

number of eigenvalues to compute.

nstates

Number of active states on which all computations and estimations are done

nstates_obs

observation_probabilities

returns the output probability matrix

pi

The stationary distribution on the MSM states

reversible

Returns whether the MSM is reversible

sparse

Returns whether the MSM is sparse

stationary_distribution

The stationary distribution on the MSM states

stationary_distribution_obs

timestep_model

Physical time corresponding to one transition matrix step, e.g.

transition_matrix

The transition matrix on the active set.

P

The transition matrix on the active set.

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

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 [NDD11], dynamical neutron scattering [LYP13], …), first compute the mean values of your experimental observable \(a\) by MSM state:

\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]

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 accurate (i.e. without statistical error) autocorrelation function of \(f(x_t)\) given the Markov model is computed by correlation(a), and the accurate cross-correlation function is computed by correlation(a,b). This is done by evaluating the equation

\[\begin{split}acf_a(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{a} \\ ccf_{a,b}(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{b}\end{split}\]

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 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 ((n,) 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 5 times the slowest relaxation time of the MSM, because after this time the signal is almost constant.

  • b ((n,) 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 computed 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)  # doctest: +SKIP

References

NDD11

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.

LYP13

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.

dt_model

Description of the physical time corresponding to the lag.

eigenvalues(k=None)

Compute the transition matrix eigenvalues

Parameters

k (int) – number of eigenvalues to be returned. By default will return all available eigenvalues

Returns

ts – transition matrix eigenvalues \(\lambda_i, i = 1, ..., k\)., sorted by descending norm.

Return type

ndarray(k,)

eigenvectors_left(k=None)

Compute the left transition matrix eigenvectors

Parameters

k (int) – number of eigenvectors to be returned. By default all available eigenvectors.

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)

Compute the right transition matrix eigenvectors

Parameters

k (int) – number of eigenvectors to be computed. By default all available eigenvectors.

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)

expectation(a)

Equilibrium expectation value of a given observable.

Parameters

a ((n,) ndarray) – Observable vector on the MSM state space

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 \pi_i a_i\]

\(\pi=(\pi_i)\) is the stationary vector of the transition matrix \(P\).

fingerprint_correlation(a, b=None, k=None, ncv=None)

Dynamical fingerprint for equilibrium time-correlation experiment.

Parameters
  • a ((n,) ndarray) – Observable, represented as vector on MSM state space

  • b ((n,) 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 computed 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 ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix

  • amplitudes ((k,) ndarray) – Amplitudes for the correlation experiment

References

Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].

NDD11

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

  • 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 ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix

  • amplitudes ((k,) ndarray) – Amplitudes for the relaxation experiment

References

Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].

NDD11

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.

get_model_params(deep=True)

Get parameters for this model.

Parameters

deep (boolean, optional) – If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns

params – Parameter names mapped to their values.

Return type

mapping of string to any

is_reversible

Returns whether the MSM is reversible

is_sparse

Returns whether the MSM is sparse

lifetimes

Lifetimes of states of the hidden transition matrix

Returns

l – state lifetimes in units of the input trajectory time step, defined by \(-\tau / ln \mid p_{ii} \mid, i = 1,...,nstates\), where \(p_{ii}\) are the diagonal entries of the hidden transition matrix.

Return type

ndarray(nstates)

classmethod load(file_name, model_name='default')

Loads a previously saved PyEMMA object from disk.

Parameters
  • file_name (str or file like object (has to provide read method)) – The file like object tried to be read for a serialized object.

  • model_name (str, default='default') – if multiple models are contained in the file, these can be accessed by their name. Use pyemma.list_models() to get a representation of all stored models.

Returns

obj

Return type

the de-serialized object

metastable_assignments

Computes the assignment to metastable sets for observable states

Notes

This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!

Returns

For each observable state, the metastable state it is located in.

Return type

ndarray((n) ,dtype=int)

metastable_distributions
Returns the output probability distributions. Identical to

observation_probabilities()

Returns

Pout – output probability matrix from hidden to observable discrete states

Return type

ndarray (m,n)

metastable_memberships
Computes the memberships of observable states to metastable sets by

Bayesian inversion as described in [1]_.

Returns

M – A matrix containing the probability or membership of each observable state to be assigned to each metastable or hidden state. The row sums of M are 1.

Return type

ndarray((n,m))

References

1

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_sets
Computes the metastable sets of observable states within each

metastable set

Notes

This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!

Returns

sets – A list of length equal to metastable states. Each element is an array with observable state indexes contained in it

Return type

list of int-arrays

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

n_metastable

Number of states chosen for PCCA++ computation.

neig

number of eigenvalues to compute.

nstates

Number of active states on which all computations and estimations are done

observation_probabilities

returns the output probability matrix

Returns

Pout – output probability matrix from hidden to observable discrete states

Return type

ndarray (m,n)

pcca(m)

Runs PCCA++ [1]_ to compute a metastable decomposition of MSM states

After calling this method you can access metastable_memberships(), metastable_distributions(), metastable_sets() and metastable_assignments().

Parameters

m (int) – Number of metastable sets

Returns

pcca_obj – An object containing all PCCA quantities. However, you can also ignore this return value and instead retrieve the quantities of your interest with the following MSM functions: metastable_memberships(), metastable_distributions(), metastable_sets() and metastable_assignments().

Return type

PCCA

Notes

If you coarse grain with PCCA++, the order of the obtained memberships might not be preserved. This also applies for metastable_memberships(), metastable_distributions(), metastable_sets(), metastable_assignments()

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

pi

The stationary distribution on the MSM states

propagate(p0, k)

Propagates the initial distribution p0 k times

Computes the product

\[p_k = p_0^T P^k\]

If the lag time of transition matrix \(P\) is \(\tau\), this will provide the probability distribution at time \(k \tau\).

Parameters
  • p0 (ndarray(n)) – Initial distribution. Vector of size of the active set.

  • k (int) – Number of time steps

Returns

pk – Distribution after k steps. Vector of size of the active set.

Return type

ndarray(n)

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:

\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]

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 accurate (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

\[E_a(k\tau) = \mathbf{p_0}^{\top} \mathbf{P(\tau)}^k \mathbf{a}\]

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 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) – 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 5 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

reversible

Returns whether the MSM is reversible

save(file_name, model_name='default', overwrite=False, save_streaming_chain=False)

saves the current state of this object to given file and name.

Parameters
  • file_name (str) – path to desired output file

  • model_name (str, default='default') – creates a group named ‘model_name’ in the given file, which will contain all of the data. If the name already exists, and overwrite is False (default) will raise a RuntimeError.

  • overwrite (bool, default=False) – Should overwrite existing model names?

  • save_streaming_chain (boolean, default=False) – if True, the data_producer(s) of this object will also be saved in the given file.

Examples

>>> import pyemma, numpy as np
>>> from pyemma.util.contexts import named_temporary_file
>>> m = pyemma.msm.MSM(P=np.array([[0.1, 0.9], [0.9, 0.1]]))
>>> with named_temporary_file() as file: # doctest: +SKIP
...    m.save(file, 'simple') # doctest: +SKIP
...    inst_restored = pyemma.load(file, 'simple') # doctest: +SKIP
>>> np.testing.assert_equal(m.P, inst_restored.P) # doctest: +SKIP
set_model_params(P=None, pobs=None, pi=None, reversible=None, dt_model='1 step', neig=None)
Parameters
  • P (ndarray(m,m)) – coarse-grained or hidden transition matrix

  • pobs (ndarray (m,n)) – observation probability matrix from hidden to observable discrete states

  • pi (ndarray(m), optional, default=None) – stationary or distribution. Can be optionally given in case if it was already computed, e.g. by the estimator.

  • reversible (bool, optional, default=None) – whether P is reversible with respect to its stationary distribution. If None (default), will be determined from P

  • dt_model (str, optional, default='1 step') – time step of the model

  • neig (int or None) – The number of eigenvalues / eigenvectors to be kept. If set to None, all eigenvalues will be used.

simulate(N, start=None, stop=None, dt=1)

Generates a realization of the Hidden Markov Model

Parameters
  • N (int) – trajectory length in steps of the lag time

  • start (int, optional, default = None) – starting hidden state. If not given, will sample from the stationary distribution of the hidden transition matrix.

  • stop (int or int-array-like, optional, default = None) – stopping hidden set. If given, the trajectory will be stopped before N steps once a hidden state of the stop set is reached

  • dt (int) – trajectory will be saved every dt time steps. Internally, the dt’th power of P is taken to ensure a more efficient simulation.

Returns

  • htraj ((N/dt, ) ndarray) – The hidden state trajectory with length N/dt

  • otraj ((N/dt, ) ndarray) – The observable state discrete trajectory with length N/dt

sparse

Returns whether the MSM is sparse

stationary_distribution

The stationary distribution on the MSM states

submodel(states=None, obs=None)

Returns a HMM with restricted state space

Parameters
  • states (None or int-array) – Hidden states to restrict the model to (if not None).

  • obs (None, str or int-array) – Observed states to restrict the model to (if not None).

Returns

hmm – The restricted HMM.

Return type

HMM

timescales(k=None)

The relaxation timescales corresponding to the eigenvalues

Parameters

k (int) – number of timescales to be returned. By default all available eigenvalues, minus 1.

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_model

Physical time corresponding to one transition matrix step, e.g. ‘10 ps’

transition_matrix

The transition matrix on the active set.

transition_matrix_obs(k=1)

Computes the transition matrix between observed states

Transition matrices for longer lag times than the one used to parametrize this HMSM can be obtained by setting the k option. Note that a HMSM is not Markovian, thus we cannot compute transition matrices at longer lag times using the Chapman-Kolmogorow equality. I.e.:

\[P (k \tau) \neq P^k (\tau)\]

This function computes the correct transition matrix using the metastable (coarse) transition matrix \(P_c\) as:

\[P (k \tau) = {\Pi}^-1 \chi^{\top} ({\Pi}_c) P_c^k (\tau) \chi\]

where \(\chi\) is the output probability matrix, \(\Pi_c\) is a diagonal matrix with the metastable-state (coarse) stationary distribution and \(\Pi\) is a diagonal matrix with the observable-state stationary distribution.

Parameters

k (int, optional, default=1) – Multiple of the lag time for which the By default (k=1), the transition matrix at the lag time used to construct this HMSM will be returned. If a higher power is given,

update_model_params(**params)

Update given model parameter if they are set to specific values