Lecture-5-HDF5-h5py-PyTables.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('./..')\n",
    "import data as tutorial_data\n",
    "file_name = tutorial_data.get_file('BYxRM')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to HDF5 in Python\n",
    "\n",
    "based on [Giacomo Debidda's notebook](https://github.com/jackdbd/hdf5-pydata-munich/blob/master/hdf5_in_python.ipynb)\n",
    "\n",
    "\n",
    "## HDF5: a filesystem in a file\n",
    "\n",
    "**HDF** stands for *Hierarchical Data Format* and is a data model, library, and file format for storing and managing big and complex data.\n",
    "\n",
    "An HDF5 file can be thought of as a container (or group) that holds a variety of heterogeneous data objects (or datasets). The datasets can be almost anything: images, tables, graphs, or even documents, such as PDF or Excel.\n",
    "\n",
    "- Datasets (i.e. files in a filesystem)\n",
    "- Groups (i.e. directories in a filesystem)\n",
    "- Attributes (i.e. metadata of file/directory)\n",
    "\n",
    "\n",
    "![hdf5_structure](images/hdf5_structure.jpg \"HDF5 structure\")\n",
    "\n",
    "Every object in an HDF5 file has a name, and they’re arranged in a POSIX-style hierarchy with `/`-separators.  \n",
    "The “folders” in this system are called groups. The **File** object we created is itself a *group*, in this case the *root* group, named `/`  \n",
    "\n",
    "**Groups work like dictionaries, and datasets work like NumPy arrays**\n",
    "\n",
    "**/** root group (every HDF5 file has a root group)\n",
    "\n",
    "**/foo** member of the root group called foo\n",
    "\n",
    "**/foo/bar** member of the group foo called bar\n",
    "\n",
    "## HDF5 in the Python data stack\n",
    "\n",
    "![hdf5_in_python](images/h5py-pytables-refactor.png \"HDF5 in Python data stack\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## h5py\n",
    "\n",
    "The h5py package is a Pythonic interface to the HDF5 binary data format.\n",
    "\n",
    "- Thin, pythonic wrapper around the HDF5 C API\n",
    "- Written in Cython\n",
    "- Tries to expose the entire HDF5 C API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='w') as f:\n",
    "    f.create_dataset(name='my_dataset', data=[1.0, 2.7, 3.7, 4.5])\n",
    "    f.create_dataset(name='another_dataset', data=[1, 2, 3, 4])\n",
    "    f.create_dataset(name='yet_another_dataset', data=[1, 2, 3, 4], dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<HDF5 dataset \"my_dataset\": shape (4,), type \"<f8\">\n",
      "[1.  2.7 3.7 4.5]\n",
      "[1.  2.7 3.7 4.5]\n"
     ]
    }
   ],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='r') as f:\n",
    "    # the array is just a proxy object\n",
    "    print(f['my_dataset'])\n",
    "    # the actual data can be accessed with these 2 syntaxes\n",
    "    print(f['my_dataset'][:])\n",
    "    print(f['my_dataset'][...])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preallocation on disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='w') as f:\n",
    "    dataset = f.create_dataset(name='my_dataset', shape=(8, 1))\n",
    "    dataset[0] = 5.2\n",
    "    dataset[1] = 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='a') as f:\n",
    "    dataset = f['my_dataset']\n",
    "    dataset[2] = 3.9\n",
    "    dataset[3] = 8.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pick the correct HDF5 datatype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0   1 254 255   0 255 254]\n"
     ]
    }
   ],
   "source": [
    "arr = np.array([0, 1, 254, 255, 256, -1, -2], dtype='uint8')\n",
    "print(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0   1 254 255 255   0   0]\n"
     ]
    }
   ],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='w') as f:\n",
    "    f.create_dataset(name='my_dataset', shape=(7,), dtype=h5py.h5t.STD_U8BE)\n",
    "    f['my_dataset'][0:8] = [0, 1, 254, 255, 123456, -1, -2]\n",
    "    print(f[\"my_dataset\"][:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='w') as f:\n",
    "    f.create_group(name='group1')\n",
    "    group2 = f.create_group(name='group2')\n",
    "    group2.create_group(name='group3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<HDF5 group \"/group2\" (1 members)>\n"
     ]
    }
   ],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='r') as f:\n",
    "    group3 = f['group2/group3']\n",
    "    print(group3.parent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='w') as f:\n",
    "    dataset = f.create_dataset(name='my_dataset', data=[1, 2, 3, 4])\n",
    "    dataset.attrs['Unit'] = 4\n",
    "    gr = f.create_group(name='my_group')\n",
    "    gr.attrs['Created'] = '18/12/2017'\n",
    "    gr.attrs.create(name='Versions', data=np.array([123, 456]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Traverse a HDF5 file with h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "my_dataset\n",
      "my_group\n"
     ]
    }
   ],
   "source": [
    "with h5py.File(name='data/my_h5py_file.h5', mode='r') as f:\n",
    "    f.visit(print)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HDF5 Command Line Tools\n",
    "\n",
    "[Here](https://support.hdfgroup.org/products/hdf5_tools/#h5dist) you can find the command line tools developed by the HDF Group. You don't need h5py or PyTables to use them.\n",
    "\n",
    "If you are on Ubuntu, you can install them with `sudo apt install hdf5-tools` but they are also installed when you install `h5py` or `PyTables` with `conda`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/                        Group\r\n",
      "/my_dataset              Dataset {4}\r\n",
      "/my_group                Group\r\n"
     ]
    }
   ],
   "source": [
    "# -r stands for 'recursive' \n",
    "!h5ls -r 'data/my_h5py_file.h5'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HDF5 \"data/my_h5py_file.h5\" {\r\n",
      "GROUP \"/\" {\r\n",
      "   DATASET \"my_dataset\" {\r\n",
      "      DATATYPE  H5T_STD_I64LE\r\n",
      "      DATASPACE  SIMPLE { ( 4 ) / ( 4 ) }\r\n",
      "      DATA {\r\n",
      "      (0): 1, 2, 3, 4\r\n",
      "      }\r\n",
      "      ATTRIBUTE \"Unit\" {\r\n",
      "         DATATYPE  H5T_STD_I64LE\r\n",
      "         DATASPACE  SCALAR\r\n",
      "         DATA {\r\n",
      "         (0): 4\r\n",
      "         }\r\n",
      "      }\r\n",
      "   }\r\n",
      "   GROUP \"my_group\" {\r\n",
      "      ATTRIBUTE \"Created\" {\r\n",
      "         DATATYPE  H5T_STRING {\r\n",
      "            STRSIZE H5T_VARIABLE;\r\n",
      "            STRPAD H5T_STR_NULLTERM;\r\n",
      "            CSET H5T_CSET_UTF8;\r\n",
      "            CTYPE H5T_C_S1;\r\n",
      "         }\r\n",
      "         DATASPACE  SCALAR\r\n",
      "         DATA {\r\n",
      "         (0): \"18/12/2017\"\r\n",
      "         }\r\n",
      "      }\r\n",
      "      ATTRIBUTE \"Versions\" {\r\n",
      "         DATATYPE  H5T_STD_I64LE\r\n",
      "         DATASPACE  SIMPLE { ( 2 ) / ( 2 ) }\r\n",
      "         DATA {\r\n",
      "         (0): 123, 456\r\n",
      "         }\r\n",
      "      }\r\n",
      "   }\r\n",
      "}\r\n",
      "}\r\n"
     ]
    }
   ],
   "source": [
    "!h5dump 'data/my_h5py_file.h5'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading existing HDF5 file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = h5py.File(file_name, 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['genotype', 'phenotype']"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(f.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The root object contains two children (*genotype* and *phenotype*).  \n",
    "To find out what type (dataset or group) of object the *phenotype* is we can run: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<HDF5 group \"/phenotype\" (3 members)>"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f['phenotype']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the *phenotype* object itself is a group with 3 children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['col_header', 'matrix', 'row_header']"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(f['phenotype'].keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An easy way to iterate over all elements of a HDF5 file is to use the visit method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "genotype: <HDF5 group \"/genotype\" (3 members)> \n",
      "genotype/col_header: <HDF5 group \"/genotype/col_header\" (3 members)> \n",
      "genotype/col_header/alleles: <HDF5 dataset \"alleles\": shape (11623, 2), type \"|S37\"> \n",
      "genotype/col_header/chrom: <HDF5 dataset \"chrom\": shape (11623,), type \"<i8\"> \n",
      "genotype/col_header/pos: <HDF5 dataset \"pos\": shape (11623,), type \"<i8\"> \n",
      "genotype/matrix: <HDF5 dataset \"matrix\": shape (1008, 11623), type \"|u1\"> \n",
      "genotype/row_header: <HDF5 group \"/genotype/row_header\" (1 members)> \n",
      "genotype/row_header/sample_ID: <HDF5 dataset \"sample_ID\": shape (1008,), type \"|S99\"> \n",
      "phenotype: <HDF5 group \"/phenotype\" (3 members)> \n",
      "phenotype/col_header: <HDF5 group \"/phenotype/col_header\" (1 members)> \n",
      "phenotype/col_header/phenotype_ID: <HDF5 dataset \"phenotype_ID\": shape (46,), type \"|S22\"> \n",
      "phenotype/matrix: <HDF5 dataset \"matrix\": shape (1008, 46), type \"<f8\"> \n",
      "phenotype/row_header: <HDF5 group \"/phenotype/row_header\" (1 members)> \n",
      "phenotype/row_header/sample_ID: <HDF5 dataset \"sample_ID\": shape (1008,), type \"|S22\"> \n"
     ]
    }
   ],
   "source": [
    "def printname(name):\n",
    "    print(\"%s: %s \" % (name, f[name]))\n",
    "f.visit(printname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every object in an HDF5 file can have custom attributes which can be checked like any other dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "list(f['phenotype']['matrix'].attrs.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us examine the *phenotype/matrix* dataset as a `Dataset` object.   \n",
    "The object we obtained isn’t an array, but an *HDF5* dataset. Like *NumPy* arrays, datasets have both a *shape* and a *data type*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape: (1008,46)\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "dset = f['phenotype/matrix']\n",
    "print('shape: (%s,%s)' % dset.shape)\n",
    "print('dtype: %s' % dset.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They also support array-style slicing. This is how you read and write data from a dataset in the file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.24597261,  2.6897812 , -0.46048117,  0.03086018, -0.29291799,\n",
       "       -1.09117648,  0.70117087, -0.86986989, -0.13682744, -0.28441332])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# retrieve all rows between 10 and 20 for the the second column\n",
    "dset[10:20, 2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When filtering/slicing/indexing rows and columns of a dataset, not the entire dataset is loaded into the main memory. This is very useful if you need to retrieve a subset of the data from a terabyte big HDF5 file.  \n",
    "To read the entire dataset into the memory run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-7.32351971e+00,  2.79992827e-01,  3.13118166e-01, ...,\n",
       "         8.90841948e-01,  4.11837231e+00,  8.59281836e+00],\n",
       "       [-8.09823582e+00, -2.06326076e-01, -5.34843783e-01, ...,\n",
       "         6.06164117e-03,  6.65150029e-02, -4.22047646e+00],\n",
       "       [ 7.60571968e+00, -1.27959825e-01, -3.11102424e-01, ...,\n",
       "         1.72315729e+00,  5.71408803e+00, -6.50651895e+00],\n",
       "       ...,\n",
       "       [            nan,  3.21844466e-01,             nan, ...,\n",
       "         1.80142187e+00,             nan, -4.99069190e-01],\n",
       "       [            nan, -1.00150686e+00,             nan, ...,\n",
       "        -1.74752575e+00,             nan, -7.43596147e+00],\n",
       "       [            nan,             nan,             nan, ...,\n",
       "                    nan,             nan,             nan]])"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dset[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ndarray'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-7.32351971e+00,  2.79992827e-01,  3.13118166e-01, ...,\n",
       "         8.90841948e-01,  4.11837231e+00,  8.59281836e+00],\n",
       "       [-8.09823582e+00, -2.06326076e-01, -5.34843783e-01, ...,\n",
       "         6.06164117e-03,  6.65150029e-02, -4.22047646e+00],\n",
       "       [ 7.60571968e+00, -1.27959825e-01, -3.11102424e-01, ...,\n",
       "         1.72315729e+00,  5.71408803e+00, -6.50651895e+00],\n",
       "       ...,\n",
       "       [            nan,  3.21844466e-01,             nan, ...,\n",
       "         1.80142187e+00,             nan, -4.99069190e-01],\n",
       "       [            nan, -1.00150686e+00,             nan, ...,\n",
       "        -1.74752575e+00,             nan, -7.43596147e+00],\n",
       "       [            nan,             nan,             nan, ...,\n",
       "                    nan,             nan,             nan]])"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr = dset[:]\n",
    "print(type(arr))\n",
    "arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excercise: Create h5py file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "Write a function, <code>create_h5py_hdf5</code>, which takes an array of `participants` and their `reaction_times` as parameters and stores them in an HDF5 file using the h5py library. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "participants = np.load(\"data/experiment_participants.npy\")\n",
    "reaction_times = np.load(\"data/experiment_data.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_h5py_hdf5(participants, reaction_times):\n",
    "    \"\"\"Create a HDF5 file using h5py containing a dataset reaction_times and the participants as an attribute.\n",
    "       Make sure to:\n",
    "\n",
    "      * Store reaction_times as a datasets inside /experiments\n",
    "      * Datasets reaction_times should be compressed with zlib (5)\n",
    "      * Return the filename of the generated HDF5 file\n",
    "      * Use float32 as a data type for the reaction_times\n",
    "      * Store the mean reaction times as an attribute 'ameanvg' on the reaction_times dataset\n",
    "      * Sore the participants array as an attriubute \"participants\" on the /experiments group\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    participants : numpy array with shape (n,),\n",
    "    reaction_times : numpy array with shape (n,)\n",
    "\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    # YOUR CODE HERE\n",
    "    return file_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from numpy.testing import assert_array_equal, assert_almost_equal\n",
    "file_name = create_h5py_hdf5(participants,reaction_times)\n",
    "assert os.path.exists(file_name)\n",
    "with h5py.File(file_name, 'r') as hf:\n",
    "    assert len(hf.keys()) == 1\n",
    "    assert len(hf['experiments'].keys()) == 1\n",
    "    assert_array_equal(participants, np.char.decode(hf['experiments'].attrs['participants'][:]))\n",
    "    assert_almost_equal(reaction_times, hf['experiments']['reaction_times'][:],decimal=3)\n",
    "    assert hf['experiments']['reaction_times'].compression == 'gzip'\n",
    "    assert hf['experiments']['reaction_times'].compression_opts == 5\n",
    "    print(\"Success!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyTables\n",
    "\n",
    "- Higher abstraction over HDF5 (it's more \"battery included\")\n",
    "- Searches are faster than in h5py\n",
    "- Does not depend on h5py (at the moment)\n",
    "\n",
    "![pytables logo](images/pytables-logo.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tables as tb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    f.create_array(where='/', \n",
    "                   name='my_array',\n",
    "                   title='My PyTables Array',\n",
    "                   obj=[1, 2, 3, 4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PyTables has a feature called \"Natural Naming\": nodes (i.e. datasets and groups in the HDF5 file) can be accessed with the dot notation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/my_array (Array(4,)) 'My PyTables Array'\n"
     ]
    }
   ],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='r') as f:\n",
    "    print(f.root.my_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    f.create_group(where='/', name='my_group')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Attributes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    f.create_array(where=f.root, name='my_array', obj=[1, 2, 3, 4], title='My PyTables Array')\n",
    "    f.set_node_attr(where='/my_array', attrname='SomeAttribute', attrvalue='SomeValue')\n",
    "    f.create_group(where='/', name='my_group')\n",
    "    f.set_node_attr(where='/my_group', attrname='SomeOtherAttribute', attrvalue=123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HDF5 datasets have many abstractions in PyTables\n",
    "\n",
    "Homogenous dataset:\n",
    "\n",
    "- Array\n",
    "- CArray\n",
    "- EArray\n",
    "- VLArray\n",
    "\n",
    "Heterogenous dataset:\n",
    "\n",
    "- Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_rows = 1000000  # 1 million\n",
    "gaussian = np.random.normal(loc=0, scale=1, size=num_rows).astype('float32')\n",
    "uniform = np.random.uniform(low=100, high=150, size=num_rows).astype('uint8')\n",
    "num_columns = 5\n",
    "matrix = np.random.random((num_rows, num_columns)).astype('float32')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Array (again!)\n",
    "[Docs](http://www.pytables.org/usersguide/libref/homogenous_storage.html#the-array-class)\n",
    "\n",
    "- Fastest I/O speed\n",
    "- Must fit in memory\n",
    "- Not compressible\n",
    "- Not enlargeable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.26 ms, sys: 9.18 ms, total: 12.4 ms\n",
      "Wall time: 14.3 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    f.create_array(where='/', name='gaussian', obj=gaussian)\n",
    "    f.create_array(where='/', name='uniform', obj=uniform)\n",
    "    f.create_array(where='/', name='matrix', obj=matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r--  1 uemit.seren  staff    24M Nov 14 18:22 data/my_pytables_file.h5\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lh data/my_pytables_file.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CArray\n",
    "[Docs](http://www.pytables.org/usersguide/libref/homogenous_storage.html#carrayclassdescr)\n",
    "\n",
    "- Compressible\n",
    "- Not enlargeable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "filters = tb.Filters(complevel=5, complib='zlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tips on how to use compression (from the PyTables docs)\n",
    "\n",
    "- A mid-level (5) compression is sufficient. No need to go all the way up (9)\n",
    "- Use zlib if you must guarantee complete portability\n",
    "- Use blosc all other times (it is optimized for HDF5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 627 ms, sys: 22.9 ms, total: 650 ms\n",
      "Wall time: 652 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "#     f.create_carray(where='/', name='gaussian', obj=gaussian)\n",
    "#     f.create_carray(where='/', name='uniform', obj=uniform)\n",
    "#     f.create_carray(where='/', name='matrix', obj=matrix)\n",
    "    f.create_carray(where='/', name='gaussian', obj=gaussian, filters=filters)\n",
    "    f.create_carray(where='/', name='uniform', obj=uniform, filters=filters)\n",
    "    f.create_carray(where='/', name='matrix', obj=matrix, filters=filters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r--  1 uemit.seren  staff    20M Nov 14 18:22 data/my_pytables_file.h5\r\n"
     ]
    }
   ],
   "source": [
    "!ls -lh data/my_pytables_file.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EArray\n",
    "[Docs](http://www.pytables.org/usersguide/libref/homogenous_storage.html#earrayclassdescr)\n",
    "\n",
    "- Enlargeable on one dimension (append)\n",
    "- Compressible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One (and only one) of the shape dimensions *must* be 0.\n",
    "# The dimension being 0 means that the resulting EArray object can be extended along it.\n",
    "# Multiple enlargeable dimensions are not supported right now.\n",
    "shape = (0, num_columns)\n",
    "\n",
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    # you can create an EArray and fill it later, but you need to specify atom and shape\n",
    "    f.create_earray(where='/', name='my_earray', atom=tb.Float32Atom(), shape=shape, filters=filters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='a') as f:\n",
    "    earray = f.root.my_earray\n",
    "    earray.append(sequence=matrix[0:1, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='a') as f:\n",
    "    earray = f.root.my_earray\n",
    "    earray.append(sequence=matrix[1:5, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table\n",
    "[Docs](http://www.pytables.org/usersguide/libref/structured_storage.html?highlight=table#tableclassdescr)\n",
    "\n",
    "- Data can be heterogeneous (i.e. different shapes and different dtypes)\n",
    "- The structure of a table is declared by its description\n",
    "- multi-column searches\n",
    "\n",
    "In order to emulate in Python records mapped to HDF5 C structs PyTables implements a special class so as to easily define all its fields and other properties. It's called `IsDescription`.\n",
    "\n",
    "A *description* defines the table structure (basically, the *schema* of your table)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Particle(tb.IsDescription):\n",
    "    identity = tb.StringCol(itemsize=22, dflt=' ', pos=0)  # character String\n",
    "    idnumber = tb.Int16Col(dflt=1, pos=1)  # short integer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'identity': StringCol(itemsize=22, shape=(), dflt=b' ', pos=0), 'idnumber': Int16Col(shape=(), dflt=1, pos=1)}\n"
     ]
    }
   ],
   "source": [
    "print(Particle.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    table = f.create_table(where='/', name='my_table', description=Particle)\n",
    "\n",
    "    for i in range(100):\n",
    "        table.row['identity'] = 'I am {identity}'.format(identity=i)\n",
    "        table.row['idnumber'] = i\n",
    "        table.row.append()\n",
    "    table.flush()  # Flush the table buffers to release memory and make sure are written to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    table = f.create_table(where='/', name='my_table', description=Particle)\n",
    "    \n",
    "    for i in range(100):\n",
    "        table.append([('I am {identity}'.format(identity=i), i)])\n",
    "        \n",
    "    table.flush()  # Flush the table buffers to release memory and make sure are written to disk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Traverse a HDF5 file with PyTables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/my_table\n"
     ]
    }
   ],
   "source": [
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='r') as f:\n",
    "    for node in f.walk_nodes('/', classname='Table'):\n",
    "        print('{}'.format(node._v_pathname))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Expressions with NumExpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.uniform(low=1, high=5, size=num_rows).astype('float32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 23.7 ms, sys: 10.1 ms, total: 33.8 ms\n",
      "Wall time: 33 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    carray = f.create_carray(where='/', name='carray_without_numexpr', atom=tb.Float32Atom(), shape=x.shape)\n",
    "    carray[:] = x**3 + 0.5*x**2 - x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 13.6 ms, sys: 12.7 ms, total: 26.3 ms\n",
      "Wall time: 16 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "with tb.open_file(filename='data/my_pytables_file.h5', mode='w') as f:\n",
    "    carray = f.create_carray(where='/', name='carray_with_numexpr', atom=tb.Float32Atom(), shape=x.shape)\n",
    "    ex = tb.Expr('x**3 + 0.5*x**2 - x')\n",
    "    ex.set_output(carray) # output will got to the CArray on disk\n",
    "    ex.eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NYC Yellow Taxi Dataset (2015)\n",
    "\n",
    "On a regular laptop (Thinkpad X220 i5 10GB RAM) it took roughly:\n",
    "\n",
    " - **40 minutes** to read/store a **single CSV**\n",
    " - **8 hours** to read/store an entire **year**\n",
    " \n",
    " ![still-waiting-meme](images/still-waiting-meme.jpg \"Still waiting\")\n",
    " \n",
    " \n",
    " To see how to improve the runtime with PyTables check out [Giacomo Debidda's notebook](https://github.com/jackdbd/hdf5-pydata-munich/blob/master/hdf5_in_python.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excercise: Create PyTables file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "Write a function, <code>create_pytables_hdf5</code>, which takes an array of `participants` and their `reaction_times` as parameters and stores them in an HDF5 file as a PyTables table. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    You can use <a href=\"https://numpy.org/doc/stable/reference/generated/numpy.core.records.fromarrays.html\">Numpy Structured arrays</a> to store the data in the PyTable table. Also make sure that you calculate the mean reaction times for each participant (the axis is important when running the numpy's mean function)\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "participants = np.load(\"data/experiment_participants.npy\")\n",
    "reaction_times = np.load(\"data/experiment_data.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_pytables_hdf5(participants, reaction_times):\n",
    "    \"\"\"Create a HDF5 file using PyTables that contains a PyTable table.\n",
    "       Make sure to:\n",
    "\n",
    "      * table should be called 'experiment'\n",
    "      * First column are participants\n",
    "      * Second column are the mean reaction_time for each participant\n",
    "      * Return the filename of the generated HDF5 file\n",
    "      * Create an index on the participants column\n",
    "      * Reaction_times column should be of type float32\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    participants : numpy array with shape (n,),\n",
    "    reaction_times : numpy array with shape (n,)\n",
    "\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    # YOUR CODE HERE\n",
    "    return file_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from numpy.testing import assert_array_equal, assert_almost_equal\n",
    "file_name = create_pytables_hdf5(participants,reaction_times)\n",
    "assert os.path.exists(file_name)\n",
    "with tb.open_file(file_name, 'r') as hf:\n",
    "    assert len(hf.list_nodes('/')) == 1\n",
    "    assert  hf.root.experiment is not None\n",
    "    recarray = hf.root.experiment.read()\n",
    "    assert hf.root.experiment.cols.participant.is_indexed\n",
    "    assert hf.root.experiment.cols.mean_reaction_time.type == 'float32'\n",
    "    assert_array_equal(participants, np.char.decode(hf.root.experiment.read(field='participant')))\n",
    "    assert_almost_equal(np.mean(reaction_times,axis=1), hf.root.experiment.read(field='mean_reaction_time'),decimal=3)\n",
    "    print(\"Success!\")"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}