02-Variation1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variation 1: Handling an Unknown Amount of Data \n",
    "\n",
    "**Table of contents:**\n",
    " - [Problem description](#Problem-description)\n",
    " - [Passing arguments to our program](#Passing-arguments-to-our-program)\n",
    " - [Creating sample paths](#Creating-sample-paths)\n",
    " - [Writing HDF5 files](#Writing-HDF5-files)\n",
    " - [Discussion](#Discussion)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem description\n",
    "\n",
    "In the original problem formulation, we assumed that via program arguments we have the exact number of paths and time steps, and that we can generate the full data set in memory before writing it to storage. In many applications, these assumptions are not practical. We often don't know the total amount of data to be written, and we may not have enough memory to buffer it. In fact, we might simultaneously acquire data across different channels to be written to different HDF5 datasets. In circumstances like these, we must write data incrementally, which is what we will explore in this section of the tutorial.\n",
    "\n",
    "In terms of our model problem, we will drop the assumption that the number of time steps is constant for all paths. Let's instead assume that it is random, and that, rather than creating all `path_count` sample paths in one go, we can only produce them in batches of up to `batch_size` paths.\n",
    "\n",
    "Our revised problem formulation is thus: Store the floating-point values of `path_count` series of random length, produced in batches not exceeding `batch_size` series. In addition, store the following floating-point calibrations, which are the same for all series\n",
    "- the time step `dt`\n",
    "- the long-term process mean `mu`\n",
    "- the reversion rate to the mean `theta`\n",
    "- the volatility of the process `sigma`.\n",
    "\n",
    "**Approach:** Clearly, this data set does no longer conveniently fit into a rectilinear two-dimensional array of floating-point numbers. Instead, a two-dimensional __[ragged (or jagged) array](https://en.wikipedia.org/wiki/Jagged_array#:~:text=In%20computer%20science%2C%20a%20jagged,edges%20when%20visualized%20as%20output.)__ comes to mind. That would make a lot of sense in a situation where we would modify individual arrays (rows) with some frequency. For this discussion, we assume that we are dealing with a WORM (write once read many) scenario, in other words, we do not intend to update (or never change the length of) individual rows. With that assumption, we could store the data points of all paths in a batch in a contiguous one-dimensional array. But there is a problem: It would be like text without periods (`.`), i.e., we couldn't tell where a sentence starts or ends. To solve this problem, we maintain a second one-dimensional array in which we record the positions of the periods, or, equivalently, the positions of the first character of the first word of each sentence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Passing arguments to our program\n",
    "\n",
    "We need two helper functions to pass and parse the relevant arguments to our programs.\n",
    "\n",
    "We use Pranav Srinivas Kumar's [`argparse` module](https://github.com/p-ranav/argparse), which is the \"gold standard\" in the C++ world."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/parse_arguments1.hpp\n",
    "#ifndef PARSE_ARGUMENTS1_HPP\n",
    "#define PARSE_ARGUMENTS1_HPP\n",
    "\n",
    "#include \"argparse.hpp\"\n",
    "\n",
    "// Sets the options for which we are looking\n",
    "extern void set_options1(argparse::ArgumentParser& program);\n",
    "\n",
    "// Tests the options and retrieves the arguments\n",
    "extern int get_arguments1\n",
    "(\n",
    "    const argparse::ArgumentParser& program,\n",
    "    size_t&                         path_count,\n",
    "    size_t&                         batch_size,\n",
    "    double&                         dt,\n",
    "    double&                         theta,\n",
    "    double&                         mu,\n",
    "    double&                         sigma\n",
    ");\n",
    "\n",
    "#endif"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementation of our functions is rather unsurprising and speaks to the quality (usability) of `argparse`.\n",
    "\n",
    "In `set_options1()`, we inform `argparse` about the program options we want to support. For each argument, we tell `argparse` the short and long forms of the command line argument, a synopsis from which a help message will be created, and the default type and value of the argument. This information guides how `argparse` processes the `main` function's arguments, `argc` and `argv`.\n",
    "\n",
    "In `get_arguments1()`, we check if an option was specificed and what value was provided, or use the default value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/parse_arguments1.cpp\n",
    "#include \"parse_arguments1.hpp\"\n",
    "#include <cfloat>\n",
    "#include <iostream>\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "void set_options1(argparse::ArgumentParser& program)\n",
    "{\n",
    "    program.add_argument(\"-p\", \"--paths\")  // command line arguments\n",
    "    .help(\"chooses the numnber of paths\")  // synopsis\n",
    "    .default_value(size_t{100})            // default value\n",
    "    .scan<'u', size_t>();                  // expected type\n",
    "\n",
    "    program.add_argument(\"-b\", \"--batch\")\n",
    "    .help(\"chooses the batch size\")\n",
    "    .default_value(size_t{100})\n",
    "    .scan<'u', size_t>();\n",
    "\n",
    "    program.add_argument(\"-d\", \"--dt\")\n",
    "    .help(\"chooses the time step\")\n",
    "    .default_value(double{0.01})\n",
    "    .scan<'f', double>();\n",
    "\n",
    "    program.add_argument(\"-t\", \"--theta\")\n",
    "    .help(\"chooses the rate of reversion to the mean\")\n",
    "    .default_value(double{1.0})\n",
    "    .scan<'f', double>();\n",
    "\n",
    "    program.add_argument(\"-m\", \"--mu\")\n",
    "    .help(\"chooses the long-term mean of the process\")\n",
    "    .default_value(double{0.0})\n",
    "    .scan<'f', double>();\n",
    "\n",
    "    program.add_argument(\"-g\", \"--sigma\")\n",
    "    .help(\"chooses the volatility of the process\")\n",
    "    .default_value(double{0.1})\n",
    "    .scan<'f', double>();\n",
    "}\n",
    "\n",
    "int get_arguments1\n",
    "(\n",
    "    const argparse::ArgumentParser& program,\n",
    "    size_t&                         path_count,\n",
    "    size_t&                         batch_size,\n",
    "    double&                         dt,\n",
    "    double&                         theta,\n",
    "    double&                         mu,\n",
    "    double&                         sigma\n",
    ")\n",
    "{\n",
    "    path_count = program.get<size_t>(\"--paths\");\n",
    "    if (path_count == 0) {\n",
    "        cerr << \"Number of paths must be greater than zero\" << endl;\n",
    "        return -1;\n",
    "    }\n",
    "    batch_size = program.get<size_t>(\"--batch\");\n",
    "    if (batch_size == 0) {\n",
    "        cerr << \"Batch size must be greater than zero\" << endl;\n",
    "        return -1;\n",
    "    }\n",
    "    dt = program.get<double>(\"--dt\");\n",
    "    if (dt < DBL_MIN) {\n",
    "        cerr << \"Time step must be greater than zero\" << endl;\n",
    "        return -1;\n",
    "    }\n",
    "    theta = program.get<double>(\"--theta\");\n",
    "    if (theta < DBL_MIN) {\n",
    "        cerr << \"Reversion rate must be greater than zero\" << endl;\n",
    "        return -1;\n",
    "    }\n",
    "    mu = program.get<double>(\"--mu\");\n",
    "    sigma = program.get<double>(\"--sigma\");\n",
    "    if (sigma < DBL_MIN) {\n",
    "        cerr << \"Volatility must be greater than zero\" << endl;\n",
    "        return -1;\n",
    "    }\n",
    "\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating sample paths\n",
    "\n",
    "The sampler function needs only minor modifications. Instead of `path_count`, we pass the `batch_size` as a parameter. The only real change is the random generation of the path length (or step count) for each path. This occurs in line 26 of the implementation of `ou_sampler1()` where we sample a uniform distribution of integers between 1 and 65,535."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/ou_sampler1.hpp\n",
    "#ifndef OU_SAMPLER1_HPP\n",
    "#define OU_SAMPLER1_HPP\n",
    "\n",
    "#include \"hdf5.h\"\n",
    "#include <vector>\n",
    "\n",
    "// Creates `batch_size` sample paths of random length with parameters\n",
    "// `dt`, `theta`, `mu`, and `sigma`\n",
    "extern void ou_sampler1\n",
    "(\n",
    "    std::vector<double>&  ou_process,\n",
    "    std::vector<hsize_t>& offset,\n",
    "    const size_t&         batch_size,\n",
    "    const double&         dt,\n",
    "    const double&         theta,\n",
    "    const double&         mu,\n",
    "    const double&         sigma\n",
    ");\n",
    "\n",
    "#endif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/ou_sampler1.cpp\n",
    "\n",
    "#include \"ou_sampler1.hpp\"\n",
    "#include <random>\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "void ou_sampler1\n",
    "(\n",
    "    vector<double>&  ou_process,\n",
    "    vector<hsize_t>& offset,\n",
    "    const size_t&    batch_size,\n",
    "    const double&    dt,\n",
    "    const double&    theta,\n",
    "    const double&    mu,\n",
    "    const double&    sigma\n",
    ")\n",
    "{\n",
    "    // Store sample paths in one contiguous buffer\n",
    "    ou_process.clear();\n",
    "    offset.clear();\n",
    "    offset.push_back(0);\n",
    "\n",
    "    random_device rd;\n",
    "    mt19937 generator(rd());\n",
    "    uniform_int_distribution<unsigned int> path_len_dist(1, USHRT_MAX);  // path length is between 1 and 65535\n",
    "    normal_distribution<double> dist(0.0, sqrt(dt));  // N(0, dt)\n",
    "\n",
    "    size_t pos = 0;  // offset into the ou_process vector\n",
    "\n",
    "    // Generate a batch of paths and offsets\n",
    "    for (size_t i = 0; i < batch_size; ++i)\n",
    "    {\n",
    "        // Generate random path length\n",
    "        size_t step_count = path_len_dist(generator);\n",
    "        // Resize the vector to make room for the new path\n",
    "        ou_process.resize(pos + step_count);  \n",
    "            \n",
    "        // Generate the path\n",
    "        ou_process[pos] = 0; // start at x = 0\n",
    "        for (size_t j = 1; j < step_count; ++j)\n",
    "        {\n",
    "            auto dW = dist(generator);\n",
    "            ++pos;  // advance the offset\n",
    "            ou_process[pos] = ou_process[pos - 1] + theta * (mu - ou_process[pos - 1]) * dt + sigma * dW;\n",
    "        }\n",
    "\n",
    "        // This is the offset of the next path\n",
    "        offset.push_back(offset.back() + (hsize_t)step_count);  \n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing HDF5 files\n",
    "\n",
    "As mentioned in the [problem formulation](#problem-description) we use two one-dimensional arrays to store the sample paths in memory. The first array stores the path data points (`X`) in a batch and the second array stores the offsets of the first data point in a path. We map both arrays into one-dimensional HDF5 datasets in the file, `/paths/data` and `/paths/descr`, respectively.\n",
    "\n",
    "Notice that after processing the command line arguments, we know the number of sample paths to be generated, in other words, we know the extent of the dataset `/paths/descr`, which will be `path_count`. However, because we will create the sample paths in batches, we will write the offsets in batches as well. This will be the first example of selections and partial I/O.\n",
    "\n",
    "We do *not* know the final extent of the dataset `/paths/data`, because the length of an individual path is selected at random. Hence, in the creation of the dataset `/paths/data`, we must convey to the HDF5 library that we are dealing with a dataset of indefinite extent. We do this on lines 31-38. The initial extent of the dataspace is 0 elements, and the maximal extend is unlimited (`H5S_UNLIMITED`). On the storage side, we tell the library to grow the dataset, as needed, in increments or chunks of 131,072 elements. (See line 34.) This could be a function of the batch size, but we chose not to do that here.\n",
    "\n",
    "The loop starting on line 56 controls the creation of batches of sample paths. After populating the `ou_process` and `offset` arrays, we write the corresponding datasets `/paths/data` and `/paths/descr`. The main difference is that before writing to `/paths/data`, we must extend the dataset by `ou_process.size()` elements. (See lines 66-74.)\n",
    "\n",
    "In both cases, before calling `H5Dwrite`, we select the regions in the datasets that will receive the array elements using so-called hyperslab selections. A hyperslab selection is a regular pattern in a multi-dimnsional rectilinear grid that can be described by four parameters: start, stride, count, and block. __[NumPy slices](https://numpy.org/doc/stable/user/basics.indexing.html)__ are a special cases of hyperslabs.\n",
    "\n",
    "In the figure below, two examples of how flexible this space/selection-driven I/O can be. On the left, the pattern of incrementally appending batches of elements is shown. On the right, a more complex scenario is shown. Firstly, the ranks (dimensionality) and extents of the source and destination dataspaces are different. Secondly, only the subregions colored in black are written or read.\n",
    "\n",
    "<img src=\"img/HDF5_spaces.png\" title=\"Mapping between HDF5 dataspaces.\" />\n",
    "\n",
    "Finally, as in the model setting, we decorate the `/paths` group with the four attributes (lines 102-115)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/ou_hdf5.1.cpp\n",
    "#include \"parse_arguments1.hpp\"\n",
    "#include \"ou_sampler1.hpp\"\n",
    "\n",
    "#include \"hdf5.h\"\n",
    "#include <vector>\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "int main(int argc, char *argv[])\n",
    "{\n",
    "    size_t path_count, batch_size;\n",
    "    double dt, theta, mu, sigma;\n",
    "\n",
    "    argparse::ArgumentParser program(\"ou_hdf5.1\");\n",
    "    set_options1(program);\n",
    "    program.parse_args(argc, argv);\n",
    "    get_arguments1(program, path_count, batch_size, dt, theta, mu, sigma);\n",
    "\n",
    "    cout << \"Running with parameters:\"\n",
    "         << \" paths=\" << path_count << \" batch=\" << batch_size\n",
    "         << \" dt=\" << dt << \" theta=\" << theta << \" mu=\" << mu << \" sigma=\" << sigma << endl;\n",
    "\n",
    "    auto file = H5Fcreate(\"ou_process.1.h5\", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);\n",
    "    hid_t paths, descr;\n",
    "\n",
    "    auto lcpl = H5Pcreate(H5P_LINK_CREATE);\n",
    "    H5Pset_create_intermediate_group(lcpl, 1);\n",
    "\n",
    "    { // create the extendible `paths/data` dataset\n",
    "        hsize_t dimsf[] = {0, H5S_UNLIMITED};\n",
    "        auto space = H5Screate_simple(1, dimsf, &dimsf[1]);\n",
    "        auto dcpl = H5Pcreate(H5P_DATASET_CREATE);\n",
    "        hsize_t cdims[] = {128 * 1024};\n",
    "        H5Pset_chunk(dcpl, 1, cdims);\n",
    "        paths = H5Dcreate(file, \"/paths/data\", H5T_NATIVE_DOUBLE, space, lcpl, dcpl, H5P_DEFAULT);\n",
    "        H5Pclose(dcpl);\n",
    "        H5Sclose(space);\n",
    "    }\n",
    "    { // create the fixed size descriptors (= offsets into paths dataset) dataset `paths/descr`\n",
    "        hsize_t dimsf[] = { (hsize_t) path_count };\n",
    "        auto space = H5Screate_simple(1, dimsf, NULL);\n",
    "        descr = H5Dcreate(file, \"/paths/descr\", H5T_NATIVE_HSIZE, space, lcpl, H5P_DEFAULT, H5P_DEFAULT);\n",
    "        H5Sclose(space);    \n",
    "    }\n",
    "    \n",
    "    H5Pclose(lcpl);\n",
    "\n",
    "    // vectors to store the paths and the descriptors in a batch\n",
    "    vector<double> ou_process;\n",
    "    vector<hsize_t> offset;\n",
    "\n",
    "    // track the global (=across batches) offset\n",
    "    hsize_t global_pos = 0;  // \n",
    "\n",
    "    for (size_t p = 0; p < path_count; p += batch_size)\n",
    "    {\n",
    "        if (p + batch_size > path_count)  // last batch\n",
    "            batch_size = path_count - p;\n",
    "        cout << \"Generating paths \" << p << \" to \" << p + batch_size << endl;\n",
    "        \n",
    "        // Generate a batch of paths and offsets\n",
    "        ou_sampler1(ou_process, offset, batch_size, dt, theta, mu, sigma);\n",
    "\n",
    "        { // write the paths\n",
    "            auto path_space = H5Dget_space(paths);\n",
    "            hsize_t path_elem_count = (hsize_t) H5Sget_simple_extent_npoints(path_space);  // 1D dataset\n",
    "            H5Sclose(path_space);\n",
    "            hsize_t path_dims[] = {path_elem_count + (hsize_t) ou_process.size()};\n",
    "\n",
    "            // make room for more data\n",
    "            H5Dset_extent(paths, path_dims);   \n",
    "            // get the updated dataspace to make the correct selection(!)\n",
    "            path_space = H5Dget_space(paths);\n",
    "\n",
    "            hsize_t start[] = {path_elem_count};\n",
    "            hsize_t count[] = {(hsize_t) ou_process.size()};\n",
    "            H5Sselect_hyperslab(path_space, H5S_SELECT_SET, start, NULL, count, NULL);\n",
    "            auto mem_space = H5Screate_simple(1, count, NULL);\n",
    "            H5Sselect_all(mem_space);  // we want to write the whole vector\n",
    "            H5Dwrite(paths, H5T_NATIVE_DOUBLE, mem_space, path_space, H5P_DEFAULT, ou_process.data());\n",
    "            H5Sclose(mem_space);\n",
    "            H5Sclose(path_space);\n",
    "        }\n",
    "        { // write the path descriptors\n",
    "            auto descr_space = H5Dget_space(descr);\n",
    "            hsize_t start[] = { (hsize_t) p };\n",
    "            hsize_t count[] = { (hsize_t) offset.size()-1};  // offset has one extra element\n",
    "            H5Sselect_hyperslab(descr_space, H5S_SELECT_SET, start, NULL, count, NULL);\n",
    "            auto mem_space = H5Screate_simple(1, count, NULL);\n",
    "            H5Sselect_all(mem_space);\n",
    "            // the offsets are 0-based and we must correct this for the global offset\n",
    "            std::for_each(offset.begin(), offset.end(), [&](hsize_t &n){ n+=global_pos; });\n",
    "            global_pos = offset.back();\n",
    "            H5Dwrite(descr, H5T_NATIVE_HSIZE, mem_space, descr_space, H5P_DEFAULT, offset.data());\n",
    "            H5Sclose(mem_space);\n",
    "            H5Sclose(descr_space);\n",
    "        }\n",
    "    }\n",
    "    \n",
    "    { // make the file self-describing by adding a few attributes to `paths`\n",
    "        auto scalar = H5Screate(H5S_SCALAR);\n",
    "        auto acpl = H5Pcreate(H5P_ATTRIBUTE_CREATE);\n",
    "        H5Pset_char_encoding(acpl, H5T_CSET_UTF8);\n",
    "        auto set_attribute = [&](const string& name, const double& value) {\n",
    "            auto attr = H5Acreate_by_name(file, \"paths\", name.c_str(), H5T_NATIVE_DOUBLE, scalar, acpl, H5P_DEFAULT, H5P_DEFAULT);\n",
    "            H5Awrite(attr, H5T_NATIVE_DOUBLE, &value);\n",
    "            H5Aclose(attr);\n",
    "        };\n",
    "        set_attribute(\"dt\", dt);\n",
    "        set_attribute(\"θ\", theta);\n",
    "        set_attribute(\"μ\", mu);\n",
    "        set_attribute(\"σ\", sigma);\n",
    "        H5Pclose(acpl);\n",
    "        H5Sclose(scalar);\n",
    "    }\n",
    "\n",
    "    H5Dclose(descr);\n",
    "    H5Dclose(paths);\n",
    "    H5Fclose(file);\n",
    "\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g++ -std=c++17 -Wall -pedantic -I/usr/include/hdf5/serial -L/usr/lib/x86_64-linux-gnu -I./include  ./src/ou_hdf5.1.cpp ./src/parse_arguments1.cpp ./src/ou_sampler1.cpp -o ./build/ou_hdf5.1 -lhdf5_serial\n",
    "./build/ou_hdf5.1 -p 256\n",
    "ls -iks ou_process.1.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "In this section, we showed how to accomodate an unknown amount of data through incremental dataset extension and partial I/O. Because of the WORM assumption, we were able to replace the dense two-dimensional `data` array with two dense one-dimensional arrays. More general solutions are possible, but require careful considerations of all assumptions and specific performance goals.\n",
    "\n",
    "The H5Web plugin is not suitable to visualize individual paths, because we must first consult `/paths/descr` to identify the slice of `/paths/data` that represents a particular sample path. However, the visualization of `/paths/descr` is a useful spot check for correctness. It is necessary (but not sufficient) that the line graph shows a strictly monotonically increasing curve.\n",
    "\n",
    "On the other hand, it is easy to modify the Python script to pick out an visualize a specific sample path as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "f = h5py.File(\"ou_process.1.h5\")\n",
    "data = f[\"/paths/data\"]\n",
    "descr = f[\"/paths/descr\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_no = 42\n",
    "arr = data[descr[path_no]:descr[path_no+1]]\n",
    "print(f\"min: {arr.min():.2f}, max: {arr.max():.2f}, mean: {arr.mean():.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('_mpl-gallery')\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(np.arange(0,len(arr)), arr, linewidth=2.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hdf5-tutorial",
   "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.12.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}