03-Variation2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variation 2: MPI Parallel HDF5\n",
    "\n",
    "**Table of contents:**\n",
    " - [Dividing the work](#Dividing-the-work)\n",
    " - [Passing arguments to our program](#Passing-arguments-to-our-program)\n",
    " - [Writing HDF5 files](#Writing-HDF5-files)\n",
    " - [Discussion](#Discussion)\n",
    "\n",
    "In this section of the tutorial, we learn how to use the __[Message Massing Interface (MPI)](https://en.wikipedia.org/wiki/Message_Passing_Interface)__ to parallelize the sample path generation, and how multiple processes can simultaneously write to a shared HDF5 file. I turns out that apart from MPI boilerplate and the code that partitions the workload, the original code remains largely unchanged.\n",
    "\n",
    "This is not an introduction to MPI. You can think of it as a framework that allows us to spawn multiple (many!) processes in a single machine or across a cluster of machines, which can then exchange messages process-to-process, or collectively among groups of processes. Given `nranks` $\\geq 1$ processes, each process is identified by its `rank`, where $0\\leq$ `rank` $<$ `nranks`.\n",
    "\n",
    "\n",
    "Our goal is to store the floating-point values of `path_count` series with `step_count` values each in a two-dimensional array of `path_count` rows and `step_count` columns. Instead of a single process, we have `nranks` processes at our disposal for the sample path generation. The first point of discussion is how we divide the work among `nranks` processes. Since passing arguments to the program(s) is identical to the original program, there isn't much to new to say about that. The only other novelty is to explore how multiple processes write different parts of a single dataset (`/data`). In the section dealing with an unknown amount of data, we saw selections and partial I/O as a means of incrementally writing a dataset. We will use a similar technique when writing in parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dividing the work\n",
    "\n",
    "Given the task of creating `path_count` sample paths (of length `step_count`), and having `nranks` processes to accomplish that task, it is reasonable to evenly divide the work among processes. Since `path_count` may not be divisible by `nranks`, there will be a few ranks to which we assign an extra path. The function `partition_work()`, shown below, implements this idea. Given a total of `path_count` paths and `nranks` processes, for each process `rank`, it computes the range of rows [`start`, `stop`) for which the process will generate sample paths. In other words, process `rank` will generate and write the row-block `[start,stop) x step_count` of the two-dimensional dataset `/data` with dimensions `path_count x step_count`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/partitioner.hpp\n",
    "\n",
    "#ifndef PARTITIONER_HPP\n",
    "#define PARTITIONER_HPP\n",
    "\n",
    "#include <cstddef>\n",
    "\n",
    "extern void partition_work(std::size_t path_count, int rank, int nranks, std::size_t &start, std::size_t &stop);\n",
    "\n",
    "#endif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/partitioner.cpp\n",
    "#include \"partitioner.hpp\"\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "void partition_work(size_t path_count, int rank, int nranks, size_t &start, size_t &stop)\n",
    "{\n",
    "  size_t count = path_count / nranks;\n",
    "  size_t remainder = path_count % nranks;\n",
    "\n",
    "  if ((size_t)rank < remainder) {\n",
    "    // The first 'remainder' ranks get 'count + 1' tasks each\n",
    "    start = rank * (count + 1);\n",
    "    stop = start + count;\n",
    "  } else {\n",
    "    // The remaining 'nranks - remainder' ranks get 'count' task each\n",
    "    start = rank * count + remainder;\n",
    "    stop = start + (count - 1);\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Passing arguments to our program\n",
    "\n",
    "By now, this should be a familiar theme. We need to account for an extra argument, `--use_subfiling`, that tells us whether to use [sub-filing](https://docs.hdfgroup.org/hdf5/rfc/RFC_VFD_subfiling_200424.pdf), i.e., we use multiple sub-files instead of a single shared file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/parse_arguments2.hpp\n",
    "#ifndef PARSE_ARGUMENTS2_HPP\n",
    "#define PARSE_ARGUMENTS2_HPP\n",
    "\n",
    "#include \"argparse.hpp\"\n",
    "\n",
    "// Sets the options for which we are looking\n",
    "extern void set_options2(argparse::ArgumentParser& program);\n",
    "\n",
    "// Tests the options and retrieves the arguments\n",
    "extern int get_arguments2\n",
    "(\n",
    "    const argparse::ArgumentParser& program,\n",
    "    size_t&                         path_count,\n",
    "    size_t&                         step_count,\n",
    "    double&                         dt,\n",
    "    double&                         theta,\n",
    "    double&                         mu,\n",
    "    double&                         sigma,\n",
    "    bool&                           use_subfiling\n",
    ");\n",
    "\n",
    "#endif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/parse_arguments2.cpp\n",
    "#include \"parse_arguments2.hpp\"\n",
    "#include <cfloat>\n",
    "#include <iostream>\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "void set_options2(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(\"-s\", \"--steps\")\n",
    "    .help(\"chooses the number of steps\")\n",
    "    .default_value(size_t{1000})\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",
    "    program.add_argument(\"--use_subfiling\");\n",
    "}\n",
    "\n",
    "int get_arguments2\n",
    "(\n",
    "    const argparse::ArgumentParser& program,\n",
    "    size_t&                         path_count,\n",
    "    size_t&                         step_count,\n",
    "    double&                         dt,\n",
    "    double&                         theta,\n",
    "    double&                         mu,\n",
    "    double&                         sigma,\n",
    "    bool&                           use_subfiling\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",
    "    step_count = program.get<size_t>(\"--steps\");\n",
    "    if (step_count == 0) {\n",
    "        cerr << \"Number of steps 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",
    "    use_subfiling = program.is_used(\"--use_subfiling\");\n",
    "\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing HDF5 files\n",
    "\n",
    "Lines 16-37 are just MPI boilerplate to initialize MPI and to determine the total number of MPI ranks, `nprocs`, and the rank of the current process, `myid`.\n",
    "Each process obtains its workpackage from the `partition_work()` function (line 72). We know already how to create sample paths and can reuse the existing `ou_sampler()` function. (Line 75)\n",
    "\n",
    "In lines 98-106, we construct the hyperslab selections (\"NumPy slices\"), in-memory and in-file, and are ready to call `H5Dwrite()` in line 112.\n",
    "\n",
    "The attribute decoration (lines 123-143) works exactly the same as in the original code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile src/ou_hdf5_mpi.cpp\n",
    "#include \"parse_arguments.hpp\"\n",
    "#include \"parse_arguments2.hpp\"\n",
    "#include \"partitioner.hpp\"\n",
    "#include \"docstring.hpp\"\n",
    "#include \"ou_sampler.hpp\"\n",
    "\n",
    "#include \"hdf5.h\"\n",
    "#include <iostream>\n",
    "#include <vector>\n",
    "\n",
    "using namespace std;\n",
    "\n",
    "int main(int argc, char *argv[])\n",
    "{\n",
    "    // <ALTERNATIVE>:: The standard MPI IO File Driver only requires:\n",
    "    //\n",
    "    // MPI_Init(&argc, &argv);\n",
    "    //\n",
    "    // However, the subfiling VFD requires MPI_Init_thread with MPI_THREAD_MULTIPLE\n",
    "\n",
    "    int mpi_thread_required = MPI_THREAD_MULTIPLE;\n",
    "    int mpi_thread_provided = 0;\n",
    "\n",
    "    MPI_Init_thread(&argc, &argv, mpi_thread_required, &mpi_thread_provided);\n",
    "    if (mpi_thread_provided < mpi_thread_required)\n",
    "    {\n",
    "      cout << \"MPI_THREAD_MULTIPLE not supported\" << endl;\n",
    "      MPI_Abort(MPI_COMM_WORLD, -1);\n",
    "    }\n",
    "\n",
    "    // Get the number of processes\n",
    "    int nprocs;\n",
    "    MPI_Comm_size(MPI_COMM_WORLD, &nprocs);\n",
    "    // Get the rank of the process\n",
    "    int myid;\n",
    "    MPI_Comm_rank(MPI_COMM_WORLD, &myid);\n",
    "\n",
    "    size_t path_count, step_count;\n",
    "    double dt, theta, mu, sigma;\n",
    "#ifdef H5_HAVE_SUBFILING_VFD    \n",
    "    bool subfiling;\n",
    "#endif    \n",
    "\n",
    "    argparse::ArgumentParser program(\"ou_hdf5_mpi\");\n",
    "#ifdef H5_HAVE_SUBFILING_VFD\n",
    "    set_options2(program);\n",
    "#else\n",
    "    set_options(program);\n",
    "#endif\n",
    "    program.parse_args(argc, argv);\n",
    "#ifdef H5_HAVE_SUBFILING_VFD\n",
    "    get_arguments2(program, path_count, step_count, dt, theta, mu, sigma, subfiling);\n",
    "#else\n",
    "    get_arguments(program, path_count, step_count, dt, theta, mu, sigma);\n",
    "#endif     \n",
    "\n",
    "    if (myid == 0)\n",
    "    {\n",
    "        cout << \"Running on \" << nprocs << \" MPI ranks with parameters:\"\n",
    "            << \" paths=\" << path_count << \" steps=\" << step_count\n",
    "            << \" dt=\" << dt << \" theta=\" << theta << \" mu=\" << mu << \" sigma=\" << sigma\n",
    "#ifdef H5_HAVE_SUBFILING_VFD\n",
    "            << \" subfiling=\" << subfiling\n",
    "#endif\n",
    "            << endl;\n",
    "    }\n",
    "\n",
    "    vector<double> ou_process;\n",
    "    \n",
    "    size_t start, stop;\n",
    "    partition_work(path_count, myid, nprocs, start, stop);\n",
    "    size_t my_path_count = stop - start + 1;\n",
    "\n",
    "    ou_sampler(ou_process, my_path_count, step_count, dt, theta, mu, sigma);\n",
    "    \n",
    "    // Use the Subfiling or MPI-IO driver\n",
    "    hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);\n",
    "#ifdef H5_HAVE_SUBFILING_VFD\n",
    "    if(subfiling)\n",
    "      H5Pset_fapl_subfiling(fapl, NULL);\n",
    "    else\n",
    "#endif\n",
    "      H5Pset_fapl_mpio(fapl, MPI_COMM_WORLD, MPI_INFO_NULL);\n",
    "\n",
    "    //\n",
    "    // Write the sample paths to an HDF5 file using the HDF5 C-API!\n",
    "    //\n",
    "    auto file = H5Fcreate(\"ou_process.2.h5\", H5F_ACC_TRUNC, H5P_DEFAULT, fapl);\n",
    "\n",
    "    add_docstring(file, \".\", \"source\", \"https://github.com/HDFGroup/hdf5-tutorial\");\n",
    "\n",
    "    { // create & write the dataset\n",
    "        hsize_t dimsf[] = {(hsize_t)path_count, (hsize_t)step_count};\n",
    "        auto filespace = H5Screate_simple(2, dimsf, NULL);\n",
    "        auto dataset = H5Dcreate(file, \"/dataset\", H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);\n",
    "        \n",
    "        // Define, by rank, a selection in memory and write it to a hyperslab in the file.\n",
    "        hsize_t count[]  = {my_path_count, step_count};\n",
    "        hsize_t offset[] = {start, 0};\n",
    "        \n",
    "        hid_t memspace = H5Screate_simple(2, count, NULL);\n",
    "\n",
    "        // Select hyperslab in the file.\n",
    "        // hid_t filespace = H5Dget_space(dataset);\n",
    "        H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);\n",
    "\n",
    "        // <OPTIONAL> Create property list for collective dataset write.\n",
    "        hid_t dxpl = H5Pcreate(H5P_DATASET_XFER);\n",
    "        H5Pset_dxpl_mpio(dxpl, H5FD_MPIO_COLLECTIVE);\n",
    "\n",
    "        H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, dxpl, ou_process.data());\n",
    "\n",
    "        // housekeeping\n",
    "        H5Pclose(dxpl);\n",
    "        H5Sclose(memspace);\n",
    "        H5Dclose(dataset);\n",
    "        H5Pclose(fapl);\n",
    "        H5Sclose(filespace);\n",
    "    }\n",
    "\n",
    "    { // make the file self-describing by adding a few attributes to `dataset`\n",
    "        add_docstring(file, \"dataset\", \"comment\", \"This dataset contains sample paths of an Ornstein-Uhlenbeck process.\");\n",
    "        add_docstring(file, \"dataset\", \"Wikipedia\", \"https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process\");\n",
    "        add_docstring(file, \"dataset\", \"rows\", \"path\");\n",
    "        add_docstring(file, \"dataset\", \"columns\", \"time\");\n",
    "\n",
    "        auto scalar = H5Screate(H5S_SCALAR);\n",
    "        auto acpl = H5Pcreate(H5P_ATTRIBUTE_CREATE);\n",
    "        H5Pset_char_encoding(acpl, H5T_CSET_UTF8);\n",
    "        \n",
    "        auto set_attribute = [&](const string& name, const double& value) {\n",
    "            auto attr = H5Acreate_by_name(file, \"dataset\", 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",
    "\n",
    "        H5Pclose(acpl);\n",
    "        H5Sclose(scalar);\n",
    "    }\n",
    "\n",
    "    H5Fclose(file);\n",
    "\n",
    "    MPI_Finalize();\n",
    "\n",
    "    return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "mpicxx -std=c++17 -Wall -pedantic -I/usr/include/hdf5/openmpi -L/usr/lib/x86_64-linux-gnu -I./include  ./src/ou_hdf5_mpi.cpp ./src/parse_arguments.cpp ./src/parse_arguments2.cpp ./src/docstring.cpp ./src/ou_sampler.cpp ./src/partitioner.cpp -o ./build/ou_hdf5_mpi -lhdf5_openmpi -lmpi\n",
    "mpiexec --host localhost:2 -n 2 ./build/ou_hdf5_mpi -p 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -iks ou_process.2.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "In this section, we saw a simple example of multiple MPI processes simutaneously writing to the same dataset in a shared HDF5 file. This divide-and-conquer pattern is common practice in High-Performance Computing to efficiently read and write extremely large datasets across many thousands of processes. "
   ]
  }
 ],
 "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
}