r.neighbors: implement parallelization with OpenMP (#1724)

OSGeo · Apr 26, 2022 · 1c7f667 · 1c7f667
1 parent 67e5157
commit 1c7f667
Show file tree

Hide file tree

Showing 16 changed files with 508 additions and 106 deletions.
diff --git a/python/grass/pygrass/modules/interface/module.py b/python/grass/pygrass/modules/interface/module.py
@@ -348,30 +348,30 @@ class Module(object):
     >>> neighbors.inputs.size = 5
     >>> neighbors.inputs.quantile = 0.5
     >>> neighbors.get_bash()
-    'r.neighbors input=mapA size=5 method=average weighting_function=none quantile=0.5 output=mapB'
+    'r.neighbors input=mapA size=5 method=average weighting_function=none quantile=0.5 nprocs=1 memory=300 output=mapB'
 
     >>> new_neighbors1 = copy.deepcopy(neighbors)
     >>> new_neighbors1.inputs.input = "mapD"
     >>> new_neighbors1.inputs.size = 3
     >>> new_neighbors1.inputs.quantile = 0.5
     >>> new_neighbors1.get_bash()
-    'r.neighbors input=mapD size=3 method=average weighting_function=none quantile=0.5 output=mapB'
+    'r.neighbors input=mapD size=3 method=average weighting_function=none quantile=0.5 nprocs=1 memory=300 output=mapB'
 
     >>> new_neighbors2 = copy.deepcopy(neighbors)
     >>> new_neighbors2(input="mapD", size=3, run_=False)
     Module('r.neighbors')
     >>> new_neighbors2.get_bash()
-    'r.neighbors input=mapD size=3 method=average weighting_function=none quantile=0.5 output=mapB'
+    'r.neighbors input=mapD size=3 method=average weighting_function=none quantile=0.5 nprocs=1 memory=300 output=mapB'
 
     >>> neighbors = Module("r.neighbors")
     >>> neighbors.get_bash()
-    'r.neighbors size=3 method=average weighting_function=none'
+    'r.neighbors size=3 method=average weighting_function=none nprocs=1 memory=300'
 
     >>> new_neighbors3 = copy.deepcopy(neighbors)
     >>> new_neighbors3(input="mapA", size=3, output="mapB", run_=False)
     Module('r.neighbors')
     >>> new_neighbors3.get_bash()
-    'r.neighbors input=mapA size=3 method=average weighting_function=none output=mapB'
+    'r.neighbors input=mapA size=3 method=average weighting_function=none nprocs=1 memory=300 output=mapB'
 
     >>> mapcalc = Module("r.mapcalc", expression="test_a = 1",
     ...                  overwrite=True, run_=False)

diff --git a/raster/r.neighbors/Makefile b/raster/r.neighbors/Makefile
@@ -2,8 +2,9 @@ MODULE_TOPDIR = ../..
 
 PGM = r.neighbors
 
-LIBES = $(STATSLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB)
+LIBES = $(STATSLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(OMPLIB)
 DEPENDENCIES = $(STATSDEP) $(RASTERDEP) $(GISDEP)
+EXTRA_CFLAGS = $(OMPCFLAGS)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

diff --git a/raster/r.neighbors/benchmark/benchmark_r_neighbors_memory.py b/raster/r.neighbors/benchmark/benchmark_r_neighbors_memory.py
@@ -0,0 +1,58 @@
+"""Benchmarking of r.neighbors
+
+@author Aaron Saw Min Sern
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+
+    reference = "r_neighbors_reference_map"
+    generate_map(rows=10000, cols=10000, fname=reference)
+    # Users can add more or modify existing reference maps
+    benchmark(0, "r.neighbors_0MB", results, reference)
+    benchmark(5, "r.neighbors_5MB", results, reference)
+    benchmark(10, "r.neighbors_10MB", results, reference)
+    benchmark(100, "r.neighbors_100MB", results, reference)
+    benchmark(300, "r.neighbors_300MB", results, reference)
+
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    bm.nprocs_plot(results, filename="r_neighbors_benchmark_memory.svg")
+
+
+def benchmark(memory, label, results, reference):
+    output = "benchmark_r_neighbors_nprocs"
+
+    module = Module(
+        "r.neighbors",
+        input=reference,
+        output=output,
+        size=9,
+        memory=memory,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=16, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=output)
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1)
+    # Generate using r.random.surface if r.surf.fractal fails
+    try:
+        print("Generating reference map using r.surf.fractal...")
+        Module("r.surf.fractal", output=fname)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/raster/r.neighbors/benchmark/benchmark_r_neighbors_nprocs.py b/raster/r.neighbors/benchmark/benchmark_r_neighbors_nprocs.py
@@ -0,0 +1,57 @@
+"""Benchmarking of r.neighbors
+
+@author Aaron Saw Min Sern
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+
+    # Users can add more or modify existing reference maps
+    benchmark(7071, "r.neighbors_50M", results)
+    benchmark(10000, "r.neighbors_100M", results)
+    benchmark(14142, "r.neighbors_200M", results)
+    benchmark(20000, "r.neighbors_400M", results)
+
+    bm.nprocs_plot(results, filename="r_neighbors_benchmark_nprocs.svg")
+
+
+def benchmark(size, label, results):
+    reference = "r_neighbors_reference_map"
+    output = "benchmark_r_neighbors_nprocs"
+
+    generate_map(rows=size, cols=size, fname=reference)
+    module = Module(
+        "r.neighbors",
+        input=reference,
+        output=output,
+        size=9,
+        memory=300,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=16, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    Module("g.remove", quiet=True, flags="f", type="raster", name=output)
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1)
+    # Generate using r.random.surface if r.surf.fractal fails
+    try:
+        print("Generating reference map using r.surf.fractal...")
+        Module("r.surf.fractal", output=fname)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/raster/r.neighbors/benchmark/benchmark_r_neighbors_size.py b/raster/r.neighbors/benchmark/benchmark_r_neighbors_size.py
@@ -0,0 +1,57 @@
+"""Benchmarking of r.neighbors
+
+@author Anna Petrasova
+"""
+
+from grass.exceptions import CalledModuleError
+from grass.pygrass.modules import Module
+from subprocess import DEVNULL
+
+import grass.benchmark as bm
+
+
+def main():
+    results = []
+
+    reference = "r_neighbors_reference_map"
+    generate_map(rows=10000, cols=10000, fname=reference)
+    # Users can add more or modify existing reference maps
+    benchmark(3, "r.neighbors_3x3", results, reference)
+    benchmark(9, "r.neighbors_9x9", results, reference)
+    benchmark(15, "r.neighbors_15x15", results, reference)
+    benchmark(31, "r.neighbors_31x31", results, reference)
+
+    Module("g.remove", quiet=True, flags="f", type="raster", name=reference)
+    bm.nprocs_plot(results, filename="r_neighbors_benchmark_size.svg")
+
+
+def benchmark(size, label, results, reference):
+    output = "benchmark_r_neighbors_nprocs"
+
+    module = Module(
+        "r.neighbors",
+        input=reference,
+        output=output,
+        size=size,
+        memory=300,
+        run_=False,
+        stdout_=DEVNULL,
+        overwrite=True,
+    )
+    results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=16, repeat=3))
+    Module("g.remove", quiet=True, flags="f", type="raster", name=output)
+
+
+def generate_map(rows, cols, fname):
+    Module("g.region", flags="p", s=0, n=rows, w=0, e=cols, res=1)
+    # Generate using r.random.surface if r.surf.fractal fails
+    try:
+        print("Generating reference map using r.surf.fractal...")
+        Module("r.surf.fractal", output=fname)
+    except CalledModuleError:
+        print("r.surf.fractal fails, using r.random.surface instead...")
+        Module("r.random.surface", output=fname)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/raster/r.neighbors/bufs.c b/raster/r.neighbors/bufs.c
@@ -12,31 +12,34 @@
 
 int allocate_bufs(void)
 {
-    int i;
+    int i, t;
     int bufsize;
 
     bufsize = (Rast_window_cols() + 2 * ncb.dist) * sizeof(DCELL);
 
-    ncb.buf = (DCELL **) G_malloc(ncb.nsize * sizeof(DCELL *));
-    for (i = 0; i < ncb.nsize; i++) {
-	ncb.buf[i] = (DCELL *) G_malloc(bufsize);
-	Rast_set_d_null_value(ncb.buf[i], Rast_window_cols() + 2 * ncb.dist);
+    ncb.buf = G_malloc(ncb.threads * sizeof(DCELL **));
+    for (t = 0; t < ncb.threads; t++) {
+        ncb.buf[t] = (DCELL **) G_malloc(ncb.nsize * sizeof(DCELL *));
+        for (i = 0; i < ncb.nsize; i++) {
+            ncb.buf[t][i] = (DCELL *) G_malloc(bufsize);
+            Rast_set_d_null_value(ncb.buf[t][i], Rast_window_cols() + 2 * ncb.dist);
+        }
     }
 
     return 0;
 }
 
-int rotate_bufs(void)
+int rotate_bufs(int thread_id)
 {
     DCELL *temp;
     int i;
 
-    temp = ncb.buf[0];
+    temp = ncb.buf[thread_id][0];
 
     for (i = 1; i < ncb.nsize; i++)
-	ncb.buf[i - 1] = ncb.buf[i];
+	ncb.buf[thread_id][i - 1] = ncb.buf[thread_id][i];
 
-    ncb.buf[ncb.nsize - 1] = temp;
+    ncb.buf[thread_id][ncb.nsize - 1] = temp;
 
     return 0;
 }
diff --git a/raster/r.neighbors/gather.c b/raster/r.neighbors/gather.c
@@ -45,7 +45,7 @@ void weights_mask(void)
 	    ncb.mask[i][j] = ncb.weights[i][j] != 0;
 }
 
