r.sim.water: Add a benchmark for the number of threads

python/grass/benchmark/__init__.py

@@ -11,6 +11,9 @@
 """Benchmarking for GRASS modules
 
 This subpackage of the grass package is experimental and the API can change anytime.
+When a benchmark of a tool needs more than the subpackage offers, for example an
+option of a plot, the subpackage is changed together with the benchmark, typically
+in the same pull request.
 The API of the package is defined by what is imported in the top-level ``__init__.py``
 file of the subpackage.
 

python/grass/benchmark/plots.py

@@ -32,7 +32,41 @@ def get_pyplot(to_file):
     return plt
 
 
-def nprocs_plot(results, filename=None, title=None, metric="time"):
+def nprocs_ranges(result, metric):
+    """Return the lowest and highest values of a metric for each nprocs.
+
+    The range of time is the fastest and the slowest run. The range of
+    speedup spans the ratios of every serial run to every run with the given
+    nprocs; efficiency is that divided by nprocs. The serial runs are those
+    with nprocs 1, or with the smallest nprocs when there is no 1. At that
+    nprocs, speedup and efficiency are 1 by definition, so the range is
+    empty there.
+
+    Returns None when the result has no *all_times*.
+    """
+    all_times = getattr(result, "all_times", None)
+    if not all_times:
+        return None
+    if metric == "time":
+        return [min(t) for t in all_times], [max(t) for t in all_times]
+    serial_index = result.nprocs.index(min(result.nprocs))
+    serial = all_times[serial_index]
+    lows, highs = [], []
+    for i, (nprocs, times) in enumerate(zip(result.nprocs, all_times, strict=True)):
+        if i == serial_index:
+            low = high = 1
+        else:
+            low, high = min(serial) / max(times), max(serial) / min(times)
+        if metric == "efficiency":
+            low, high = low / nprocs, high / nprocs
+        lows.append(low)
+        highs.append(high)
+    return lows, highs
+
+
+def nprocs_plot(
+    results, filename=None, title=None, metric="time", figsize=None, dpi=None
+):
     """Plot results from a multiple nprocs (thread) benchmarks.
 
     *results* is a list of individual results from separate benchmarks.
@@ -48,31 +82,39 @@ def nprocs_plot(results, filename=None, title=None, metric="time"):
 
     Optionally, result can have an *all_times* attribute which is a list
     of lists. One sublist is all times recorded for each value of nprocs.
+    With it, the range of the runs is shaded around each line in the line's
+    color, see :func:`nprocs_ranges`.
 
     Each result can come with a different list of nprocs, i.e., benchmarks
     which used different values for nprocs can be combined in one plot.
