chainreader comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Demonstrate discrepancy with trjcat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T14:07:12.653190Z",
     "start_time": "2018-06-30T14:07:12.321308Z"
    }
   },
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import numpy as np\n",
    "import MDAnalysis as mda\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T14:07:12.664342Z",
     "start_time": "2018-06-30T14:07:12.654470Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.18.1-dev'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mda.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T14:07:12.715016Z",
     "start_time": "2018-06-30T14:07:12.666080Z"
    }
   },
   "outputs": [],
   "source": [
    "def build_trajectories(*sequences, printout=False):\n",
    "    \"\"\"\n",
    "    Expose trjcat behavior for a given scenario.\n",
    "    \n",
    "    A scenario is given as a series of time sequences. The result is\n",
    "    returned as a list of time and origin of the frame. Each element of that\n",
    "    result list correspond to a frame of the conatenated trajectory, The\n",
    "    first element of each tuple is the time it corresponds to, the second\n",
    "    element is the index of the input the frame comes from.\n",
    "    \"\"\"\n",
    "    template = 'trjcat_test{}.xtc'\n",
    "    final_name = 'trjcat_test_cat.xtc'\n",
    "    \n",
    "    # Use an empty universe to have a topology\n",
    "    utop = mda.Universe.empty(1, trajectory=True)\n",
    "    \n",
    "    # Create synthetic trajectories. The times come from the user input,\n",
    "    # the coordinates indicate the index of the input.\n",
    "    for index, subseq in enumerate(sequences):\n",
    "        coords = np.zeros((len(subseq), 1, 3), dtype=np.float32) + index\n",
    "        u = mda.Universe(utop._topology, coords)\n",
    "        out_traj = mda.coordinates.XTC.XTCWriter(template.format(index),\n",
    "                                                 n_atoms=len(u.atoms))\n",
    "        with out_traj:\n",
    "            for ts, time in zip(u.trajectory, subseq):\n",
    "                # The time step needs a box to avoid a warning\n",
    "                ts.dimensions = [10, 10, 10, 90, 90, 90]\n",
    "                # The time is set from the user input\n",
    "                ts.time = time\n",
    "                out_traj.write(ts)\n",
    "    \n",
    "    flist = [template.format(index) for index, _ in enumerate(sequences)]\n",
    "    subprocess.call(['gmx', 'trjcat', '-o', final_name, '-f'] + flist)\n",
    "    u = mda.Universe(utop._topology, final_name)\n",
    "    final_seq = [(ts.time, int(ts.positions[0, 0])) for ts in u.trajectory]\n",
    "    \n",
    "    ucontinuous = mda.Universe(utop._topology, flist, continuous=True)\n",
    "    continuous_seq = [(ts.time, int(ts.positions[0, 0])) for ts in ucontinuous.trajectory]\n",
    "    \n",
    "    # Clean the files\n",
    "    del u  # Maybe we can try to close the file before deleting it?\n",
    "    for index, _ in enumerate(sequences):\n",
    "        os.remove(template.format(index))\n",
    "    os.remove(final_name)\n",
    "    # Clean the offset to avoid a warning\n",
    "    os.remove('.{}_offsets.npz'.format(final_name))\n",
    "    \n",
    "    return final_seq == continuous_seq, final_seq, continuous_seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T14:07:12.726357Z",
     "start_time": "2018-06-30T14:07:12.717028Z"
    }
   },
   "outputs": [],
   "source": [
    "test_ranges = [range(0, 6), range(2, 5), range(2, 5), range(2, 5), range(3, 8)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-30T14:07:12.800246Z",
     "start_time": "2018-06-30T14:07:12.728622Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jon/dev/mdanalysis/package/MDAnalysis/coordinates/XDR.py:195: UserWarning: Reload offsets from trajectory\n",
      " ctime or size or n_atoms did not match\n",
      "  warnings.warn(\"Reload offsets from trajectory\\n \"\n",
      "/home/jon/dev/mdanalysis/package/MDAnalysis/coordinates/chain.py:312: UserWarning: Missing frame in continuous chain\n",
      "  warnings.warn(\"Missing frame in continuous chain\", UserWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(False,\n",
       " [(0.0, 0),\n",
       "  (1.0, 0),\n",
       "  (2.0, 3),\n",
       "  (3.0, 4),\n",
       "  (4.0, 4),\n",
       "  (5.0, 4),\n",
       "  (6.0, 4),\n",
       "  (7.0, 4)],\n",
       " [(0.0, 0),\n",
       "  (1.0, 0),\n",
       "  (2.0, 1),\n",
       "  (2.0, 2),\n",
       "  (2.0, 3),\n",
       "  (3.0, 4),\n",
       "  (4.0, 4),\n",
       "  (5.0, 4),\n",
       "  (6.0, 4),\n",
       "  (7.0, 4)])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "build_trajectories(*[list(r) for r in test_ranges])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "apricot (py3.6.2)",
   "language": "python",
   "name": "apricot"
  },
  "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.3"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}