-int gather(DCELL *values, int offset)
+int gather(DCELL *values, int offset, int thread_id)
 {
     int row, col;
     int n = 0;
@@ -58,7 +58,7 @@ int gather(DCELL *values, int offset)
 	    if (ncb.mask && !ncb.mask[row][col])
 		continue;
 
-	    values[n] = ncb.buf[row][offset + col];
+	    values[n] = ncb.buf[thread_id][row][offset + col];
 
 	    n++;
 	}
@@ -67,7 +67,7 @@ int gather(DCELL *values, int offset)
     return n;
 }
 
-int gather_w(DCELL *values, DCELL (*values_w)[2], int offset)
+int gather_w(DCELL *values, DCELL (*values_w)[2], int offset, int thread_id)
 {
     int row, col;
     int n = 0;
@@ -77,7 +77,7 @@ int gather_w(DCELL *values, DCELL (*values_w)[2], int offset)
 
     for (row = 0; row < ncb.nsize; row++) {
 	for (col = 0; col < ncb.nsize; col++) {
-	    values[n] = values_w[n][0] = ncb.buf[row][offset + col];
+	    values[n] = values_w[n][0] = ncb.buf[thread_id][row][offset + col];
 	    values_w[n][1] = ncb.weights[row][col];
 
 	    n++;

diff --git a/raster/r.neighbors/local_proto.h b/raster/r.neighbors/local_proto.h
@@ -1,15 +1,17 @@
+#include <grass/gis.h>
+
 /* bufs.c */
 extern int allocate_bufs(void);
-extern int rotate_bufs(void);
+extern int rotate_bufs(int);
 
 /* gather */
 extern void circle_mask(void);
 extern void weights_mask(void);
-extern int gather(DCELL *, int);
-extern int gather_w(DCELL *, DCELL(*)[2], int);
+extern int gather(DCELL *, int, int);
+extern int gather_w(DCELL *, DCELL(*)[2], int, int);
 
 /* readcell.c */
-extern int readcell(int, int, int, int);
+extern int readcell(int, int, int, int, int);
 
 /* divr_cats.c */
 extern int divr_cats(void);