notebooks/01-data-io-and-featurization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 01 - Data-I/O and featurization\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 load (and visualize) data with PyEMMA.\n",
    "We are going to extract different features (collective variables)\n",
    "and compare how well these are suited for Markov state model building.\n",
    "Further, we will introduce the concept of streaming data, which is mandatory to work with large data sets.\n",
    "\n",
    "As with the other notebooks, we employ multiple examples.\n",
    "The idea is, first, to highlight the fundamental ideas with a non-physical test system (diffusion in a 2D double-well potential)\n",
    "and, second, to demonstrate real-world applications with molecular dynamics data.\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 numpy as np\n",
    "import mdshare\n",
    "import pyemma\n",
    "# for visualization of molecular structures:\n",
    "import nglview\n",
    "import mdtraj\n",
    "from threading import Timer\n",
    "from nglview.player import TrajectoryPlayer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 1: preprocessed data (toy model)\n",
    "\n",
    "In the most convenient case, we already have preprocessed time series data available in some kind of archive which can be read using numpy:"
   ]
  },
  {
   "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",
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the archive contains multiple arrays\n",
    "and we store the array which belongs to the `trajectory` key into the variable `data`.\n",
    "To see which arrays/keys are available, you can `print` the output of `fh.keys()`.\n",
    "\n",
    "Please note that mdshare is a package designed for distributing our tutorial data:\n",
    "the `mdshare.fetch()` call ensures that the file `data/hmm-doublewell-2d-100k.npz` exists\n",
    "and passes this relative path to the variable `file`.\n",
    "In practice, we could have written\n",
    "```Python\n",
    "file = 'data/hmm-doublewell-2d-100k.npz'\n",
    "```\n",
    "provided that this file actually exists.\n",
    "\n",
    "Once we have the data in memory,\n",
    "we can use one of PyEMMA's plotting functions to visualize what we have loaded:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyemma.plots.plot_feature_histograms(data, feature_labels=['$x$', '$y$']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `plot_feature_histograms()` function visualizes the distributions of all degrees of freedom\n",
    "assuming that the columns of `data` represent different features and the rows represent different time steps.\n",
    "\n",
    "While `plot_feature_histograms()` can handle arbitrary numbers of features,\n",
    "we have additional plotting functions for 2D projections.\n",
    "First, we visualize the sample density in the $x/y$-plane..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax, misc = pyemma.plots.plot_density(*data.T)\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": [
    "... then, we show the corresponding pseudo free energy surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax, misc = pyemma.plots.plot_free_energy(*data.T, 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": [
    "Both functions make a two-dimensional histogram for the given features;\n",
    "the free energy surface is defined by the negative logarithm of the probability computed from the histogram counts.\n",
    "\n",
    "⚠️ Please note that these functions visualize the density and free energy of the **sampled data**,\n",
    "not the equilibrium distribution of the underlying system.\n",
    "To account for nonequiblibrium data, you can supply frame-wise weights using the `weights` parameter.\n",
    "This will be covered in the follow-up [Notebook 04 ➜ 📓](04-msm-analysis.ipynb).\n",
    "\n",
    "## Case 2: loading `*.xtc` files (alanine dipeptide)\n",
    "\n",
    "To load molecular dynamics data from one of the standard file formats ( `*.xtc`),\n",
    "we need not only the actual simulation data, but a topology file, too.\n",
    "This might differ for other formats though."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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",
    "print(pdb)\n",
    "print(files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a look at the structure with the aid of NGLView.\n",
    "We load the PDB file into memory with mdtraj and visualize it.\n",
    "The widget will auto-close after 30 second;\n",
    "if you want to watch it again, please execute the cell below again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "widget = nglview.show_mdtraj(mdtraj.load(pdb))\n",
    "p = TrajectoryPlayer(widget)\n",
    "widget.add_ball_and_stick()\n",
    "p.spin = True\n",
    "def stop_spin():\n",
    "    p.spin = False\n",
    "    widget.close()\n",
    "Timer(30, stop_spin).start()\n",
    "widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In PyEMMA, the `featurizer` is a central object that incorporates the system's topology.\n",
    "We start by creating it using the topology file.\n",
    "Features are now easily computed by adding the target feature.\n",
    "If no feature is added, the featurizer will extract Cartesian coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we pass the featurizer to the load function to extract the Cartesian coordinates from the trajectory files into memory.\n",
    "For real world examples one would prefer the `source()` function,\n",
    "because usually one has more data available than main memory in the workstation.\n",
    "\n",
    "The warning about **plain coordinates** is triggered,\n",
    "because these coordinates will include diffusion as a dynamical process,\n",
    "which might not be what one is interested in.\n",
    "If the molecule of interest has been aligned to a reference prior the analysis,\n",
    "it is fine to use these coordinates, but we will see that there are better choices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "print('type of data:', type(data))\n",
    "print('lengths:', len(data))\n",
    "print('shape of elements:', data[0].shape)\n",
    "print('n_atoms:', feat.topology.n_atoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we start adding features which we want to extract from the simulation data.\n",
    "Here, we want to load the backbone torsions,\n",
    "because these angles are known to describe all flexibility in the system.\n",
    "Since this feature is two dimensional, it is also easier to visualize.\n",
    "Please note that, in complex systems, it is not trivial to visualize plain input features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat.add_backbone_torsions(periodic=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️ Please note that the trajectories have been aligned to a reference structure before.\n",
    "Since in that case we loose track of the periodic box,\n",
    "we have to switch off the `periodic` flag for the torsion angle computations.\n",
    "By default PyEMMA assumes your simulation uses periodic boundary conditions.\n",
    "\n",
    "We can always call the featurizer's `describe()` method to show which features are requested.\n",
    "You might have noticed that you can combine arbitrary features by having multiple calls to `add_` methods of the featurizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "print('type of data:', type(data))\n",
    "print('lengths:', len(data))\n",
    "print('shape of elements:', data[0].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we have selected all desired features,\n",
    "we can call the `load()` function to load all features into memory or,\n",
    "alternatively, the `source()` function to create a streamed feature reader.\n",
    "For now, we will use `load()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(feat.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently, we have loaded a list of three two-dimensional `numpy.ndarray` objects from our three trajectory files.\n",
    "We can visualize these features using the aforementioned plotting functions,\n",
    "but to do so we have to concatenate the three individual trajectories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_concatenated = np.concatenate(data)\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now measure the quantity of kinetic variance of the just selected feature by computing a VAMP-2 score <a id=\"ref-1\" href=\"#cite-vamp-preprint\">wu-17</a>.\n",
    "This enables us to distinguish features on how well they might be suited for MSM building.\n",
    "The minimum value of this score is $1$, which corresponds to the invariant measure or the observed equilibrium.\n",
    "\n",
    "With the dimension parameter we specify the amount of dynamic processes that we want to score.\n",
    "This is of importance later on, when we want to compare different input features\n",
    "([Notebook 02 ➜ 📓](02-dimension-reduction-and-discretization.ipynb)).\n",
    "If we did not fix this number, we would not have an upper bound for the score.\n",
    "\n",
    "⚠️ Please also note that we split our available data into training and test sets,\n",
    "where we leave out the last file in training and then use it as test.\n",
    "This is an important aspect in practice to avoid over-fitting the score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_phi_psi = pyemma.coordinates.vamp(data[:-1], dim=2).score(\n",
    "        test_data=data[-1],\n",
    "        score_method='VAMP2')\n",
    "print('VAMP2-score backbone torsions: {:.2f}'.format(score_phi_psi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The score of $1.5$ means that we have the constant of $1$ plus a total contribution of $0.5$ from the other dynamic processes.\n",
    "\n",
    "We now use PyEMMA's `plot_density()` and `plot_free_energy()` functions to create Ramachandran plots of our system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharex=True, sharey=True)\n",
    "# the * operator used in a function call is used to unpack\n",
    "# the iterable variable into its single elements. \n",
    "pyemma.plots.plot_density(*data_concatenated.T, ax=axes[0])\n",
    "pyemma.plots.plot_free_energy(*data_concatenated.T, ax=axes[1], legacy=False)\n",
    "for ax in axes.flat:\n",
    "    ax.set_xlabel('$\\Phi$')\n",
    "    ax.set_aspect('equal')\n",
    "axes[0].set_ylabel('$\\Psi$')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that the distribution in backbone torsion space contains several basins that will be assigned to metastable states in follow-up notebooks ([Notebook 05 ➜ 📓](05-pcca-tpt.ipynb), [Notebook 07 ➜ 📓](07-hidden-markov-state-models.ipynb)).\n",
    "Again, the free energy plot only depicts a pseudo free energy surface of the sampled data and was not re-weighted to equilibrium.\n",
    "\n",
    "Let us look at a different featurization example and load the positions of all heavy atoms instead.\n",
    "We create a new featurizer object and use its `add_selection()` method to request the positions of a given selection of atoms.\n",
    "For this selection, we can use the `select_Heavy()` method which returns the indices of all heavy atoms.\n",
    "\n",
    "Again, we load the data into memory and show what we loaded using the `describe()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat.add_selection(feat.select_Heavy())\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "feat.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️ Please note that PyEMMA has flattened the $x, y$ and $z$ coordinates into an array that will be used for further analysis.\n",
    "\n",
    "We visualize the distributions of all loaded features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(np.concatenate(data), feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we have a look at the VAMP-2 score of the heavy atom coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_heavy_atoms = pyemma.coordinates.vamp(data[:-1], dim=2).score(\n",
    "        test_data=data[:-1],\n",
    "        score_method='VAMP2')\n",
    "print('VAMP2-score backbone torsions: {:f}'.format(score_phi_psi))\n",
    "print('VAMP2-score xyz: {:f}'.format(score_heavy_atoms))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see, the score for the heavy atom positions is much higher as the one for the $\\phi/\\psi$ torsion angles.\n",
    "The feature with a higher score should be favored for further analysis,\n",
    "because it means that this feature contains more information about slow processes.\n",
    "If you are already digging deeper into your system of interest,\n",
    "you can of course restrict the analysis to a set of features you already know describes your processes of interest,\n",
    "regardless of its VAMP score.\n",
    "\n",
    "Another featurization that is interesting especially for proteins are pairwise heavy atom distances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "heavy_atom_distance_pairs = feat.pairs(feat.select_Heavy())\n",
    "feat.add_distances(heavy_atom_distance_pairs, periodic=False)\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "feat.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us compare the score of heavy atom distance pairs to the other scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_pair_heavy_atom_dists = pyemma.coordinates.vamp(data[:-1], dim=2).score(\n",
    "    test_data=data[-1],\n",
    "    score_method='VAMP2')\n",
    "print('VAMP2-score: {:.2f}'.format(score_pair_heavy_atom_dists))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently, the heavy atom distance pairs cover an amount of kinetic variance which is comparable to the coordinates themselves.\n",
    "However we probably would not use much information by using this internal degree of freedom,\n",
    "while avoiding the need to align our trajectories first.\n",
    "\n",
    "### `load()` versus `source()`\n",
    "\n",
    "Using `load()`, we put the full data into memory.\n",
    "This is possible for all examples in this tutorial.\n",
    "\n",
    "Many real world applications, though, require more memory than your workstation might provide.\n",
    "For these cases, you should use the `source()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reader = pyemma.coordinates.source(files, features=feat)\n",
    "print(reader)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function creates a reader, wich allows to stream the underlying data in chunks instead of the full set.\n",
    "Most of the functions in the `coordinates` sub-package accept data in memory as well as readers.\n",
    "However, some plotting functions require the data to be in memory.\n",
    "To load a (sub-sampled) subset into memory, we can use the `get_output()` method with a stride parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_output = reader.get_output(stride=5)\n",
    "len(data_output)\n",
    "print('number of frames in first file: {}'.format(reader.trajectory_length(0)))\n",
    "print('number of frames after striding: {}'.format(len(data_output[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have loaded every fifth frame into memory.\n",
    "Again, we can visualize the (concatenated) features with a feature histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(\n",
    "    np.concatenate(data_output), feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing your progress\n",
    "\n",
    "In the remainder of this notebook, you will find short exercises to test your newly learned skills.\n",
    "The exercises are announced by the keyword **Exercise** and followed by an incomplete cell.\n",
    "Missing parts are indicated by\n",
    "```python\n",
    "#FIXME\n",
    "```\n",
    "Exercise cells come with a button (**Show Solution**) to reveal the solution.\n",
    "\n",
    "#### Exercise 1: inverse heavy atom distances\n",
    "\n",
    "Please fix the following code block such that the inverse distances between all heavy atoms are loaded and visualized.\n",
    "\n",
    "**Hint**: try to use the auto-complete feature on the feat object to gain some insight.\n",
    "Also take a look at the previous demonstrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "pairs = feat.pairs(# FIXME)\n",
    "feat. #FIXME\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(\n",
    "    np.concatenate(data), feature_labels=feat, ax=ax)\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_inverse_distances(pairs)\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(\n",
    "    np.concatenate(data), feature_labels=feat, ax=ax)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Excercise 2: compare the inverse distances feature to the previously computed ones\n",
    "\n",
    "Compute and discuss a cross-validated VAMP score for the inverse pairwise heavy atom distances and plot the result. What do you observe and which feature would you choose for further processing?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "score_inv_dist = pyemma.coordinates. #FIXME\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "score_mapping = dict(score_heavy_atoms=score_heavy_atoms,\n",
    "                     score_phi_psi=score_phi_psi,\n",
    "                     score_pair_heavy_atom_dists=score_pair_heavy_atom_dists,\n",
    "                     score_inv_dist=score_inv_dist)\n",
    "lbl = []\n",
    "for i, (key, value) in enumerate(sorted(score_mapping.items(), key=lambda x: x[1])):\n",
    "    ax.bar(i, height=value)\n",
    "    lbl.append(key)\n",
    "ax.set_xticks(np.arange(0, len(score_mapping), 1))\n",
    "ax.set_xticklabels(lbl)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "###### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "score_inv_dist = pyemma.coordinates.vamp(\n",
    "    data[:-1], dim=2).score(test_data=data[-1])\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "score_mapping = dict(score_heavy_atoms=score_heavy_atoms,\n",
    "                     score_phi_psi=score_phi_psi,\n",
    "                     score_pair_heavy_atom_dists=score_pair_heavy_atom_dists,\n",
    "                     score_inv_dist=score_inv_dist)\n",
    "lbl = []\n",
    "for i, (key, value) in enumerate(sorted(score_mapping.items(), key=lambda x: x[1])):\n",
    "    ax.bar(i, height=value)\n",
    "    lbl.append(key)\n",
    "ax.set_xticks(np.arange(0, len(score_mapping), 1))\n",
    "ax.set_xticklabels(lbl)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case 3: loading `*.xtc` files (pentapeptide)\n",
    "\n",
    "Once we have obtained the raw data files..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdb = mdshare.fetch('pentapeptide-impl-solv.pdb', working_directory='data')\n",
    "files = mdshare.fetch('pentapeptide-*-500ns-impl-solv.xtc', working_directory='data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and had a quick look at the structure again..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "widget = nglview.show_mdtraj(mdtraj.load(pdb))\n",
    "p = TrajectoryPlayer(widget)\n",
    "widget.add_ball_and_stick()\n",
    "p.spin = True\n",
    "def stop_spin():\n",
    "    p.spin = False\n",
    "    widget.close()\n",
    "Timer(30, stop_spin).start()\n",
    "widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... we can load a selection of features into memory.\n",
    "Here, we want the $\\cos/\\sin$ transformations of the backbone and $\\chi_1$ sidechain torsions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat.add_backbone_torsions(cossin=True, periodic=False)\n",
    "feat.add_sidechain_torsions(which='chi1', cossin=True, periodic=False)\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "print(feat.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we visualize the (concatenated) features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises: feature selection and visualization\n",
    "\n",
    "#### Exercise 3\n",
    "\n",
    "Complete the following code block to load/visualize the distances between all $\\text{C}_\\alpha$ carbon atoms.\n",
    "\n",
    "**Hint**: You might find the `add_distances_ca()` method of the featurizer object helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat. #FIXME\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\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",
    "feat.add_distances_ca(periodic=False)\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 4\n",
    "\n",
    "Complete the following code block to load/visualize the minimal distances between all residues (excluding first and second neighbors).\n",
    "\n",
    "**Hint**: You might find the `add_residue_mindist()` method of the featurizer object helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat. #FIXME\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\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",
    "feat.add_residue_mindist(periodic=False)\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 5\n",
    "\n",
    "Complete the following code block to load/visualize the position of all backbone atoms.\n",
    "\n",
    "**Hint**: You might find the `select_Backbone()` method of the featurizer object helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat. #FIXME\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\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",
    "feat.add_selection(feat.select_Backbone())\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 12))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 6\n",
    "\n",
    "Complete the following code block to load/visualize the position of all $\\text{C}_\\alpha$ atoms.\n",
    "\n",
    "**Hint**: You might find the `select_Ca()` method of the featurizer object helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "format": "row",
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "feat = pyemma.coordinates.featurizer(pdb)\n",
    "feat. #FIXME\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\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",
    "feat.add_selection(feat.select_Ca())\n",
    "\n",
    "data = pyemma.coordinates.load(files, features=feat)\n",
    "data_concatenated = np.concatenate(data)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "pyemma.plots.plot_feature_histograms(data_concatenated, feature_labels=feat, ax=ax)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wrapping up\n",
    "\n",
    "In this notebook, we have learned how to load and visualize molecular simulation data with PyEMMA.\n",
    "In detail, we have used\n",
    "- `pyemma.coordinates.featurizer()` to define a selection of features we want to extract,\n",
    "- `pyemma.coordinates.load()` to load data into memory, and\n",
    "- `pyemma.coordinates.source()` to create a streamed feature reader in case the data does not fit into memory.\n",
    "\n",
    "After loading the data into memory, we have used\n",
    "- `pyemma.coordinates.vamp().score()` to score the quality of the features,\n",
    "- `pyemma.plots.plot_feature_histograms()` to show the distributions of all loaded features,\n",
    "- `pyemma.plots.plot_density()` to visualize the sample density, and\n",
    "- `pyemma.plots.plot_free_energy()` to visualize the free energy surface of two selected features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "<a id=\"cite-vamp-preprint\"/><sup><a href=#ref-1>[^]</a></sup>Wu, H. and Noé, F.. 2017. _Variational approach for learning Markov processes from time series data_. [URL](https://arxiv.org/pdf/1707.04659.pdf)\n",
    "\n"
   ]
  }
 ],
 "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
}