03_parallelization_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Parallelization in GRASS GIS\n",
    "The goal of parallelization is to speed up computation by using multiple cores. This notebooks introduces parallelization concepts, existing parallelized tools, and approaches to parallelizing user scripts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start GRASS to run examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import GRASS packages\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "session = gj.init(\"~/data/grassdata\", \"nc_basic_spm_grass7\", \"user1\")\n",
    "\n",
    "# Set computational region to the elevation raster.\n",
    "gs.run_command(\"g.region\", raster=\"elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: most examples assume we are already in an active GRASS session."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Using already parallelized tools\n",
    "\n",
    "There are many tools in GRASS GIS that are already parallelized, [see the list](https://grass.osgeo.org/grass-stable/manuals/keywords.html#parallel). Many tools in GRASS Addons are parallelized as well.\n",
    "\n",
    "Generally, there are two types of implementation in GRASS GIS.\n",
    "Multithreading in C tools:\n",
    "   * Threads have low overhead, so they can be spawned more efficiently.\n",
    "   * Tools use OpenMP API. One of the advantages of OpenMP for software distribution is that code works (compiles and runs in serial) also without OpenMP library present on the system.\n",
    "   * Memory is shared, so programmer needs to be cautious about race conditions (e.g., writing into the same variable).\n",
    "   \n",
    "Multiprocessing in Python tools:\n",
    "   * There are multiple ways to implement it, typically tools use `subprocess` and `multiprocessing` package.\n",
    "   * Python tools are often wrappers around GRASS tools implemented in C. For example, tool [r.sun.daily](https://grass.osgeo.org/grass-stable/manuals/addons/r.sun.daily.html) runs [r.sun](https://grass.osgeo.org/grass-stable/manuals/r.sun.html) for multiple days in parallel.\n",
    "   \n",
    "Parallelized tools have `nprocs` parameter to specify number of cores to use. For C tools using OpenMP, GRASS GIS needs to be compiled with OpenMP support to take advantage of it. Both implementations work well on a single machine, but can't be scaled to a distributed system. Scaling to a distributed system is covered at the end of this tutorial.\n",
    "   \n",
    "Example of calling a parallelized tool in Bash and comparing the time with using 1 and 4 cores: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "time r.neighbors --q input=elevation output=elevation_smooth size=25 method=average nprocs=1\n",
    "time r.neighbors --q input=elevation output=elevation_smooth size=25 method=average nprocs=4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The speedup (processing time with 1 core / processing time with N cores) typically does not increase linearly with  the number of cores and parallel efficiency (speedup / N cores) decreases when adding cores. See, e.g.,  [benchmarks for r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html#performance). This behavior is due to the serial parts of the code (see [Amdahl's law](https://en.wikipedia.org/wiki/Amdahl%27s_law)) and computation overhead. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parallelization of workflows\n",
    "In a geoprocessing workflow, there are often multiple independent tasks that can be executed in parallel.\n",
    "The strategy how to divide the workflow into parallel tasks generally falls under either data-based or task-based parallelization.\n",
    "Task-based parallelism partitions tasks so that independent tasks can be completed simultaneously.\n",
    "With data-based parallelization, the spatial domain is partitioned for concurrent computations of individual spatial units \n",
    "and once processed, spatial units are mosaicked back into the initial spatial domain (if applicable)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data-based parallelization\n",
    "Data-based parallelism involves spatial domain decomposition, a process of splitting data into smaller datasets that can be processed in parallel.\n",
    "As part of [GRASS GIS Python API](https://grass.osgeo.org/grass-stable/manuals/libpython/index.html), [GridModule](https://grass.osgeo.org/grass-stable/manuals/libpython/pygrass.modules.grid.html) decomposes input data into rectangular tiles, executes a given tool for each tile in parallel, and merges the results. Effectively, tiling is applied virtually (using computational region), determining the spatial extent of analysis for each parallel process. In some cases such as moving window analysis, tiles need to overlap to get correct results. Note that GridModule can be fairly inefficient due to the overhead from merging results back and is therefore best suited for computationally-itensive tasks such as interpolation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example shows IDW interpolation split into 4 tiles. In this case, specifying an overlap is needed to get correct results without edge artifacts. Here, the number and size of tiles is automatically derived from the number of cores, but can be specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region vector=elev_points res=1 -p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%python\n",
    "import time\n",
    "from grass.pygrass.modules.grid import GridModule\n",
    "import time\n",
    "\n",
    "start = time.time()\n",
    "grid = GridModule(\n",
    "    \"v.surf.idw\",\n",
    "    input=\"elev_points\",\n",
    "    output=\"idw\",\n",
    "    column=\"value\",\n",
    "    npoints=25,\n",
    "    processes=4,\n",
    "    overlap=5,\n",
    "    quiet=True,\n",
    ")\n",
    "grid.run()\n",
    "print(f\"Elapsed time: {time.time() - start} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the same tool ran in serial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "time v.surf.idw --q input=elev_points column=value output=idw npoints=25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are tools that already integrate tiling. For example, addon [r.mapcalc.tiled](https://grass.osgeo.org/grass-stable/manuals/addons/r.mapcalc.tiled.html) uses the tiling concept for  raster algebra computation. More complex algebra expression will increase the speedup of this method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.extension r.mapcalc.tiled\n",
    "g.region raster=elevation res=1\n",
    "# parallel\n",
    "time r.mapcalc.tiled expression=\"rescaled_elevation = graph(elevation,60,1,80,10,100,100,120,100,140,1000,160,1000)\" nprocs=4\n",
    "# serial\n",
    "time r.mapcalc expression=\"rescaled_elevation = graph(elevation,60,1,80,10,100,100,120,100,140,1000,160,1000)\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task-based parallelization\n",
    "With task-based parallelism, we identify independent tasks and execute them concurrently.\n",
    "Tasks are typically GRASS processing tools executed as separate processes. Processes, unlike threads, do not share memory. When tasks are limited by disk I/O, parallel processing may have large overhead.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examples in Python\n",
    "There are multiple ways to execute tasks in parallel using Python, for example, there are libraries `multiprocessing` and `concurrent.futures`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following example viewsheds from different coordinates are computed in parallel using `multiprocessing.Pool` class. To avoid issues when using multiprocessing from Jupyter Notebook (multiprocessing.Pool does not work with interactive interpreters), we will first write a Python script with main function and then execute it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile example.py\n",
    "import os\n",
    "from multiprocessing import Pool\n",
    "import grass.script as gs\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    gs.run_command(\"r.viewshed\", input=\"elevation\", output=f\"viewshed_{cat}\", coordinates=(x, y))\n",
    "    return f\"viewshed_{cat}\"\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    viewpoints = [(633709, 225663, 1),\n",
    "                  (639432, 222826, 2),\n",
    "                  (640385, 220502, 3),\n",
    "                  (636521, 219353, 4)]\n",
    "    with Pool(processes=4) as pool:\n",
    "        maps = pool.map(viewshed, viewpoints)\n",
    "    print(maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region raster=elevation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run example.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examples in Bash\n",
    "In a simplest case, they can be executed in parallel from a command line shell by running a geoprocessing task in the background (by appending `&`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.region vector=schools,firestations res=30\n",
    "v.kernel --q input=schools output=kernel_schools radius=10000 &\n",
    "v.kernel --q input=firestations output=kernel_firestations radius=10000 &\n",
    "wait"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Larger number of tasks can be scheduled to run in parallel by tools such as [GNU Parallel](https://www.gnu.org/software/parallel/) and xargs.\n",
    "In this simple example, we use a loop to write commands into a file and execute those commands in parallel, using 2 cores. \n",
    "Whenever a task is finished, a next one is picked from the queued tasks.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for VECTOR in schools firestations hospitals\n",
    "do\n",
    "    echo v.kernel --q input=${VECTOR} output=kernel_${VECTOR} radius=10000 >> commands.sh\n",
    "done\n",
    "parallel -j 2 < commands.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See manual pages of GNU Parallel or xargs for more advanced uses. GNU Parallel can be configured to distribute jobs across multiple machines. In that case, use `--exec` interface described below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Safe execution of parallel tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While you can execute tasks in parallel within a single mapset, it is *not safe* when your tasks:\n",
    " \n",
    " * write output maps/files with identical names (common mistake, but easy to fix)\n",
    " * modify computational region\n",
    " * modify vector attribute database\n",
    " * modify MASK\n",
    " * use [r.reclass](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) to reclassify from the same base map\n",
    "\n",
    "The following sections provide solutions, starting with the option to execute tools in separate mapsets, which addresses all of the issues above, except r.reclass."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Executing processes in separate mapsets\n",
    "To execute tasks in separate mapsets, we can use `--exec` [interface](https://grass.osgeo.org/grass-stable/manuals/grass.html)\n",
    "that allows GRASS tools and user scripts to be executed in a GRASS GIS non-interactive session.\n",
    "This also enables parallelization on distributed systems.\n",
    "\n",
    "For example, here is a simple call to list all available vectors in PERMANENT mapset:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass ~/data/grassdata/nc_basic_spm_grass7/PERMANENT --exec g.list type=vector mapset=PERMANENT -t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the previous examples that was running within GRASS session in a single mapset can be rewritten so that each task runs in a newly created mapset. Note that by default newly created mapsets use default computational region for that GRASS location (you can use `g.region -s` to modify it). For raster computations, you need to change the computational region for each new mapset if the default one is not desired."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "for VECTOR in schools firestations hospitals\n",
    "do\n",
    "    # first create a new mapset with -c flag and set computational region based on the input vector\n",
    "    grass -c ~/data/grassdata/nc_basic_spm_grass7/${VECTOR} --exec g.region vector=${VECTOR} res=30\n",
    "    # write the command executing v.kernel in the newly created mapset to a file\n",
    "    echo grass ~/data/grassdata/nc_basic_spm_grass7/${VECTOR} --exec v.kernel --q input=${VECTOR} output=kernel_${VECTOR} radius=10000 >> exec_commands.sh\n",
    "done\n",
    "parallel -j 2 < exec_commands.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some cases, only a temporary mapset or location is needed, see [examples](https://grass.osgeo.org/grass-stable/manuals/grass.html#batch-jobs-with-the-exec-interface).\n",
    "Besides individual tools, the `--exec` interface can execute an entire script to enable more complex workflows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Safely modifying computational region in a single mapset\n",
    "\n",
    "Sometimes modifying computational region in a script is needed. It is a good practice to not change the global computational region, which effectively modifies a file in a mapset,\n",
    "but only change the environment variable `GRASS_REGION`.\n",
    "Here, we modified the previous viewshed example to compute in parallel viewsheds with different extents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile example.py\n",
    "import os\n",
    "from multiprocessing import Pool\n",
    "import grass.script as gs\n",
    "\n",
    "def viewshed(point):\n",
    "    x, y, cat = point\n",
    "    # copy current environment, modify and pass it to r.viewshed\n",
    "    env = os.environ.copy()\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(e=x + 300, w=x - 300, n=y + 300, s=y - 300, align=\"elevation\")\n",
    "    gs.run_command(\"r.viewshed\", input=\"elevation\", output=f\"viewshed_{cat}\", coordinates=(x, y), max_distance=300, env=env)\n",
    "    return f\"viewshed_{cat}\"\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    viewpoints = [(633709, 225663, 1),\n",
    "                  (639432, 222826, 2),\n",
    "                  (640385, 220502, 3),\n",
    "                  (636521, 219353, 4)]\n",
    "    with Pool(processes=4) as pool:\n",
    "        maps = pool.map(viewshed, viewpoints)\n",
    "    print(maps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run example.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Safely modifying vectors with attributes in a single mapset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default vector maps share a single SQLite database file, however SQLite does not support concurrent write access. That poses a problem when modifying vectors with attributes in parallel. While this can be solved by running the computations in separate mapsets, it is also possible to change the default behavior to write attributes of each vector to the vector's individual SQLite file. This behavior can be activated after a new mapset is created with:\n",
    "\n",
    "```\n",
    " db.connect driver=sqlite database='$GISDBASE/$LOCATION_NAME/$MAPSET/vector/$MAP/sqlite.db'\n",
    "```\n",
    "\n",
    "Alternatively, a PostgreSQL or another database backend can be used for attributes to offload the parallel writing to the database system."
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}