Bootstraped MSM CIsΒΆ

quadwell_evaluated
In [1]:
import numpy as np
from matplotlib import pyplot as plt
from msmbuilder.example_datasets import QuadWell, quadwell_eigs
from msmbuilder.cluster import NDGrid
from msmbuilder.msm import MarkovStateModel
from sklearn.pipeline import Pipeline
In [2]:
dataset = QuadWell(random_state=0).get()
true_eigenvalues = quadwell_eigs(200)[0]
true_timescales = -1 / np.log(true_eigenvalues[1:])
print(QuadWell.description())

loading "/Users/rmcgibbo/msmbuilder_data/quadwell/version-1_random-state-0.pkl"...
This dataset consists of 100 trajectories simulated with Brownian dynamics
on the reduced potential function

V = 4(x^8 + 0.8 exp(-80 x^2) + 0.2 exp(-80 (x-0.5)^2) + 0.5 exp(-40 (x+0.5)^2)).

The simulations are governed by the stochastic differential equation

dx_t/dt = -\nabla V(x) + \sqrt{2D} * R(t),

where R(t) is a standard normal white-noise process, and D=1e3. The timsetep
is 1e-3. Each trajectory is 10^3 steps long, and starts from a random
initial point sampled from the uniform distribution on [-1, 1].
In [3]:
def msm_timescales(trajectories, n_states):
    pipeline = Pipeline([
        ('grid', NDGrid(min=-1.2, max=1.2)),
        ('msm', MarkovStateModel(n_timescales=4, reversible_type='transpose', verbose=False))
    ])
    pipeline.set_params(grid__n_bins_per_feature=n_states)
    pipeline.fit(trajectories)
    return pipeline.named_steps['msm'].timescales_

n_states = [5, 10, 50, 100]
ts = np.array([msm_timescales(dataset.trajectories, n) for n in n_states])
In [4]:
for i, c in enumerate(['b', 'r', 'm']):
    plt.plot(n_states, ts[:, i], c=c, marker='x')
    plt.axhline(true_timescales[i], ls='--', c=c, lw=2)

plt.xlabel('Number of states')
plt.ylabel('Timescale (steps)')
plt.show()
In [5]:

(quadwell.ipynb; quadwell_evaluated.ipynb; quadwell.py)

Versions