Cluster Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster Analysis of Molecular Dynamics Trajectory Data\n",
    "\n",
    "*A tutorial/workshop created by [Charlie Laughton](mailto:charles.laughton@nottingham.ac.uk)  for [CCP-BioSim](http://www.ccpbiosim.ac.uk)*\n",
    "\n",
    "### Introduction\n",
    "By the \"jumping amongst minima\" model of MD dynamics, we expect large ensembles/long trajectories to show clustering in conformational space. Cluster centres correspond to local minima conformations, and the sizes and shapes of the clusters tell us about the shape of the free energy surface.\n",
    "\n",
    "The theory and practice of clustering is a massive subject that we can't cover much here; for a bit of an introduction, see [this Towards Data Science article](https://towardsdatascience.com/the-5-clustering-algorithms-data-scientists-need-to-know-a36d136ef68) and for a molecular-dynamics focussed discussion, [this article](https://pubs.acs.org/doi/10.1021/ct700119m) from the Cheatham group.\n",
    "\n",
    "Here we illustrate how some of the clustering methods in the [scipy](https://scipy.org) Python package can be applied to MD trajectories/ensembles, and also look at another method developed by Rodriguez and Laio \n",
    "[(*Science*, 2014,  **344**, 1492-1495)](http://science.sciencemag.org/content/344/6191/1492) that seems - often - to work pretty well for MD type data.\n",
    "\n",
    "A feature of all the clustering methods used here is that they are \"agglomerative\". This means there is no up-front decision by the user as to how many clusters there should be, which makes sense for MD-type applications where you hope the data will tell you where the clusters are. However in general there is still a point in the process where a human eye has to make a decision as to whether something is a cluster or not - you will see that here.\n",
    "\n",
    "The data set we will investigate consists of 10,000 snapshots from a 100 nanosecond implicit solvent simulation of the \"alanine dipeptide\". We use PCA to help us visualise the distribution of the data in conformational space, and so how the clustering methods perform.\n",
    "\n",
    "#### Requirements\n",
    "\n",
    "If you are running this tutorial on the CCP-BioSim training cloud all required Python packages will be present. If you have cloned the repo to run the Jupyter notebook locally, you will need ('pip install'):\n",
    "\n",
    "1. mdtraj\n",
    "2. mdplus\n",
    "3. numpy\n",
    "4. scipy\n",
    "5. rlcluster\n",
    "6. matplotlib\n",
    "\n",
    "-------\n",
    "### Step 1: import the required libraries:\n",
    "\n",
    "1. MDTraj to load and manipulate the trajectory data.\n",
    "2. Matplotlib for the graphs\n",
    "3. MDPlus for PCA\n",
    "4. Scipy for clustering methods\n",
    "5. RLCluster for Rodrigues-Laio clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mdtraj as mdt\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "from mdplus.pca import PCA\n",
    "from scipy.spatial.distance import pdist, cdist, squareform\n",
    "from scipy.cluster.hierarchy import ward, fcluster, single, complete\n",
    "import rlcluster\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Load the trajectory data and strip down to heavy atoms only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = mdt.load('ala2.md.nc', top='ala2.prmtop')\n",
    "heavy = t.topology.select('mass > 2.0')\n",
    "t = t.atom_slice(heavy)\n",
    "print(t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Find a way of visualizing the data that will help you evaluate the performance of your clustering methods. \n",
    "Here we will use Principal Component Analysis of the snapshots to create a low-dimensional representation of their distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = PCA()\n",
    "scores = p.fit_transform(t.xyz)\n",
    "plt.bar(range(1,11), p.eigenvalues[:10])\n",
    "plt.xlabel('eigenvector #')\n",
    "plt.ylabel('explained variance (nm**2)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scores` matrix generated by the PCA process contains one row per snapshot. Each row contains the coordinates of the structure in the 30-dimensional PCA space. From the eigenvalue analysis above, we see that only the first 3 or so values in each row are \"interesting\" in that they show significant variance across the data set.\n",
    "\n",
    "Create a utility function to plot data in the low-dimensional PCA space (it will save a lot of repeated code later). Create three views: one in the PC1/PC2 plane, one in the PC1/PC3 plane, and one in the PC2/PC3 plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_distribution(axs, score_sets, title=None):\n",
    "    \"\"\"\n",
    "    Plot selected scores data (more than one subset may be given)\n",
    "    \"\"\"\n",
    "    if not isinstance(score_sets, list):\n",
    "        score_sets = [score_sets]\n",
    "    for scores in score_sets:\n",
    "        axs[0].plot(scores[:, 0], scores[:, 1], '.')\n",
    "        axs[1].plot(scores[:, 0], scores[:, 2], '.')\n",
    "        axs[2].plot(scores[:, 1], scores[:, 2], '.')\n",
    "    axs[0].set_xlabel('PC1')\n",
    "    axs[0].set_ylabel('PC2')\n",
    "    axs[1].set_xlabel('PC1')\n",
    "    axs[1].set_ylabel('PC3')\n",
    "    axs[2].set_xlabel('PC2')\n",
    "    axs[2].set_ylabel('PC3')\n",
    "    axs[1].set_title(title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(12, 4))\n",
    "result = plot_distribution(axs, scores, 'original data')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Compare alternative clustering methods.\n",
    "\n",
    "Here you look at the [single](https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.single.html#scipy.cluster.hierarchy.single), [complete](https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.complete.html#scipy.cluster.hierarchy.complete) and [ward](https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.ward.html#scipy.cluster.hierarchy.ward) methods in scipy, and the [Rodriguez-Laio](http://science.sciencemag.org/content/344/6191/1492) approach.\n",
    "\n",
    "#### Step 4.0: Create the distance matrix\n",
    "All clustering methods used here require the pre-calculation of the distances between all samples, i.e. the creation of a 2D distance matrix D where D[i, j] is the distance between samples i and j. Because D[i, j] = D[j, i] and D[i, i] = 0, its possible to save memory and calculation effort by just creating a \"condensed\" distance matrix: the upper triangle of the full 2D matrix. The `scipy` clustering methods can use a condensed distance matrix directly, the `rlcluster` method needs the full 2D version, but that's easy as `scipy` includes a utility function `squareform` that can convert between the two formats.\n",
    "\n",
    "For structure data like molecular conformations, RMSD makes a good distance metric - but others are possible. `MDTraj` includes a very fast routine that can calculate the RMSD between one conformation in a trajectory and a set of others, but to create the full (or condensed) 2D matrix we have to create our own routine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmsdlist = []\n",
    "for i in range(t.n_frames - 1):\n",
    "    rmsdlist.append(mdt.rmsd(t[i+1:], t[i]))\n",
    "rmsds = np.concatenate(rmsdlist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting array `rmsds` is a \"condensed\" distance matrix.\n",
    "\n",
    "#### Step 4.1: Look at RL first:\n",
    "\n",
    "The RL cluster method works on the full 2D distance matrix and really has only one adjustable parameter, `sigma`, which is related to the threshold 'distance' above which cluster merging stops. The default value of `sigma` is 5.0, this is often a good guess but typically (and this will go for other agglomerative clustering methods as well, as you will see later) there is an iterative process of trying a value, assessing the performance, and then adusting it if required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl_result = rlcluster.cluster(squareform(rmsds), sigma=2.0) # a value of 2.0 works here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RL method uses a \"decision graph\" to help assess whether the choice of sigma is appropriate (see the paper). The Python package includes a metod to generate this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rlcluster.decision_graph(rl_result, plt.axes())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that sigma=2.0 was a good choice. There seems a very clear separation (the dotted line, related to the chosen sigma value) between a few cluster centres and the rest of the data. The label assigned to each sample (the index of the cluster that it has been assigned to) is in the `.assignments` attribute of `rl_result`. We will use this to separate the clusters in the scores matrix. A value of 0 is used to mark outliers, so we skip over that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rl_labels = rl_result.assignments\n",
    "rl_clusters = [scores[rl_labels == i] for i in set(rl_labels) if i != 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4.2:  `Single` clustering:\n",
    "\n",
    "The `scipy` methods all work the same way. They take a distance matrix (in condensed form) and return a \"linkage matrix\" - see the scipy documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.linkage.html#scipy.cluster.hierarchy.linkage) for an explanation of the format of this matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zsingle = single(rmsds)\n",
    "plt.plot(Zsingle[:, 2], 'o')\n",
    "plt.plot([0.0, 10000], [0.027, 0.027], '--')\n",
    "plt.xlabel('partition number')\n",
    "plt.ylabel('distance threshold')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scipy` methods don't include a \"decision graph\" type tool, but a plot of selected data out of the linkage matrix Z does a similar job (see above). A cut at 0.027 looks reasonable, so we can then use this with the `fcluster` routine to produce a vector of labels for each sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_labels = fcluster(Zsingle, 0.027, criterion='distance')\n",
    "single_clusters = [scores[single_labels == i] for i in set(single_labels) if i != 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4.3:  `complete` clustering:\n",
    "\n",
    "The process is the same: create the linkage matrix, graph data from it, decide on the optimal clustering threshold from examination of the graph, and use `fcluster` to generate the labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zcomplete = complete(rmsds)\n",
    "plt.plot(Zcomplete[:, 2], 'o')\n",
    "plt.plot([0.0, 10000], [0.13, 0.13], '--') # we see - below - that this seems reasonable\n",
    "plt.xlabel('partition number')\n",
    "plt.ylabel('distance threshold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_labels = fcluster(Zcomplete, 0.13, criterion='distance')\n",
    "complete_clusters = [scores[complete_labels == i] for i in set(complete_labels) if i != 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4.4:  `ward` clustering:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zward = ward(rmsds)\n",
    "plt.plot(Zward[:, 2], 'o')\n",
    "plt.plot([0.0, 10000], [2.0, 2.0], '--') # this turns out to look reasonable\n",
    "plt.xlabel('partition number')\n",
    "plt.ylabel('distance threshold')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ward_labels = fcluster(Zward, 2.0, criterion='distance')\n",
    "ward_clusters = [scores[ward_labels == i] for i in set(ward_labels) if i != 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 5: Compare the results\n",
    "\n",
    "Sort each list of clusters so the biggest one comes first (so the colour coding of the plots is consistent), then plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ward_clusters.sort(key=lambda l: -len(l))\n",
    "rl_clusters.sort(key=lambda l: -len(l))\n",
    "single_clusters.sort(key=lambda l: -len(l))\n",
    "complete_clusters.sort(key=lambda l: -len(l))\n",
    "\n",
    "fig, axs = plt.subplots(4, 3, figsize=(12, 12))\n",
    "plot_distribution(axs[0], rl_clusters, 'RL')\n",
    "plot_distribution(axs[1], single_clusters, 'Single')\n",
    "plot_distribution(axs[2], complete_clusters, 'Complete')\n",
    "plot_distribution(axs[3], ward_clusters, 'Ward')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Single is obviously bad, complete looks unlikely. Ward and RL look reasonable and quite similar, but are they \"right\"?\n",
    "\n",
    "--------\n",
    "### Step 6: Evaluation of the different methods\n",
    "We have been doing all the work so far in Cartesian space, but the conformation of the alanine dipeptide is much more simply described by the values of the peptide bond phi and psi angles - there are no other significant degrees of freedom. Let's look at the data in the form of the familiar Ramachandran map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phipsi_atoms = [[1, 3, 4, 6], [3, 4, 6, 7]] # indices of the atoms defining phi and psi\n",
    "phipsi = mdt.compute_dihedrals(t, phipsi_atoms)\n",
    "\n",
    "Z, x_edges, y_edges = np.histogram2d(phipsi[:, 0], phipsi[:, 1], bins=60)\n",
    "plt.imshow(-np.log(Z.T), extent=[x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]], origin='lower')\n",
    "plt.xlabel('phi (radians)')\n",
    "plt.ylabel('psi (radians)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can identify cluster centres by eye in this map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "centres = np.array([[-1.2, -0.5], [-2.6, -0.5], [-1.4, 3], [1, 3], [-1.4, -3], [1, -3]])\n",
    "plt.imshow(-np.log(Z.T), extent=[x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]], origin='lower')\n",
    "plt.plot(centres[:, 0], centres[:, 1], 'ro')\n",
    "plt.xlabel('phi')\n",
    "plt.ylabel('psi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6.1: Use the Ramachandran map to generate a reference clustering.\n",
    "Assign each snapshot to the centre its phi and psi values are closest to. We use scipy's fast `cdist` function for this. The last two centres in our list are really the same as the two previous ones - just split by the -twopi/+twopi boundary in psi, so we reassign all points with label=5 to have label=3 and label=6 to habe label=4. Then we generate our list of clusters and sort by size, as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adist = cdist(phipsi, centres)\n",
    "\n",
    "phipsi_labels = adist.argmin(axis=1) + 1\n",
    "phipsi_labels = np.where(phipsi_labels==5, 3, phipsi_labels)\n",
    "phipsi_labels = np.where(phipsi_labels==6, 4, phipsi_labels)\n",
    "\n",
    "phipsi_clusters = [scores[phipsi_labels == i] for i in set(phipsi_labels) if i != 0]\n",
    "phipsi_clusters.sort(key=lambda l: -len(l))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 6.2: Compare the reference with Ward and RL results.\n",
    "Now we can compare the clustering done in torsional space according to the Ramachandran map with the Ward and RL clusterings done in Cartesian space (for this system, a much tougher challenge):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(3, 3, figsize=(12, 12))\n",
    "plot_distribution(axs[0], phipsi_clusters, 'PhiPsi')\n",
    "plot_distribution(axs[1], ward_clusters, 'Ward')\n",
    "plot_distribution(axs[2], rl_clusters, 'RL')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that both Ward and RL have done a good job, but the RL method seems to edge it.\n",
    "\n",
    "--------------\n",
    "\n",
    "### Sampling and clustering\n",
    "\n",
    "It will be very evident that the performance of a clustering method is going to be sensitive to the amount of data available. Let's have a bit of a look at that.\n",
    "\n",
    "We will take smaller and smaller subsets of the full trajectory, and see how that affects clustering using the RL approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "divisors = [1, 2, 4, 8, 16]\n",
    "\n",
    "#subrmsds = [squareform(rmsds)[::d][:,::d] for d in divisors]\n",
    "rmsdsq = squareform(rmsds)\n",
    "subrmsds = [rmsdsq[::d, ::d] for d in divisors]\n",
    "\n",
    "all_rl_clusters = []\n",
    "for i, d in enumerate(divisors):\n",
    "    result = rlcluster.cluster(subrmsds[i], sigma=2.0)\n",
    "    labels = result.assignments\n",
    "    clusters = [scores[::d][labels == i] for i in set(labels) if i != 0]\n",
    "    clusters.sort(key=lambda l: -len(l))\n",
    "    all_rl_clusters.append(clusters)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nd = len(divisors)\n",
    "fig, axs = plt.subplots(nd, 3, figsize=(12, 4*nd))\n",
    "for i, clusters in enumerate(all_rl_clusters):\n",
    "    plot_distribution(axs[i], clusters, 'sample size= {}'.format(10000/divisors[i]))\n",
    "    print('sample size: {:6.0f}  number of clusters: {}'.format(10000/divisors[i], len(clusters)))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary:\n",
    "\n",
    "Clustering is tricky, and subjective. Some key points:\n",
    "\n",
    "1. Compare different methods.\n",
    "2. There is always a bit of subjectivity - be transparent and always show how you decided on, for example, what clustering distance threshold to use.\n",
    "3. Project your data into the \"right\" space and use an appropriate distance metric, and your life will be easier. Here we saw that distributions and clustering that were quite tricky when attempted in Cartesian space with RMSD as the distance metric became almost trivial if torsion angle space was used instead. Unfortunately in a more complex system this can be much harder to do; the good news is that even with a generic RMSD-based approach, methods like Ward and Rodriguez-Laio can work well.\n",
    "4. No clustering method will work well if there isn't enough data. Conversely, when adding more data doesn't change the clustering, that may be useful evidence of convergence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Follow up work\n",
    "\n",
    "This notebook has been presented more as a demonstration/tutorial than as an interactive workshop, but there are plenty of ways you may wish to extend it:\n",
    "\n",
    "1. Explore what happens if you make different decisions about the cutoffs used to control the cluster merging process.\n",
    "2. Explore alternative clustering algorithms - most of the scipy methods have a very similar calling signature and you can 'swap out' methods here for alternative ones quite easily.\n",
    "3. Explore your own data - the Ramachandran analysis based section is very specific to the alanine dipeptide test case used here, but for the other sections of the notebook it should be straightforward to substitute in your own trajectory data and analyse that.\n",
    "\n",
    "## Footnote: fast 2D RMSD calculations\n",
    "\n",
    "The RMSD calculation routines in MDTraj are very fast, but as you will have noticed when you ran this step, it's still a time-consuming process, and it scales as the square of the number of samples so can eventually become prohibitive.\n",
    "\n",
    "If you can accept a small degree of approximation in the RMSD calculation (and most clustering applications will be fine with this), there is a PCA method that can be several times faster. Its based on the fact that in the full 3N-6 dimensional PCA space, the Cartesian distance between points is very close to their RMSD, multiplied by the square root of the number of atoms. Numpy contains a function `pdist` that can calculate Cartesian distance matrices very fast, so you can get an approximate 2D RMSD matrix this way.\n",
    "\n",
    "In the code below we test the PCA approach agains the MDTRaj approach, checking accuracy and speed, for larger and larger data sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmsd2d_mdt(t):\n",
    "    \"\"\"\n",
    "    Returns the condensed 2D RMSD matrix for the structures in trajectory t\n",
    "    \"\"\"\n",
    "    rmsdlist = []\n",
    "    for i in range(t.n_frames - 1):\n",
    "        rmsdlist.append(mdt.rmsd(t[i+1:], t[i]))\n",
    "    rmsds = np.concatenate(rmsdlist)\n",
    "    return rmsds\n",
    "\n",
    "def rmsd2d_pca(t):\n",
    "    \"\"\"\n",
    "    Returns an approximated condensed 2D distance matrix for the structures in trajectory t\n",
    "    \"\"\"\n",
    "    p = PCA()\n",
    "    scores = p.fit_transform(t.xyz)\n",
    "    scores /= np.sqrt(t.n_atoms)\n",
    "    return pdist(scores)\n",
    "\n",
    "def mue(d1, d2):\n",
    "    \"\"\"\n",
    "    The mean unsigned error between two sets of distance/rmsd calculations\n",
    "    \"\"\"\n",
    "    return np.abs(d1 -d2).mean()\n",
    "\n",
    "\n",
    "for n_frames in range(1000, 11000, 1000):\n",
    "    t_truncated = t[:n_frames]\n",
    "    start = time.time()\n",
    "    rmsd_mdt = rmsd2d_mdt(t_truncated)\n",
    "    mdt_time = time.time() - start\n",
    "\n",
    "    start = time.time()\n",
    "    rmsd_pca = rmsd2d_pca(t_truncated)\n",
    "    pca_time = time.time() - start\n",
    "\n",
    "    print('n_frames:{:6d} mdt time: {:5.2f} pca time: {:5.2f} mue: {:7.5f}'.format(\n",
    "           n_frames, mdt_time, pca_time, mue(rmsd_mdt, rmsd_pca)))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The performance you observe will be very data set and machine-dependent, but you might see the PCA method being up to ten times faster, for all trajectory sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}