Skip to content

msm

Create Markov State Model.

This submodule contains all methods related to estimate the Markov state model.

estimate_markov_model(trajs, lagtime)

Estimates Markov State Model.

This method estimates the MSM based on the transition count matrix.

Parameters:

  • trajs (StateTraj or list or ndarray or list of ndarray) –

    State trajectory/trajectories used to estimate the MSM.

  • lagtime (int) –

    Lag time for estimating the markov model given in [frames].

Returns:

  • T ( ndarray ) –

    Transition probability matrix \(T_{ij}\), containing the transition probability transition from state \(i o j\).

  • states ( ndarray ) –

    Array holding states corresponding to the columns of \(T_{ij}\).

Source code in src/msmhelper/msm/msm.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def estimate_markov_model(trajs, lagtime):
    """Estimates Markov State Model.

    This method estimates the MSM based on the transition count matrix.

    Parameters
    ----------
    trajs : StateTraj or list or ndarray or list of ndarray
        State trajectory/trajectories used to estimate the MSM.
    lagtime : int
        Lag time for estimating the markov model given in [frames].

    Returns
    -------
    T : ndarray
        Transition probability matrix $T_{ij}$, containing the transition
        probability transition from state $i\to j$.
    states : ndarray
        Array holding states corresponding to the columns of $T_{ij}$.

    """
    trajs = StateTraj(trajs)
    return trajs.estimate_markov_model(lagtime)

row_normalize_matrix(mat)

Row normalize the given 2d matrix.

Parameters:

  • mat (ndarray) –

    Matrix to be row normalized.

Returns:

  • mat ( ndarray ) –

    Normalized matrix.

Source code in src/msmhelper/msm/msm.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
@numba.njit
def row_normalize_matrix(mat):
    """Row normalize the given 2d matrix.

    Parameters
    ----------
    mat : ndarray
        Matrix to be row normalized.

    Returns
    -------
    mat : ndarray
        Normalized matrix.

    """
    row_sum = np.sum(mat, axis=1)
    if not row_sum.all():
        row_sum[row_sum == 0] = 1

    # due to missing np.newaxis row_sum[:, np.newaxis] becomes # noqa: SC100
    return mat / row_sum.reshape(mat.shape[0], 1)

equilibrium_population(tmat, allow_non_ergodic=True)

Calculate equilibirum population.

If there are non ergodic states, their population is set to zero.

Parameters:

  • tmat (ndarray) –

    Quadratic transition matrix, needs to be ergodic.

  • allow_non_ergodic (bool, default: True ) –

    If True only the largest ergodic subset will be used. Otherwise it will throw an error if not ergodic.

Returns:

  • peq ( ndarray ) –

    Equilibrium population of input matrix.

Source code in src/msmhelper/msm/msm.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
@decorit.alias('peq')
def equilibrium_population(tmat, allow_non_ergodic=True):
    """Calculate equilibirum population.

    If there are non ergodic states, their population is set to zero.

    Parameters
    ----------
    tmat : ndarray
        Quadratic transition matrix, needs to be ergodic.
    allow_non_ergodic : bool
        If True only the largest ergodic subset will be used. Otherwise it will
        throw an error if not ergodic.

    Returns
    -------
    peq : ndarray
        Equilibrium population of input matrix.

    """
    tmat = np.asarray(tmat)
    is_ergodic = tests.is_ergodic(tmat)
    if not allow_non_ergodic and not is_ergodic:
        raise ValueError('tmat needs to be ergodic transition matrix.')

    # calculate ev for ergodic subset
    if is_ergodic:
        _, eigenvectors = linalg.left_eigenvectors(tmat, nvals=1)
        eigenvectors = eigenvectors[0]
    else:
        mask = tests.ergodic_mask(tmat)
        _, evs_mask = linalg.left_eigenvectors(
            row_normalize_matrix(
                tmat[np.ix_(mask, mask)],
            ),
            nvals=1,
        )

        eigenvectors = np.zeros(len(tmat), dtype=tmat.dtype)
        eigenvectors[mask] = evs_mask[0]

    return eigenvectors / np.sum(eigenvectors)