+
+    *figsize* (width and height in inches) and *dpi* are passed to
+    matplotlib; the size of a saved image in pixels is their product.
     """
     ylabel = ""
     plt = get_pyplot(to_file=bool(filename))
-    _, axes = plt.subplots()
+    _, axes = plt.subplots(figsize=figsize, dpi=dpi)
 
     x_ticks = set()  # gather x values
     for result in results:
         x = result.nprocs
         x_ticks.update(x)
         if metric == "time":
-            mins = [min(i) for i in result.all_times]
-            maxes = [max(i) for i in result.all_times]
-            plt.plot(x, result.times, label=result.label)
-            plt.fill_between(x, mins, maxes, color="gray", alpha=0.3)
+            values = result.times
             ylabel = "Time [s]"
         elif metric in {"speedup", "efficiency"}:
+            values = getattr(result, metric)
             ylabel = metric.title()
-            plt.plot(x, getattr(result, metric), label=result.label)
         else:
             msg = f"Invalid metric '{metric}' in result, it should be:\
                 'time', 'speedup' or 'efficiency'"
             raise ValueError(msg)
+        (line,) = plt.plot(x, values, label=result.label)
+        ranges = nprocs_ranges(result, metric)
+        if ranges:
+            plt.fill_between(
+                x, *ranges, color=line.get_color(), alpha=0.25, linewidth=0
+            )
     plt.legend()
     # If there is not many x values, show ticks for each, but use default
     # ticks when there is a lot of x values.
@@ -88,12 +130,13 @@ def nprocs_plot(results, filename=None, title=None, metric="time"):
     plt.ylabel(ylabel)
     if title:
         plt.title(title)
-    elif metric == "times":
+    elif metric == "time":
         plt.title("Execution time by processing elements")
     elif metric in {"speedup", "efficiency"}:
         plt.title(f"{metric.title()} by processing elements")
     if filename:
         plt.savefig(filename)
+        plt.close()
     else:
         plt.show()
 

python/grass/benchmark/tests/grass_benchmark_plots_test.py

@@ -0,0 +1,74 @@
+"""Tests of the nprocs plot and its ranges in grass.benchmark"""
+
+from types import SimpleNamespace
+
+import pytest
+
+from grass.benchmark import nprocs_plot
+from grass.benchmark.plots import nprocs_ranges
+
+
+def make_result(nprocs, all_times, label="test"):
+    times = [sum(t) / len(t) for t in all_times]
+    return SimpleNamespace(
+        nprocs=nprocs,
+        all_times=all_times,
+        times=times,
+        speedup=[times[0] / t for t in times],
+        efficiency=[times[0] / (n * t) for n, t in zip(nprocs, times, strict=True)],
+        label=label,
+    )
+
+
+def test_time_range_is_min_max_of_runs():
+    result = make_result([1, 2], [[10.0, 12.0, 11.0], [6.0, 5.0, 7.0]])
+    assert nprocs_ranges(result, "time") == ([10.0, 5.0], [12.0, 7.0])
+
+
+def test_speedup_range_spans_ratios_of_serial_to_parallel_runs():
+    result = make_result([1, 2], [[10.0, 12.0], [5.0, 6.0]])
+    lows, highs = nprocs_ranges(result, "speedup")
+    assert lows == [1, 10.0 / 6.0]
+    assert highs == [1, 12.0 / 5.0]
+
+
+def test_efficiency_range_is_speedup_range_per_thread():
+    result = make_result([1, 4], [[8.0, 12.0], [2.0, 4.0]])
+    lows, highs = nprocs_ranges(result, "efficiency")
+    assert lows == [1, 8.0 / 4.0 / 4]
+    assert highs == [1, 12.0 / 2.0 / 4]
+
+
+def test_serial_is_smallest_nprocs():
+    result = make_result([2, 4], [[6.0, 6.0], [3.0, 4.0]])
+    lows, highs = nprocs_ranges(result, "speedup")
+    assert lows == [1, 6.0 / 4.0]
+    assert highs == [1, 6.0 / 3.0]
+
+
+def test_no_range_without_all_times():
+    result = SimpleNamespace(nprocs=[1, 2], times=[2.0, 1.0], label="test")
+    assert nprocs_ranges(result, "time") is None
+
+
+@pytest.mark.parametrize("metric", ["time", "speedup", "efficiency"])
+def test_plot_file_has_requested_size(tmp_path, metric):
+    pytest.importorskip("matplotlib")
+    results = [make_result([1, 2, 4], [[4.0, 4.2], [2.1, 2.3], [1.2, 1.4]])]
+    filename = tmp_path / f"{metric}.png"
+    nprocs_plot(
+        results, filename=str(filename), metric=metric, figsize=(6, 4.5), dpi=50
+    )
+    with open(filename, "rb") as f:
+        header = f.read(24)
+    # PNG width and height are big-endian integers at bytes 16 to 24.
+    assert int.from_bytes(header[16:20], "big") == 300
+    assert int.from_bytes(header[20:24], "big") == 225
+
+
+def test_plot_without_all_times(tmp_path):
+    pytest.importorskip("matplotlib")
+    results = [SimpleNamespace(nprocs=[1, 2], times=[2.0, 1.0], label="test")]
+    filename = tmp_path / "time.png"
+    nprocs_plot(results, filename=str(filename))
+    assert filename.is_file()

raster/r.sim/r.sim.water/benchmark/benchmark_r_sim_water_nprocs.py

@@ -0,0 +1,115 @@
+"""Benchmarking of r.sim.water with different numbers of threads
+
+Runs r.sim.water with nprocs from 1 to MAX_NPROCS on the elevation raster
+of the SECREF North Carolina dataset (secref_northcarolina_usa_epsg6542)
+at several resolutions of the same extent. At resolutions coarser than
+the original 1 m, the elevation is the average of the 1 m cells. Otherwise
+r.sim.water runs with its defaults, so the number of walkers (twice the
+number of cells) and the number of iterations grow with finer resolution.
+
+Run in a temporary mapset of the dataset, for example:
+
+    grass --tmp-mapset ~/grassdata/secref_northcarolina_usa_epsg6542 \
+        --exec python benchmark_r_sim_water_nprocs.py
+
+The plots and the results as JSON are written to the current directory.
+To redo the plots from saved results without running the benchmark, pass
+the JSON file as an argument.
+
+The plots show the mean of the repeated runs with the range of the runs
+shaded (see grass.benchmark.nprocs_plot).
+
+@author Vaclav Petras
+"""
+
+import sys
+from subprocess import DEVNULL
+
+from grass.pygrass.modules import Module
+
+import grass.benchmark as bm
+import grass.script as gs
+
+# Users can modify the resolutions (in meters), the number of threads and
+# the number of repeats. With the default duration, the serial run at 1 m
+# takes about two minutes.
+RESOLUTIONS = [4, 2, 1]
+MAX_NPROCS = 16
+REPEAT = 3
+
+
+def main():
+    if len(sys.argv) > 1:
+        results = bm.load_results_from_file(sys.argv[1]).results
+    else:
+        results = []
+        for resolution in RESOLUTIONS:
+            benchmark(resolution, results)
+        bm.save_results_to_file(results, "r_sim_water_benchmark.json")
+    plot(results)
+
+
+def plot(results):
+    """Plot time, speedup and efficiency with the range of the runs"""
+    for metric in ["time", "speedup", "efficiency"]:
+        # Twice the 600 pixels of the documentation, shown at half the size.
+        bm.nprocs_plot(
+            results,
+            filename=f"r_sim_water_benchmark_{metric}.png",
+            title=f"r.sim.water {metric}",
+            metric=metric,
+            figsize=(6, 4.5),
+            dpi=200,
+        )
+
+
+def benchmark(resolution, results):
+    elevation = "benchmark_r_sim_water_elevation"
+    dx = "benchmark_r_sim_water_dx"
+    dy = "benchmark_r_sim_water_dy"
+    depth = "benchmark_r_sim_water_depth"
+
+    Module("g.region", raster="elevation")
+    Module("g.region", res=resolution, flags="a")
+    # Average the 1 m cells rather than taking the nearest one.
+    Module(
+        "r.resamp.stats",
+        input="elevation",
+        output=elevation,
+        method="average",
+        overwrite=True,
+    )
+    Module("r.slope.aspect", elevation=elevation, dx=dx, dy=dy, flags="e")
+    cells = gs.region()["cells"]
+
+    module = Module(
+        "r.sim.water",
+        elevation=elevation,
+        dx=dx,
+        dy=dy,
+        depth=depth,
+        random_seed=1,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+        quiet=True,
+    )
+    results.append(
+        bm.benchmark_nprocs(
+            module,
+            label=f"{cells / 1000:.0f}k cells ({resolution} m)",
+            max_nprocs=MAX_NPROCS,
+            repeat=REPEAT,
+        )
+    )
+    Module(
+        "g.remove",
+        quiet=True,
+        flags="f",
+        type="raster",
+        pattern="benchmark_r_sim_water_*",
+    )
+
+
+if __name__ == "__main__":
+    main()

raster/r.sim/r.sim.water/r.sim.water.html

@@ -332,6 +332,38 @@ <h3>Random numbers and parallel processing</h3>
 since the order in which the walkers reach a cell depends on the
 threads. Use <b>nprocs=1</b> when results must be reproducible.
 
+<h3>Performance</h3>
+
+To enable parallel processing, the user can specify the number of
+threads to be used with the <b>nprocs</b> parameter (default 1). Figures
+below show benchmark results for the <code>elevation</code> raster map of the
+SECREF North Carolina dataset at 4 m, 2 m and 1 m resolution with the
+default number of walkers, as the mean of 3 runs on Intel&reg; Xeon&reg; W-2295
+CPU @ 3.00GHz &times; 18. See the benchmark script in the source code for
+more details.
+
+<p>
+The time step is derived from the cell size and the mean flow velocity,
+so the number of iterations, and with it the run time, depends on the
+terrain as well as on the number of cells. As a result, the benchmark
+results may vary depending on the study area.
+
+<div align="center" style="margin: 10px">
+  <a href="r_sim_water_benchmark_time.png"><img src="r_sim_water_benchmark_time.png"
+  width="600" height="450"
+  alt="time for r.sim.water with different numbers of cells" border="0"></a>
+  <a href="r_sim_water_benchmark_speedup.png"><img src="r_sim_water_benchmark_speedup.png"
+  width="600" height="450"
+  alt="speedup for r.sim.water with different numbers of cells" border="0"></a>
+  <a href="r_sim_water_benchmark_efficiency.png"><img src="r_sim_water_benchmark_efficiency.png"
+  width="600" height="450"
+  alt="efficiency for r.sim.water with different numbers of cells" border="0"></a>
+  <br>
+  <i>Figure: Benchmark shows execution time, parallel speedup and
+    efficiency for different numbers of cells (33k, 131k and 525k);
+    shading shows the range of the 3 runs.</i>
+</div>
+
 <h2>EXAMPLE</h2>
 
 This example uses the

raster/r.sim/r.sim.water/r.sim.water.md

@@ -323,6 +323,28 @@ results differ slightly between thread counts and between repeated runs,
 since the order in which the walkers reach a cell depends on the
 threads. Use **nprocs=1** when results must be reproducible.
 
+### Performance
+
+To enable parallel processing, the user can specify the number of
+threads to be used with the **nprocs** parameter (default 1). Figures
+below show benchmark results for the `elevation` raster map of the
+SECREF North Carolina dataset at 4 m, 2 m and 1 m resolution with the
+default number of walkers, as the mean of 3 runs on Intel® Xeon® W-2295
+CPU @ 3.00GHz × 18. See the benchmark script in the source code for
+more details.
+
+The time step is derived from the cell size and the mean flow velocity,
+so the number of iterations, and with it the run time, depends on the
+terrain as well as on the number of cells. As a result, the benchmark
+results may vary depending on the study area.
+
+![time for r.sim.water with different numbers of cells](r_sim_water_benchmark_time.png)
+![speedup for r.sim.water with different numbers of cells](r_sim_water_benchmark_speedup.png)
+![efficiency for r.sim.water with different numbers of cells](r_sim_water_benchmark_efficiency.png)  
+*Figure: Benchmark shows execution time, parallel speedup and efficiency
+for different numbers of cells (33k, 131k and 525k); shading shows the
+range of the 3 runs.*
+
 ## EXAMPLE
 
 This example uses the