diff --git a/raster/r.univar/Makefile b/raster/r.univar/Makefile index fd1e7fd3317..e6d2e5ba831 100644 --- a/raster/r.univar/Makefile +++ b/raster/r.univar/Makefile @@ -1,9 +1,10 @@ MODULE_TOPDIR = ../.. -LIBES2 = $(RASTERLIB) $(GISLIB) $(MATHLIB) -LIBES3 = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) +LIBES2 = $(RASTERLIB) $(GISLIB) $(MATHLIB) $(OMPLIB) +LIBES3 = $(RASTER3DLIB) $(RASTERLIB) $(GISLIB) $(MATHLIB) $(OMPLIB) DEPENDENCIES = $(RASTER3DDEP) $(GISDEP) $(RASTERDEP) +EXTRA_CFLAGS = $(OMPCFLAGS) PROGRAMS = r.univar r3.univar diff --git a/raster/r.univar/benchmark/benchmark_r_univar.py b/raster/r.univar/benchmark/benchmark_r_univar.py new file mode 100644 index 00000000000..2765b067783 --- /dev/null +++ b/raster/r.univar/benchmark/benchmark_r_univar.py @@ -0,0 +1,54 @@ +"""Benchmarking of r.univar + +@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.univar_50M", results) + benchmark(10000, "r.univar_100M", results) + benchmark(14142, "r.univar_200M", results) + benchmark(20000, "r.univar_400M", results) + + bm.nprocs_plot(results) + + +def benchmark(size, label, results): + reference = "r_univar_reference_map" + output = "benchmark_r_univar_nprocs" + + generate_map(rows=size, cols=size, fname=reference) + module = Module( + "r.univar", + map=reference, + 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.univar/globals.h b/raster/r.univar/globals.h index d234a7fc62d..d152ae05be8 100644 --- a/raster/r.univar/globals.h +++ b/raster/r.univar/globals.h @@ -56,7 +56,7 @@ typedef struct typedef struct { struct Option *inputfile, *zonefile, *percentile, *output_file, - *separator; + *separator, *nprocs; struct Flag *shell_style, *extended, *table, *use_rast_region; } param_type; diff --git a/raster/r.univar/r.univar.html b/raster/r.univar/r.univar.html index 399fd495397..33e34169086 100644 --- a/raster/r.univar/r.univar.html +++ b/raster/r.univar/r.univar.html @@ -48,6 +48,28 @@

NOTES

map and uploads statistics to new attribute columns, see v.rast.stats. +

PERFORMANCE

+ +

+r.univar suports parallel processing using OpenMP. The user +can specify the number of threads to be used with the nprocs parameter. +However, parallelization is disabled when the -e extended statistics +flag is used. + +

+Due to the differences in summation order, users may encounter small floating points +discrepancies when r.univar is run on very large raster files when different +nprocs parameters are used. However, since the work allocation among threads +is static, users should expect to have the same results when run with the same +number of threads. + +

+ benchmark for number of cells +
+ Figure: Benchmark shows execution time for different + number of cells and cores. See benchmark scripts in source code. +
+

EXAMPLES

Univariate statistics

diff --git a/raster/r.univar/r.univar_main.c b/raster/r.univar/r.univar_main.c index 1e7740c3734..436e2c80c03 100644 --- a/raster/r.univar/r.univar_main.c +++ b/raster/r.univar/r.univar_main.c @@ -15,12 +15,18 @@ * This program is a replacement for the r.univar shell script */ +#if defined(_OPENMP) +#include +#endif + #include #include +#include #include "globals.h" param_type param; zone_type zone_info; +int nprocs; /* ************************************************************************* */ /* Set up the arguments we are expecting ********************************** */ @@ -52,6 +58,8 @@ void set_params() _("Percentile to calculate (requires extended statistics flag)"); param.percentile->guisection = _("Extended"); + param.nprocs = G_define_standard_option(G_OPT_M_NPROCS); + param.separator = G_define_standard_option(G_OPT_F_SEP); param.separator->guisection = _("Formatting"); @@ -82,8 +90,8 @@ void set_params() static int open_raster(const char *infile); static univar_stat *univar_stat_with_percentiles(int map_type); -static void process_raster(univar_stat * stats, int fd, int fdz, - const struct Cell_head *region); +static void process_raster(univar_stat * stats, int *fd, + int *fdz, const struct Cell_head *region); /* *************************************************************** */ /* **** the main functions for r.univar ************************** */ @@ -96,9 +104,10 @@ int main(int argc, char *argv[]) struct GModule *module; univar_stat *stats; char **p, *z; - int fd, fdz, cell_type, min, max; + int *fd, *fdz, cell_type, min, max; struct Range zone_range; const char *mapset, *name; + int t; G_gisinit(argv[0]); @@ -107,6 +116,7 @@ int main(int argc, char *argv[]) G_add_keyword(_("statistics")); G_add_keyword(_("univariate statistics")); G_add_keyword(_("zonal statistics")); + G_add_keyword(_("parallel")); module->label = _("Calculates univariate statistics from the non-null cells of a raster map."); @@ -132,6 +142,26 @@ int main(int argc, char *argv[]) } } + /* set nprocs parameter */ + sscanf(param.nprocs->answer, "%d", &nprocs); + if (nprocs < 1) + G_fatal_error(_("<%d> is not valid number of nprocs."), nprocs); +#if defined(_OPENMP) + if (param.extended->answer) { + /* Calculation of extended statistics is not parallelized yet */ + if (nprocs > 1) + G_warning(_("Computing extending statistics is not parallelized yet. Ignoring threads setting.")); + nprocs = 1; + } + else { + omp_set_num_threads(nprocs); + } +#else + if (nprocs != 1) + G_warning(_("GRASS is compiled without OpenMP support. Ignoring threads setting.")); + nprocs = 1; +#endif + /* table field separator */ zone_info.sep = G_option_to_separator(param.separator); @@ -139,15 +169,18 @@ int main(int argc, char *argv[]) zone_info.max = 0.0 / 0.0; /* set to nan as default */ zone_info.n_zones = 0; - fdz = -1; + fdz = NULL; + fd = G_malloc(nprocs * sizeof(int)); /* open zoning raster */ if ((z = param.zonefile->answer)) { mapset = G_find_raster2(z, ""); - fdz = open_raster(z); + fdz = G_malloc(nprocs * sizeof(int)); + for (t = 0; t < nprocs; t++) + fdz[t] = open_raster(z); - cell_type = Rast_get_map_type(fdz); + cell_type = Rast_get_map_type(fdz[0]); if (cell_type != CELL_TYPE) G_fatal_error("Zoning raster must be of type CELL"); @@ -187,11 +220,12 @@ int main(int argc, char *argv[]) G_get_window(®ion); } - fd = open_raster(*p); + for (t = 0; t < nprocs; t++) + fd[t] = open_raster(*p); if (map_type != -1) { /* NB: map_type must match when doing extended stats */ - int this_type = Rast_get_map_type(fd); + int this_type = Rast_get_map_type(fd[0]); assert(this_type > -1); if (map_type < -1) { @@ -208,12 +242,15 @@ int main(int argc, char *argv[]) process_raster(stats, fd, fdz, ®ion); /* close input raster */ - Rast_close(fd); + for (t = 0; t < nprocs; t++) + Rast_close(fd[t]); } /* close zoning raster */ - if (z) - Rast_close(fdz); + if (z) { + for (t = 0; t < nprocs; t++) + Rast_close(fdz[t]); + } /* create the output */ if (param.table->answer) @@ -238,6 +275,7 @@ static int open_raster(const char *infile) } fd = Rast_open_old(infile, mapset); + G_free((void *)mapset); return fd; } @@ -266,132 +304,214 @@ static univar_stat *univar_stat_with_percentiles(int map_type) } static void -process_raster(univar_stat * stats, int fd, int fdz, +process_raster(univar_stat * stats, int *fd, int *fdz, const struct Cell_head *region) { /* use G_window_rows(), G_window_cols() here? */ const unsigned int rows = region->rows; const unsigned int cols = region->cols; - const RASTER_MAP_TYPE map_type = Rast_get_map_type(fd); + const RASTER_MAP_TYPE map_type = Rast_get_map_type(fd[0]); const size_t value_sz = Rast_cell_size(map_type); + int t; unsigned int row; - void *raster_row; - CELL *zoneraster_row = NULL; + void **raster_row; + CELL **zoneraster_row = NULL; int n_zones = zone_info.n_zones; + int n_alloc = n_zones ? n_zones : 1; - raster_row = Rast_allocate_buf(map_type); - if (n_zones) - zoneraster_row = Rast_allocate_c_buf(); - - for (row = 0; row < rows; row++) { - void *ptr; - CELL *zptr = NULL; - unsigned int col; - - Rast_get_row(fd, raster_row, row, map_type); + raster_row = G_malloc(nprocs * sizeof(void *)); + if (n_zones) { + zoneraster_row = G_malloc(nprocs * sizeof(CELL *)); + } + for (t = 0; t < nprocs; t++) { + raster_row[t] = Rast_allocate_buf(map_type); if (n_zones) { - Rast_get_c_row(fdz, zoneraster_row, row); - zptr = zoneraster_row; + zoneraster_row[t] = Rast_allocate_c_buf(); } + } - ptr = raster_row; +#if defined(_OPENMP) + omp_lock_t *minmax = G_malloc(n_alloc * sizeof(omp_lock_t)); - for (col = 0; col < cols; col++) { - double val; - int zone = 0; + for (int i = 0; i < n_alloc; i++) { + omp_init_lock(&(minmax[i])); + } +#endif + int computed = 0; + +#pragma omp parallel if(nprocs > 1) + { + int i; + int t_id = 0; + +#if defined(_OPENMP) + if (!param.extended->answer) { + t_id = omp_get_thread_num(); + } +#endif + size_t *n = G_calloc(n_alloc, sizeof(size_t)); + double *sum = G_calloc(n_alloc, sizeof(double)); + double *sumsq = G_calloc(n_alloc, sizeof(double)); + double *sum_abs = G_calloc(n_alloc, sizeof(double)); + size_t *size = G_calloc(n_alloc, sizeof(size_t)); + double *min = G_malloc(n_alloc * sizeof(double)); + double *max = G_malloc(n_alloc * sizeof(double)); + + for (i = 0; i < n_alloc; i++) { + max[i] = DBL_MIN; + min[i] = DBL_MAX; + } + +#pragma omp for schedule(static) + for (row = 0; row < rows; row++) { + void *ptr; + CELL *zptr = NULL; + unsigned int col; + + Rast_get_row(fd[t_id], raster_row[t_id], row, map_type); if (n_zones) { - /* skip NULL cells in zone map */ - if (Rast_is_c_null_value(zptr)) { + Rast_get_c_row(fdz[t_id], zoneraster_row[t_id], row); + zptr = zoneraster_row[t_id]; + } + + ptr = raster_row[t_id]; + + for (col = 0; col < cols; col++) { + double val; + int zone = 0; + + if (n_zones) { + /* skip NULL cells in zone map */ + if (Rast_is_c_null_value(zptr)) { + ptr = G_incr_void_ptr(ptr, value_sz); + zptr++; + continue; + } + zone = *zptr - zone_info.min; + } + + /* count all including NULL cells in input map */ + size[zone]++; + + /* can't do stats with NULL cells in input map */ + if (Rast_is_null_value(ptr, map_type)) { ptr = G_incr_void_ptr(ptr, value_sz); - zptr++; + if (n_zones) + zptr++; continue; } - zone = *zptr - zone_info.min; - } - /* count all including NULL cells in input map */ - stats[zone].size++; + if (param.extended->answer) { + /* check allocated memory */ + /* parallelization is disabled, local variable reflects global state */ + if (stats[zone].n + n[zone] >= stats[zone].n_alloc) { + stats[zone].n_alloc += 1000; + size_t msize; + + switch (map_type) { + case DCELL_TYPE: + msize = stats[zone].n_alloc * sizeof(DCELL); + stats[zone].dcell_array = + (DCELL *) G_realloc((void *)stats[zone]. + dcell_array, msize); + stats[zone].nextp = + (void *)&(stats[zone].dcell_array[stats[zone].n + n[zone]]); + break; + case FCELL_TYPE: + msize = stats[zone].n_alloc * sizeof(FCELL); + stats[zone].fcell_array = + (FCELL *) G_realloc((void *)stats[zone]. + fcell_array, msize); + stats[zone].nextp = + (void *)&(stats[zone].fcell_array[stats[zone].n + n[zone]]); + break; + case CELL_TYPE: + msize = stats[zone].n_alloc * sizeof(CELL); + stats[zone].cell_array = + (CELL *) G_realloc((void *)stats[zone]. + cell_array, msize); + stats[zone].nextp = + (void *)&(stats[zone].cell_array[stats[zone].n + n[zone]]); + break; + default: + break; + } + } + /* put the value into stats->XXXcell_array */ + memcpy(stats[zone].nextp, ptr, value_sz); + stats[zone].nextp = + G_incr_void_ptr(stats[zone].nextp, value_sz); + } + + val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr) + : (map_type == FCELL_TYPE) ? *((FCELL *) ptr) + : *((CELL *) ptr)); + + sum[zone] += val; + sumsq[zone] += val * val; + sum_abs[zone] += fabs(val); + + if (val > max[zone]) + max[zone] = val; + if (val < min[zone]) + min[zone] = val; - /* can't do stats with NULL cells in input map */ - if (Rast_is_null_value(ptr, map_type)) { ptr = G_incr_void_ptr(ptr, value_sz); if (n_zones) zptr++; - continue; - } + n[zone]++; - if (param.extended->answer) { - /* check allocated memory */ - if (stats[zone].n >= stats[zone].n_alloc) { - stats[zone].n_alloc += 1000; - size_t msize; - - switch (map_type) { - case DCELL_TYPE: - msize = stats[zone].n_alloc * sizeof(DCELL); - stats[zone].dcell_array = - (DCELL *) G_realloc((void *)stats[zone]. - dcell_array, msize); - stats[zone].nextp = - (void *)&(stats[zone].dcell_array[stats[zone].n]); - break; - case FCELL_TYPE: - msize = stats[zone].n_alloc * sizeof(FCELL); - stats[zone].fcell_array = - (FCELL *) G_realloc((void *)stats[zone]. - fcell_array, msize); - stats[zone].nextp = - (void *)&(stats[zone].fcell_array[stats[zone].n]); - break; - case CELL_TYPE: - msize = stats[zone].n_alloc * sizeof(CELL); - stats[zone].cell_array = - (CELL *) G_realloc((void *)stats[zone].cell_array, - msize); - stats[zone].nextp = - (void *)&(stats[zone].cell_array[stats[zone].n]); - break; - default: - break; - } - } - /* put the value into stats->XXXcell_array */ - memcpy(stats[zone].nextp, ptr, value_sz); - stats[zone].nextp = - G_incr_void_ptr(stats[zone].nextp, value_sz); + } /* end column loop */ + if (!(param.shell_style->answer)) { +#pragma omp atomic update + computed++; + G_percent(computed, rows, 2); } - - val = ((map_type == DCELL_TYPE) ? *((DCELL *) ptr) - : (map_type == FCELL_TYPE) ? *((FCELL *) ptr) - : *((CELL *) ptr)); - - stats[zone].sum += val; - stats[zone].sumsq += val * val; - stats[zone].sum_abs += fabs(val); - - if (stats[zone].first) { - stats[zone].max = val; - stats[zone].min = val; - stats[zone].first = FALSE; + } /* end row loop */ + + for (i = 0; i < n_alloc; i++) { +#pragma omp atomic update + stats[i].n += n[i]; +#pragma omp atomic update + stats[i].size += size[i]; +#pragma omp atomic update + stats[i].sum += sum[i]; +#pragma omp atomic update + stats[i].sumsq += sumsq[i]; +#pragma omp atomic update + stats[i].sum_abs += sum_abs[i]; + +#if defined(_OPENMP) + omp_set_lock(&(minmax[i])); +#endif + if (stats[i].max < max[i] || + (stats[i].max != stats[i].max && max[i] != DBL_MIN)) { + stats[i].max = max[i]; } - else { - if (val > stats[zone].max) - stats[zone].max = val; - if (val < stats[zone].min) - stats[zone].min = val; + if (stats[i].min > min[i] || + (stats[i].min != stats[i].min && min[i] != DBL_MAX)) { + stats[i].min = min[i]; } +#if defined(_OPENMP) + omp_unset_lock(&(minmax[i])); +#endif + } - ptr = G_incr_void_ptr(ptr, value_sz); - if (n_zones) - zptr++; - stats[zone].n++; + } /* end parallel region */ + + for (t = 0; t < nprocs; t++) { + G_free(raster_row[t]); + } + G_free(raster_row); + if (n_zones) { + for (t = 0; t < nprocs; t++) { + G_free(zoneraster_row[t]); } - if (!(param.shell_style->answer)) - G_percent(row, rows, 2); + G_free(zoneraster_row); } if (!(param.shell_style->answer)) - G_percent(rows, rows, 2); /* finish it off */ + G_percent(rows, rows, 2); } diff --git a/raster/r.univar/r_univar_benchmark_size.png b/raster/r.univar/r_univar_benchmark_size.png new file mode 100644 index 00000000000..a9ca8b377af Binary files /dev/null and b/raster/r.univar/r_univar_benchmark_size.png differ diff --git a/raster/r.univar/testsuite/test_r_univar.py b/raster/r.univar/testsuite/test_r_univar.py index 3241c10dd87..34437261408 100644 --- a/raster/r.univar/testsuite/test_r_univar.py +++ b/raster/r.univar/testsuite/test_r_univar.py @@ -47,7 +47,16 @@ def test_1(self): sum=1547100""" self.assertRasterFitsUnivar( - raster="map_a", reference=univar_string, precision=3 + raster="map_a", reference=univar_string, precision=6 + ) + self.assertModuleKeyValue( + module="r.univar", + map="map_a", + flags="g", + nprocs=4, + reference=univar_string, + precision=6, + sep="=", ) def test_2(self): @@ -64,7 +73,16 @@ def test_2(self): self.runModule("g.region", res=10) self.assertRasterFitsUnivar( - raster="map_a", reference=univar_string, precision=3 + raster="map_a", reference=univar_string, precision=6 + ) + self.assertModuleKeyValue( + module="r.univar", + map="map_a", + flags="g", + nprocs=4, + reference=univar_string, + precision=6, + sep="=", ) def test_3(self): @@ -89,7 +107,16 @@ def test_3(self): map="map_a", flags="rg", reference=univar_string, - precision=3, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map="map_a", + flags="rg", + nprocs=4, + reference=univar_string, + precision=6, sep="=", ) @@ -110,7 +137,16 @@ def test_multiple_1(self): map=["map_a", "map_b"], flags="rg", reference=univar_string, - precision=3, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + flags="rg", + nprocs=4, + reference=univar_string, + precision=6, sep="=", ) @@ -132,7 +168,16 @@ def test_multiple_2(self): map=["map_a", "map_b"], flags="g", reference=univar_string, - precision=3, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + flags="g", + nprocs=4, + reference=univar_string, + precision=6, sep="=", ) @@ -159,13 +204,22 @@ def test_multiple_3(self): map=["map_a", "map_b"], flags="rg", reference=univar_string, - precision=3, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + flags="rg", + nprocs=4, + reference=univar_string, + precision=6, sep="=", ) def test_1_zone(self): """ - multiple maps and zone + one map and zone :return: """ @@ -198,7 +252,133 @@ def test_1_zone(self): zones="zone_map", flags="g", reference=univar_string, - precision=3, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a"], + zones="zone_map", + flags="g", + nprocs=4, + reference=univar_string, + precision=6, + sep="=", + ) + + def test_2_zone(self): + """ + multiple maps and zone + :return: + """ + + # Output of r.univar + univar_string = """zone=1; + n=3420 + null_cells=0 + cells=3420 + min=102 + max=309 + range=207 + mean=205.5 + mean_of_abs=205.5 + stddev=56.6119834192962 + variance=3204.91666666667 + coeff_var=27.5484104230152 + sum=702810 + zone=2; + n=12780 + null_cells=0 + cells=3420 + min=121 + max=380 + range=259 + mean=250.5 + mean_of_abs=250.5 + stddev=59.9576239244574 + variance=3594.91666666667 + coeff_var=23.9351792113602 + sum=3201390""" + + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + zones="zone_map", + flags="g", + reference=univar_string, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + zones="zone_map", + flags="g", + nprocs=4, + reference=univar_string, + precision=6, + sep="=", + ) + + def test_3_zone(self): + """ + multiple maps and zone + :return: + """ + + # Output of r.univar + univar_string = """zone=1; + n=3420 + null_cells=0 + cells=3420 + min=102 + max=309 + range=207 + mean=205.5 + mean_of_abs=205.5 + stddev=56.6119834192962 + variance=3204.91666666667 + coeff_var=27.5484104230152 + sum=702810 + first_quartile=155 + median=205.5 + third_quartile=255 + percentile_90=282 + zone=2; + n=12780 + null_cells=0 + cells=3420 + min=121 + max=380 + range=259 + mean=250.5 + mean_of_abs=250.5 + stddev=59.9576239244574 + variance=3594.91666666667 + coeff_var=23.9351792113602 + sum=3201390 + first_quartile=200 + median=250.5 + third_quartile=300 + percentile_90=330""" + + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + zones="zone_map", + flags="ge", + reference=univar_string, + precision=6, + sep="=", + ) + self.assertModuleKeyValue( + module="r.univar", + map=["map_a", "map_b"], + zones="zone_map", + flags="ge", + nprocs=4, + reference=univar_string, + precision=6, sep="=", )