notebooks/04-msm-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 04 - MSM analysis\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\" title='This work is licensed under a Creative Commons Attribution 4.0 International License.' align=\"right\"/></a>\n",
    "\n",
    "In this notebook, we will cover how to analyze an MSM and how the modeled processes correspond to MSM spectral properties.\n",
    "We assume that you are familiar with data loading/visualization\n",
    "([Notebook 01 ➜ 📓](01-data-io-and-featurization.ipynb)),\n",
    "dimension reduction ([Notebook 02 ➜ 📓](02-dimension-reduction-and-discretization.ipynb)), and\n",
    "the estimation and validation process ([Notebook 03 ➜ 📓](03-msm-estimation-and-validation.ipynb)).\n",
    "\n",
    "Maintainers: [@cwehmeyer](https://github.com/cwehmeyer), [@marscher](https://github.com/marscher), [@thempel](https://github.com/thempel), [@psolsson](https://github.com/psolsson)\n",
    "\n",
    "**Remember**:\n",
    "- to run the currently highlighted cell, hold <kbd>&#x21E7; Shift</kbd> and press <kbd>&#x23ce; Enter</kbd>;\n",
    "- to get help for a specific function, place the cursor within the function's brackets, hold <kbd>&#x21E7; Shift</kbd>, and press <kbd>&#x21E5; Tab</kbd>;\n",
    "- you can find the full documentation at [PyEMMA.org](http://www.pyemma.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import mdshare\n",
    "import pyemma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1: preprocessed, two-dimensional data (toy model)\n",
    "We load the two-dimensional trajectory from an archive using numpy,\n",
    "directly discretize the full space using $k$-means clustering,\n",
    "visualize the marginal and joint distributions of both components as well as the cluster centers,\n",
    "and show the implied timescale (ITS) convergence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = mdshare.fetch('hmm-doublewell-2d-100k.npz', working_directory='data')\n",
    "with np.load(file) as fh:\n",
    "    data = fh['trajectory']\n",
    "\n",
    "cluster = pyemma.coordinates.cluster_kmeans(data, k=50, max_iter=50)\n",
    "its = pyemma.msm.its(\n",
    "    cluster.dtrajs, lags=[1, 2, 3, 5, 7, 10], nits=3, errors='bayes')\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n",
    "pyemma.plots.plot_feature_histograms(data, feature_labels=['$x$', '$y$'], ax=axes[0])\n",
    "pyemma.plots.plot_density(*data.T, ax=axes[1], cbar=False, alpha=0.1)\n",
    "axes[1].scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "axes[1].set_xlabel('$x$')\n",
    "axes[1].set_ylabel('$y$')\n",
    "axes[1].set_xlim(-4, 4)\n",
    "axes[1].set_ylim(-4, 4)\n",
    "axes[1].set_aspect('equal')\n",
    "pyemma.plots.plot_implied_timescales(its, ylog=False, ax=axes[2])\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show us the marginal (left panel) and joint distributions along with the cluster centers (middle panel).\n",
    "The implied timescales are converged (right panel). \n",
    "\n",
    "Before we proceed, let's have a look at the implied timescales error bars.\n",
    "They were computed from a Bayesian MSM, as requested by the `errors='bayes'` argument of the `pyemma.msm.its()` function.\n",
    "As mentioned before, Bayesian MSMs incorporate a sample of transition matrices.\n",
    "Target properties such as implied timescales can now simply be computed from the individual matrices.\n",
    "Thereby, the posterior distributions of these properties can be estimated.\n",
    "The ITS plot shows a confidence interval that contains $95\\%$ of the Bayesian samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayesian_msm = pyemma.msm.bayesian_markov_model(cluster.dtrajs, lag=1, conf=0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For any PyEMMA method that derives target properties from MSMs, sample mean and confidence intervals (as defined by the function argument above) are directly accessible with  `sample_mean()` and `sample_conf()`.\n",
    "Further, `sample_std()` is available for computing the standard deviation.\n",
    "In the more general case, it might be interesting to extract the full sample of a function evaluation with `sample_f()`.\n",
    "The syntax is equivalent for all those functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_mean = bayesian_msm.sample_mean('timescales', k=1)\n",
    "sample_conf_l, sample_conf_r = bayesian_msm.sample_conf('timescales', k=1)\n",
    "\n",
    "print('Mean of first ITS: {:f}'.format(sample_mean[0]))\n",
    "print('Confidence interval: [{:f}, {:f}]'.format(sample_conf_l[0], sample_conf_r[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that sample mean and maximum likelihood estimates are not identical and generally do not provide numerically identical results.\n",
    "\n",
    "Now, for the sake of simplicity we proceed with the analysis of a maximum likelihood MSM.\n",
    "We estimate it at lag time $1$ step..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msm = pyemma.msm.estimate_markov_model(cluster.dtrajs, lag=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and check for disconnectivity.\n",
    "The MSM is constructed on the largest set of discrete states that are (reversibly) connected.\n",
    "The `active_state_fraction` and `active_count_fraction` show us the fraction of discrete states and transition counts from our data which are part of this largest set and, thus, used for the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('fraction of states used = {:f}'.format(msm.active_state_fraction))\n",
    "print('fraction of counts used = {:f}'.format(msm.active_count_fraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fraction is, in both cases, $1$ and, thus, we have no disconnected states (which we would have to exclude from our analysis).\n",
    "\n",
    "If there were any disconnectivities in our data (fractions $<1$),\n",
    "we could access the indices of the **active states** (members of the largest connected set) via the `active_set` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(msm.active_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this potential issue out of the way, we can extract our first (stationary/thermodynamic) property,\n",
    "the `stationary_distribution` or, as a shortcut, `pi`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(msm.stationary_distribution)\n",
    "print('sum of weights = {:f}'.format(msm.pi.sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The attribute `msm.pi` tells us, for each discrete state, the absolute probability of observing said state in global equilibrium.\n",
    "Mathematically speaking, the stationary distribution $\\pi$ is the left eigenvector of the transition matrix $\\mathbf{P}$ to the eigenvalue $1$:\n",
    "\n",
    "$$\\pi^\\top \\mathbf{P} = \\pi^\\top.$$\n",
    "\n",
    "Please note that the $\\pi$ is fundamentaly different from a normalized histogram of states:\n",
    "for the histogram of states to accurately describe the stationary distribution, the data needs to be sampled from global equilibrium, i.e, the data points need to be statistically independent.\n",
    "The MSM approach, on the other hand, only requires local equilibrium, i.e., statistical independence of state transitions.\n",
    "Thus, the MSM approach requires a much weaker and, in practice, much easier to satisfy condition than simply counting state visits.\n",
    "\n",
    "We can use the stationary distribution to, e.g., visualize the weight of the dicrete states and, thus, to highlight which areas of our feature space are most probable.\n",
    "Here, we show all data points in a two dimensional scatter plot and color/weight them according to their discrete state membership:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax, misc = pyemma.plots.plot_contour(\n",
    "    *data.T, msm.pi[cluster.dtrajs[0]],\n",
    "    cbar_label='stationary distribution',\n",
    "    method='nearest', mask=True)\n",
    "ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_xlim(-4, 4)\n",
    "ax.set_ylim(-4, 4)\n",
    "ax.set_aspect('equal')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stationary distribution can also be used to correct the `pyemma.plots.plot_free_energy()` function that we used to visualize this dataset in [Notebook 01 ➜ 📓](01-data-io-and-featurization.ipynb).\n",
    "This might be necessary if the data points are not sampled from global equilibrium.\n",
    "\n",
    "In this case, we assign the weight of the corresponding discrete state to each data point and pass this information to the plotting function via its `weights` parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax, misc = pyemma.plots.plot_free_energy(\n",
    "    *data.T,\n",
    "    weights=np.concatenate(msm.trajectory_weights()),\n",
    "    legacy=False)\n",
    "ax.set_xlabel('$x$')\n",
    "ax.set_ylabel('$y$')\n",
    "ax.set_xlim(-4, 4)\n",
    "ax.set_ylim(-4, 4)\n",
    "ax.set_aspect('equal')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will see further uses of the stationary distribution later.\n",
    "But for now, we continue the analysis of our model by visualizing its (right) eigenvectors which encode the dynamical processes.\n",
    "First, we notice that the first right eigenvector is a constant $1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvec = msm.eigenvectors_right()\n",
    "print('first eigenvector is one: {} (min={}, max={})'.format(\n",
    "    np.allclose(eigvec[:, 0], 1, atol=1e-15), eigvec[:, 0].min(), eigvec[:, 0].max()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, the higher eigenvectors can be visualized as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    pyemma.plots.plot_contour(\n",
    "        *data.T, eigvec[cluster.dtrajs[0], i + 1], ax=ax, cmap='PiYG',\n",
    "        cbar_label='{}. right eigenvector'.format(i + 2), mask=True)\n",
    "    ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_xlim(-4, 4)\n",
    "    ax.set_ylim(-4, 4)\n",
    "    ax.set_aspect('equal')\n",
    "axes[0].set_ylabel('$y$')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right eigenvectors can be used to visualize the processes governed by the corresponding implied timescales.\n",
    "The first right eigenvector (always) is $(1,\\dots,1)^\\top$ for an MSM transition matrix and it corresponds to the stationary process (infinite implied timescale).\n",
    "\n",
    "The second right eigenvector corresponds to the slowest process;\n",
    "its entries are negative for one group of discrete states and positive for the other group.\n",
    "This tells us that the slowest process happens between these two groups and that the process relaxes on the slowest ITS ($\\approx 8.5$ steps).\n",
    "\n",
    "The third and fourth eigenvectors show a larger spread of values and no clear grouping.\n",
    "In combination with the ITS convergence plot, we can safely assume that these eigenvectors contain just noise and do not indicate any resolved processes.\n",
    "\n",
    "We then continue to validate our MSM with a CK test for $2$ metastable states which are already indicated by the second right eigenvector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nstates = 2\n",
    "pyemma.plots.plot_cktest(msm.cktest(nstates));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now save the model to do more analyses with PCCA++ and TPT in [Notebook 05 ➜ 📓](05-pcca-tpt.ipynb):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster.save('nb4.pyemma', model_name='doublewell_cluster', overwrite=True)\n",
    "msm.save('nb4.pyemma', model_name='doublewell_msm', overwrite=True)\n",
    "bayesian_msm.save('nb4.pyemma', model_name='doublewell_bayesian_msm', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 2: low-dimensional molecular dynamics data (alanine dipeptide)\n",
    "\n",
    "We fetch the alanine dipeptide data set, load the backbone torsions into memory,\n",
    "directly discretize the full space using $k$-means clustering,\n",
    "visualize the margial and joint distributions of both components as well as the cluster centers,\n",
    "and show the ITS convergence to help selecting a suitable lag time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pdb = mdshare.fetch('alanine-dipeptide-nowater.pdb', working_directory='data')\n",
    "files = mdshare.fetch('alanine-dipeptide-*-250ns-nowater.xtc', working_directory='data')\n",
    "\n",
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat.add_backbone_torsions(periodic=False)\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "cluster = pyemma.coordinates.cluster_kmeans(data, k=100, max_iter=50, stride=10)\n",
    "dtrajs_concatenated = np.concatenate(cluster.dtrajs)\n",
    "\n",
    "its = pyemma.msm.its(\n",
    "    cluster.dtrajs, lags=[1, 2, 5, 10, 20, 50], nits=4, errors='bayes')\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n",
    "pyemma.plots.plot_feature_histograms(\n",
    "    np.concatenate(data), feature_labels=['$\\Phi$', '$\\Psi$'], ax=axes[0])\n",
    "pyemma.plots.plot_density(*data_concatenated.T, ax=axes[1], cbar=False, alpha=0.1)\n",
    "axes[1].scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "axes[1].set_xlabel('$\\Phi$')\n",
    "axes[1].set_ylabel('$\\Psi$')\n",
    "pyemma.plots.plot_implied_timescales(its, ax=axes[2], units='ps')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plots show us the marginal (left panel) and joint distributions along with the cluster centers (middle panel).\n",
    "The implied timescales are converged (right panel). \n",
    "\n",
    "We then estimate an MSM at lag time $10$ ps and visualize the stationary distribution by coloring all data points according to the stationary weight of the discrete state they belong to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msm = pyemma.msm.estimate_markov_model(cluster.dtrajs, lag=10, dt_traj='1 ps')\n",
    "\n",
    "print('fraction of states used = {:f}'.format(msm.active_state_fraction))\n",
    "print('fraction of counts used = {:f}'.format(msm.active_count_fraction))\n",
    "\n",
    "fig, ax, misc = pyemma.plots.plot_contour(\n",
    "    *data_concatenated.T, msm.pi[dtrajs_concatenated],\n",
    "    cbar_label='stationary_distribution',\n",
    "    method='nearest', mask=True)\n",
    "ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "ax.set_xlabel('$\\Phi$')\n",
    "ax.set_ylabel('$\\Psi$')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we visualize the first six right eigenvectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigvec = msm.eigenvectors_right()\n",
    "print('first eigenvector is one: {} (min={}, max={})'.format(\n",
    "    np.allclose(eigvec[:, 0], 1, atol=1e-15), eigvec[:, 0].min(), eigvec[:, 0].max()))\n",
    "\n",
    "fig, axes = plt.subplots(2, 3, figsize=(12, 6))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    pyemma.plots.plot_contour(\n",
    "        *data_concatenated.T, eigvec[dtrajs_concatenated, i + 1], ax=ax, cmap='PiYG',\n",
    "        cbar_label='{}. right eigenvector'.format(i + 2), mask=True)\n",
    "    ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "    ax.set_xlabel('$\\Phi$')\n",
    "    ax.set_ylabel('$\\Psi$')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we have the $(1,\\dots,1)^\\top$ first right eigenvector of the stationary process.\n",
    "\n",
    "The second to fourth right eigenvectors illustrate the three slowest processes which are (in that order):\n",
    "\n",
    "- rotation of the $\\Phi$ dihedral,\n",
    "- rotation of the $\\Psi$ dihedral when $\\Phi\\approx-2$ rad, and\n",
    "- rotation of the $\\Psi$ dihedral when $\\Phi\\approx1$ rad.\n",
    "\n",
    "Eigenvectors five, six, and seven indicate further processes which, however, relax faster than the lag time and cannot be resolved clearly.\n",
    "\n",
    "We now proceed our validation process using a Bayesian MSM with four metastable states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nstates = 4\n",
    "bayesian_msm = pyemma.msm.bayesian_markov_model(cluster.dtrajs, lag=10, dt_traj='1 ps')\n",
    "pyemma.plots.plot_cktest(bayesian_msm.cktest(nstates), units='ps');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that four metastable states are a reasonable choice for our MSM.\n",
    "\n",
    "In general, the number of metastable states is a modeler's choice; it is adjusted to map the kinetics to be modeled.\n",
    "In the current example, increasing the resolution with a higher number of metastable states or resolving only the slowest process between $2$ states would be possible.\n",
    "However, the number of states is not arbitrary as the observed processes in metastable state space need not be Markovian in general.\n",
    "A failed Chapman-Kolmogorov test can thus also hint to a bad choice of the metastable state number.\n",
    "\n",
    "In order to perform further analysis, we save the model to disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster.save('nb4.pyemma', model_name='ala2_cluster', overwrite=True)\n",
    "msm.save('nb4.pyemma', model_name='ala2_msm', overwrite=True)\n",
    "bayesian_msm.save('nb4.pyemma', model_name='ala2_bayesian_msm', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1\n",
    "Load the heavy atom distances into memory, TICA (`lag=3` and `dim=2`), discretize with  100 $k$-means centers and a stride of $10$, and show the ITS convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = #FIXME\n",
    "feat. #FIXME\n",
    "data = #FIXME\n",
    "\n",
    "tica = #FIXME\n",
    "tica_concatenated = #FIXME\n",
    "\n",
    "cluster = #FIXME\n",
    "dtrajs_concatenated = #FIXME\n",
    "\n",
    "its = #FIXME\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n",
    "pyemma.plots.plot_feature_histograms(tica_concatenated, feature_labels=['IC 1', 'IC 2'], ax=axes[0])\n",
    "pyemma.plots.plot_density(*tica_concatenated.T, ax=axes[1], cbar=False, alpha=0.3)\n",
    "axes[1].scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "axes[1].set_xlabel('IC 1')\n",
    "axes[1].set_ylabel('IC 2')\n",
    "pyemma.plots.plot_implied_timescales(its, ax=axes[2], units='ps')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "pairs = feat.pairs(feat.select_Heavy())\n",
    "feat.add_distances(pairs, periodic=False)\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "tica = pyemma.coordinates.tica(data, lag=3, dim=2)\n",
    "tica_concatenated = np.concatenate(tica.get_output())\n",
    "\n",
    "cluster = pyemma.coordinates.cluster_kmeans(tica, k=100, max_iter=50, stride=10)\n",
    "dtrajs_concatenated = np.concatenate(cluster.dtrajs)\n",
    "\n",
    "its = pyemma.msm.its(\n",
    "    cluster.dtrajs, lags=[1, 2, 5, 10, 20, 50], nits=4, errors='bayes')\n",
    "\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n",
    "pyemma.plots.plot_feature_histograms(tica_concatenated, feature_labels=['IC 1', 'IC 2'], ax=axes[0])\n",
    "pyemma.plots.plot_density(*tica_concatenated.T, ax=axes[1], cbar=False, alpha=0.3)\n",
    "axes[1].scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "axes[1].set_xlabel('IC 1')\n",
    "axes[1].set_ylabel('IC 2')\n",
    "pyemma.plots.plot_implied_timescales(its, ax=axes[2], units='ps')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 2\n",
    "Estimate an MSM at lag time $10$ ps with `dt_traj='1 ps'` and visualize the stationary distribution using a two-dimensional colored scatter plot of all data points in TICA space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "msm = #FIXME\n",
    "\n",
    "print('fraction of states used = {:f}'. #FIXME\n",
    "print('fraction of counts used = {:f}'. #FIXME\n",
    "\n",
    "fig, ax, misc = pyemma.plots.plot_contour(\n",
    "    *tica_concatenated.T, msm.pi[dtrajs_concatenated],\n",
    "    cbar_label='stationary_distribution',\n",
    "    method='nearest', mask=True)\n",
    "ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "ax.set_xlabel('IC 1')\n",
    "ax.set_ylabel('IC 2')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "msm = pyemma.msm.estimate_markov_model(cluster.dtrajs, lag=10, dt_traj='1 ps')\n",
    "\n",
    "print('fraction of states used = {:f}'.format(msm.active_state_fraction))\n",
    "print('fraction of counts used = {:f}'.format(msm.active_count_fraction))\n",
    "\n",
    "fig, ax, misc = pyemma.plots.plot_contour(\n",
    "    *tica_concatenated.T, msm.pi[dtrajs_concatenated],\n",
    "    cbar_label='stationary_distribution',\n",
    "    method='nearest', mask=True)\n",
    "ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "ax.set_xlabel('IC 1')\n",
    "ax.set_ylabel('IC 2')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 3\n",
    "Visualize the first six right eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "eigvec = #FIXME\n",
    "print('first eigenvector is one: {} (min={}, max={})'.format( #FIXME\n",
    "\n",
    "fig, axes = plt.subplots(2, 3, figsize=(12, 6))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    pyemma.plots.plot_contour( #FIXME )\n",
    "    ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "    ax.set_xlabel('IC 1')\n",
    "    ax.set_ylabel('IC 2')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "eigvec = msm.eigenvectors_right()\n",
    "print('first eigenvector is one: {} (min={}, max={})'.format(\n",
    "    np.allclose(eigvec[:, 0], 1, atol=1e-15), eigvec[:, 0].min(), eigvec[:, 0].max()))\n",
    "\n",
    "fig, axes = plt.subplots(2, 3, figsize=(12, 6))\n",
    "for i, ax in enumerate(axes.flat):\n",
    "    pyemma.plots.plot_contour(\n",
    "        *tica_concatenated.T, eigvec[dtrajs_concatenated, i + 1], ax=ax, cmap='PiYG',\n",
    "        cbar_label='{}. right eigenvector'.format(i + 2), mask=True)\n",
    "    ax.scatter(*cluster.clustercenters.T, s=15, c='C1')\n",
    "    ax.set_xlabel('IC 1')\n",
    "    ax.set_ylabel('IC 2')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can you already guess from eigenvectors two to four which the metastable states are?\n",
    "\n",
    "#### Exercise 4\n",
    "Estimate a Bayesian MSM at lag time $10$ ps and perform/show a CK test for four metastable states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "bayesian_msm = #FIXME\n",
    "\n",
    "nstates = 4\n",
    "pyemma.plots. #FIXME"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "bayesian_msm = pyemma.msm.bayesian_markov_model(cluster.dtrajs, lag=10, dt_traj='1 ps')\n",
    "\n",
    "nstates = 4\n",
    "pyemma.plots.plot_cktest(bayesian_msm.cktest(nstates), units='ps');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 5\n",
    "Save the MSM, Bayesian MSM and Cluster objects to the same file as before.\n",
    "Use the model names `ala2tica_msm`, `ala2tica_bayesian_msm` and `ala2tica_cluster`, respectively.\n",
    "Further, include the TICA object with model name `ala2tica_tica`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "#FIXME "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "cluster.save('nb4.pyemma', model_name='ala2tica_cluster', overwrite=True)\n",
    "msm.save('nb4.pyemma', model_name='ala2tica_msm', overwrite=True)\n",
    "bayesian_msm.save('nb4.pyemma', model_name='ala2tica_bayesian_msm', overwrite=True)\n",
    "tica.save('nb4.pyemma', model_name='ala2tica_tica', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrapping up\n",
    "In this notebook, we have learned how to analyze an MSM and how to extract kinetic information from the model. In detail, we have used\n",
    "- the `active_state_fraction`, `active_count_fraction`, and `active_set` attributes of an MSM object to see how much (and which parts) of our data form the largest connected set represented by the MSM,\n",
    "- the `stationary_distribution` (or `pi`) attribute of an MSM object to access its stationary vector,\n",
    "- the `eigenvectors_right()` method of an MSM object to access its (right) eigenvectors,\n",
    "\n",
    "For visualizing MSMs or kinetic networks we used\n",
    "- `pyemma.plots.plot_density()`\n",
    "- `pyemma.plots.plot_contour()` and\n",
    "- `pyemma.plots.plot_cktest()`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}