From 6af4ca5f527b0671b5e32133122ad9c76f1ee377 Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Thu, 20 Jun 2024 21:28:56 -0400 Subject: [PATCH 01/25] expand buffer and cache to multiple rows for multithreading --- raster/r.mapcalc/Makefile | 6 +- raster/r.mapcalc/evaluate.c | 101 +++++++++++++++++++++++++--------- raster/r.mapcalc/expression.h | 10 ++-- raster/r.mapcalc/globals.h | 3 +- raster/r.mapcalc/main.c | 16 +++++- raster/r.mapcalc/map.c | 32 +++++++++-- raster/r.mapcalc/map3.c | 10 +++- raster/r.mapcalc/mapcalc.h | 3 + raster/r.mapcalc/xarea.c | 13 ++++- raster/r.mapcalc/xcoor.c | 11 +++- raster/r.mapcalc/xcoor3.c | 11 +++- raster/r.mapcalc/xrowcol.c | 11 +++- 12 files changed, 182 insertions(+), 45 deletions(-) diff --git a/raster/r.mapcalc/Makefile b/raster/r.mapcalc/Makefile index 7b8bbf90d16..d8ed91eeb1c 100644 --- a/raster/r.mapcalc/Makefile +++ b/raster/r.mapcalc/Makefile @@ -12,8 +12,10 @@ r3_mapcalc_OBJS := $(filter-out map.o xcoor.o xres.o, $(AUTO_OBJS)) include $(MODULE_TOPDIR)/include/Make/Multi.make EXTRA_CFLAGS = $(READLINEINCPATH) $(PTHREADINCPATH) -LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) -LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) +LIBES2 = $(CALCLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB) +LIBES3 = $(CALCLIB) $(RASTER3DLIB) $(GISLIB) $(RASTERLIB) $(BTREELIB) $(READLINELIBPATH) $(READLINELIB) $(HISTORYLIB) $(PTHREADLIBPATH) $(PTHREADLIB) $(OPENMP_LIBPATH) $(OPENMP_LIB) +EXTRA_CFLAGS = $(OPENMP_CFLAGS) +EXTRA_INC = $(OPENMP_INCPATH) default: multi diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c index d8eadc73731..20dc2f5d325 100644 --- a/raster/r.mapcalc/evaluate.c +++ b/raster/r.mapcalc/evaluate.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include @@ -12,7 +16,8 @@ /****************************************************************************/ -int current_depth, current_row; +int current_depth; +int *current_row; int depths, rows; /* Local variables for map management */ @@ -69,10 +74,18 @@ void extract_maps(expression *e) static void allocate_buf(expression *e) { - e->buf = G_malloc(columns * Rast_cell_size(e->res_type)); + + int threads = 1; +#if defined(_OPENMP) + threads = omp_get_max_threads(); +#endif + + e->buf = (void **)G_malloc(sizeof(void *) * threads); + for (int t = 0; t < threads; t++) + e->buf[t] = G_malloc(columns * Rast_cell_size(e->res_type)); } -static void set_buf(expression *e, void *buf) +static void set_buf(expression *e, void **buf) { e->buf = buf; } @@ -102,8 +115,7 @@ static void initialize_function(expression *e) int i; allocate_buf(e); - - e->data.func.argv = G_malloc((e->data.func.argc + 1) * sizeof(void *)); + e->data.func.argv = G_malloc((e->data.func.argc + 1) * sizeof(void **)); e->data.func.argv[0] = e->buf; for (i = 1; i <= e->data.func.argc; i++) { @@ -163,9 +175,14 @@ static void end_evaluate(struct expression *e) static void evaluate_constant(expression *e) { - int *ibuf = e->buf; - float *fbuf = e->buf; - double *dbuf = e->buf; + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + + int *ibuf = e->buf[tid]; + float *fbuf = e->buf[tid]; + double *dbuf = e->buf[tid]; int i; switch (e->res_type) { @@ -195,15 +212,25 @@ static void evaluate_variable(expression *e UNUSED) static void evaluate_map(expression *e) { - get_map_row( - e->data.map.idx, e->data.map.mod, current_depth + e->data.map.depth, - current_row + e->data.map.row, e->data.map.col, e->buf, e->res_type); + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + + get_map_row(e->data.map.idx, e->data.map.mod, + current_depth + e->data.map.depth, + current_row[tid] + e->data.map.row, e->data.map.col, + e->buf[tid], e->res_type); } static void evaluate_function(expression *e) { int i; int res; + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif if (e->data.func.argc > 1 && e->data.func.func != f_eval) { for (i = 1; i <= e->data.func.argc; i++) @@ -216,8 +243,19 @@ static void evaluate_function(expression *e) for (i = 1; i <= e->data.func.argc; i++) evaluate(e->data.func.args[i]); - res = (*e->data.func.func)(e->data.func.argc, e->data.func.argt, - e->data.func.argv); + /* copy the argv in the individual thread */ + void **thread_argv = G_malloc((e->data.func.argc + 1) * sizeof(void *)); + for (i = 0; i < e->data.func.argc + 1; i++) + thread_argv[i] = e->data.func.argv[i][tid]; + + res = + (*e->data.func.func)(e->data.func.argc, e->data.func.argt, thread_argv); + + /* copy the results from thread_argv to e */ + for (i = 0; i < e->data.func.argc + 1; i++) + e->data.func.argv[i][tid] = thread_argv[i]; + + G_free(thread_argv); switch (res) { case E_ARG_LO: @@ -303,6 +341,12 @@ void execute(expr_list *ee) int verbose = isatty(2); expr_list *l; int count, n; + int threads = 1; + +#if defined(_OPENMP) + threads = omp_get_max_threads(); +#endif + current_row = (int *)G_malloc(sizeof(int) * threads); exprs = ee; G_add_error_handler(error_handler, NULL); @@ -361,27 +405,34 @@ void execute(expr_list *ee) count = rows * depths; n = 0; - G_init_workers(); - for (current_depth = 0; current_depth < depths; current_depth++) { - for (current_row = 0; current_row < rows; current_row++) { +#pragma omp parallel for default(shared) schedule(static, 1) ordered + for (int row = 0; row < rows; row++) { if (verbose) G_percent(n, count, 2); - for (l = ee; l; l = l->next) { + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + current_row[tid] = row; + for (l = ee; l != NULL; l = l->next) { expression *e = l->exp; int fd; evaluate(e); - - if (e->type != expr_type_binding) - continue; - - fd = e->data.bind.fd; - put_map_row(fd, e->buf, e->res_type); +#pragma omp ordered + { + if (e->type == expr_type_binding) { + fd = e->data.bind.fd; + put_map_row(fd, e->buf[tid], e->res_type); + } + } + } +#pragma omp critical + { + n++; } - - n++; } } diff --git a/raster/r.mapcalc/expression.h b/raster/r.mapcalc/expression.h index 9f049cc7a6a..e8ec3bea5d5 100644 --- a/raster/r.mapcalc/expression.h +++ b/raster/r.mapcalc/expression.h @@ -35,10 +35,10 @@ typedef struct expr_data_func { const char *oper; int prec; func_t *func; - int argc; - struct expression **args; - int *argt; - void **argv; + int argc; /* number of args in the whole expression */ + struct expression **args; /* array of expressions */ + int *argt; /* type of expressions */ + void ***argv; /* values in e->buf for each expression */ } expr_data_func; typedef struct expr_data_bind { @@ -50,7 +50,7 @@ typedef struct expr_data_bind { typedef struct expression { int type; int res_type; - void *buf; + void **buf; union { expr_data_const con; expr_data_var var; diff --git a/raster/r.mapcalc/globals.h b/raster/r.mapcalc/globals.h index 0d49c166198..351ebbabe36 100644 --- a/raster/r.mapcalc/globals.h +++ b/raster/r.mapcalc/globals.h @@ -6,7 +6,8 @@ extern long seed_value; extern long seeded; extern int region_approach; -extern int current_depth, current_row; +extern int current_depth; +extern int *current_row; extern int depths, rows, columns; #endif /* __GLOBALS_H_ */ diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index ea57d561d10..25926e865c3 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -11,6 +11,9 @@ * for details. * *****************************************************************************/ +#if defined(_OPENMP) +#include +#endif #include #include @@ -59,10 +62,11 @@ static expr_list *parse_file(const char *filename) int main(int argc, char **argv) { struct GModule *module; - struct Option *expr, *file, *seed, *region; + struct Option *expr, *file, *seed, *region, *nprocs; struct Flag *random, *describe; int all_ok; char *desc; + int threads = 1; G_gisinit(argv[0]); @@ -119,6 +123,8 @@ int main(int argc, char **argv) describe->key = 'l'; describe->description = _("List input and output maps"); + nprocs = G_define_standard_option(G_OPT_M_NPROCS); + if (argc == 1) { char **p = G_malloc(3 * sizeof(char *)); @@ -183,6 +189,14 @@ int main(int argc, char **argv) } pre_exec(); + threads = atoi(nprocs->answer); +#if defined(_OPENMP) + omp_set_num_threads(threads); + G_message(_("Computing in parallel, number of threads: %d"), + omp_get_max_threads()); +#else + G_message(_("Number of threads: 1")); +#endif execute(result); post_exec(); diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 118453c0610..884ab2b5520 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include @@ -60,7 +64,7 @@ struct map { struct Categories cats; struct Colors colors; BTREE btree; - struct row_cache cache; + struct row_cache *caches; #ifdef HAVE_PTHREAD_H pthread_mutex_t mutex; #endif @@ -383,13 +387,19 @@ static void translate_from_cats(struct map *m, CELL *cell, DCELL *xcell, static void setup_map(struct map *m) { int nrows = m->max_row - m->min_row + 1; - + int threads = 1; #ifdef HAVE_PTHREAD_H pthread_mutex_init(&m->mutex, NULL); #endif +#ifdef _OPENMP + threads = omp_get_max_threads(); +#endif + m->caches = + (struct row_cache *)G_malloc(threads * sizeof(struct row_cache)); if (nrows > 1 && nrows <= max_rows_in_memory) { - cache_setup(&m->cache, m->fd, nrows); + for (int i = 0; i < threads; i++) + cache_setup(&m->caches[i], m->fd, nrows); m->use_rowio = 1; } else @@ -426,8 +436,13 @@ static void read_map(struct map *m, void *buf, int res_type, int row, int col) return; } + int tid = 0; +#ifdef _OPENMP + tid = omp_get_thread_num(); +#endif + if (m->use_rowio) - cache_get(&m->cache, buf, row, res_type); + cache_get(&m->caches[tid], buf, row, res_type); else read_row(m->fd, buf, row, res_type); @@ -437,6 +452,11 @@ static void read_map(struct map *m, void *buf, int res_type, int row, int col) static void close_map(struct map *m) { + int threads = 1; +#ifdef _OPENMP + threads = omp_get_max_threads(); +#endif + if (m->fd < 0) return; @@ -458,7 +478,9 @@ static void close_map(struct map *m) } if (m->use_rowio) { - cache_release(&m->cache); + for (int i = 0; i < threads; i++) + cache_release(&m->caches[i]); + G_free(&m->caches); m->use_rowio = 0; } } diff --git a/raster/r.mapcalc/map3.c b/raster/r.mapcalc/map3.c index 374f6a76109..ed543870d66 100644 --- a/raster/r.mapcalc/map3.c +++ b/raster/r.mapcalc/map3.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include @@ -619,8 +623,12 @@ int open_output_map(const char *name, int res_type) void put_map_row(int fd, void *buf, int res_type) { void *handle = omaps[fd]; + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif - write_row(handle, buf, res_type, current_depth, current_row); + write_row(handle, buf, res_type, current_depth, current_row[tid]); } void close_output_map(int fd) diff --git a/raster/r.mapcalc/mapcalc.h b/raster/r.mapcalc/mapcalc.h index 75d6b3da020..a2c6947c641 100644 --- a/raster/r.mapcalc/mapcalc.h +++ b/raster/r.mapcalc/mapcalc.h @@ -2,6 +2,9 @@ #define _MAPCALC_H_ /****************************************************************************/ +#if defined(_OPENMP) +#include +#endif #include diff --git a/raster/r.mapcalc/xarea.c b/raster/r.mapcalc/xarea.c index 74bd1392c3c..6bbe605ba6a 100644 --- a/raster/r.mapcalc/xarea.c +++ b/raster/r.mapcalc/xarea.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include "globals.h" @@ -10,6 +14,11 @@ area() area of a cell in square meters int f_area(int argc, const int *argt, void **args) { + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + DCELL *res = args[0]; int i; static int row = -1; @@ -21,11 +30,11 @@ int f_area(int argc, const int *argt, void **args) if (argt[0] != DCELL_TYPE) return E_RES_TYPE; - if (row != current_row) { + if (row != current_row[tid]) { if (row == -1) G_begin_cell_area_calculations(); - row = current_row; + row = current_row[tid]; cell_area = G_area_of_cell_at_row(row); } diff --git a/raster/r.mapcalc/xcoor.c b/raster/r.mapcalc/xcoor.c index 5f783d7e90f..1a209a25281 100644 --- a/raster/r.mapcalc/xcoor.c +++ b/raster/r.mapcalc/xcoor.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include "globals.h" @@ -35,6 +39,11 @@ int f_x(int argc, const int *argt, void **args) int f_y(int argc, const int *argt, void **args) { + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + DCELL *res = args[0]; DCELL y; int i; @@ -45,7 +54,7 @@ int f_y(int argc, const int *argt, void **args) if (argt[0] != DCELL_TYPE) return E_RES_TYPE; - y = Rast_row_to_northing(current_row + 0.5, ¤t_region2); + y = Rast_row_to_northing(current_row[tid] + 0.5, ¤t_region2); for (i = 0; i < columns; i++) res[i] = y; diff --git a/raster/r.mapcalc/xcoor3.c b/raster/r.mapcalc/xcoor3.c index 8c8906b890e..21fa5a381b0 100644 --- a/raster/r.mapcalc/xcoor3.c +++ b/raster/r.mapcalc/xcoor3.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include "globals.h" @@ -36,6 +40,11 @@ int f_x(int argc, const int *argt, void **args) int f_y(int argc, const int *argt, void **args) { + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + RASTER3D_Region *window = ¤t_region3; DCELL *res = args[0]; DCELL y; @@ -47,7 +56,7 @@ int f_y(int argc, const int *argt, void **args) if (argt[0] != DCELL_TYPE) return E_RES_TYPE; - y = window->north - (current_row + 0.5) * window->ns_res; + y = window->north - (current_row[tid] + 0.5) * window->ns_res; for (i = 0; i < columns; i++) res[i] = y; diff --git a/raster/r.mapcalc/xrowcol.c b/raster/r.mapcalc/xrowcol.c index 5f95d63e618..4f332589e79 100644 --- a/raster/r.mapcalc/xrowcol.c +++ b/raster/r.mapcalc/xrowcol.c @@ -1,3 +1,7 @@ +#if defined(_OPENMP) +#include +#endif + #include #include #include "globals.h" @@ -32,8 +36,13 @@ int f_col(int argc, const int *argt, void **args) int f_row(int argc, const int *argt, void **args) { + int tid = 0; +#if defined(_OPENMP) + tid = omp_get_thread_num(); +#endif + CELL *res = args[0]; - int row = current_row + 1; + int row = current_row[tid] + 1; int i; if (argc > 0) From 72a0970bad36486e6af98b8b64ee9f24e8c48e2d Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Sat, 22 Jun 2024 10:15:07 -0400 Subject: [PATCH 02/25] fix the segmentfault problems --- raster/r.mapcalc/evaluate.c | 63 ++++++++++++++++++++++++++++++++----- raster/r.mapcalc/main.c | 6 +++- 2 files changed, 61 insertions(+), 8 deletions(-) diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c index 20dc2f5d325..d053edd4f79 100644 --- a/raster/r.mapcalc/evaluate.c +++ b/raster/r.mapcalc/evaluate.c @@ -90,6 +90,34 @@ static void set_buf(expression *e, void **buf) e->buf = buf; } +static void free_buf(expression *e) +{ + int threads = 1; +#if defined(_OPENMP) + threads = omp_get_max_threads(); +#endif + + for (int t = 0; t < threads; t++) { + G_free(e->buf[t]); + e->buf[t] = NULL; + } + G_free(e->buf); + e->buf = NULL; +} + +static void free_argv(expression *e) +{ + int i; + + for (i = 1; i <= e->data.func.argc; i++) { + free_buf(e->data.func.args[i]); + e->data.func.args[i]->buf = NULL; + } + + G_free(e->data.func.argv); + e->data.func.argv = NULL; +} + /****************************************************************************/ static void initialize_constant(expression *e) @@ -340,7 +368,9 @@ void execute(expr_list *ee) { int verbose = isatty(2); expr_list *l; - int count, n; + expression **exp_arr; + int count, n, i; + int num_exprs = 0; int threads = 1; #if defined(_OPENMP) @@ -367,13 +397,19 @@ void execute(expr_list *ee) G_fatal_error(_("output map <%s> exists. To overwrite, " "use the --overwrite flag"), var); + num_exprs++; } + /* Create a array of expreesion and stored it in heap */ + exp_arr = G_malloc(num_exprs * sizeof(struct expression *)); + /* Parse each expression and extract all raster maps */ - for (l = ee; l; l = l->next) { + l = ee; + for (i = 0; i < num_exprs; i++) { expression *e = l->exp; - extract_maps(e); + exp_arr[i] = e; + l = l->next; } /* Set the region from the input maps */ @@ -385,8 +421,9 @@ void execute(expr_list *ee) setup_region(); /* Parse each expression and initialize the maps, buffers and variables */ - for (l = ee; l; l = l->next) { - expression *e = l->exp; + + for (i = 0; i < num_exprs; i++) { + expression *e = exp_arr[i]; const char *var; expression *val; @@ -416,8 +453,8 @@ void execute(expr_list *ee) tid = omp_get_thread_num(); #endif current_row[tid] = row; - for (l = ee; l != NULL; l = l->next) { - expression *e = l->exp; + for (i = 0; i < num_exprs; i++) { + expression *e = exp_arr[i]; int fd; evaluate(e); @@ -472,6 +509,18 @@ void execute(expr_list *ee) } G_unset_error_routine(); + + /* Free the memory and make it unreachable */ + G_free(current_row); + for (i = 0; i < num_exprs; i++) { + expression *e = exp_arr[i]; + free_buf(e); + if (e->type == expr_type_function) + free_argv(e); + } + G_free(exp_arr); + current_row = NULL; + exp_arr = NULL; } void describe_maps(FILE *fp, expr_list *ee) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 25926e865c3..81e3ba18822 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -133,6 +133,8 @@ int main(int argc, char **argv) p[2] = NULL; argv = p; argc = 2; + G_free(p); + p = NULL; } if (G_parser(argc, argv)) @@ -195,13 +197,15 @@ int main(int argc, char **argv) G_message(_("Computing in parallel, number of threads: %d"), omp_get_max_threads()); #else - G_message(_("Number of threads: 1")); + G_message(_("Computing in serial (no OpenMP support)")); #endif execute(result); post_exec(); all_ok = 1; + G_free(nprocs->answer); + if (floating_point_exception_occurred) { G_warning(_("Floating point error(s) occurred in the calculation")); all_ok = 0; From 420cf19f5bee79e5f82b7a0c1cd5b69a85c82a66 Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Sat, 22 Jun 2024 10:16:23 -0400 Subject: [PATCH 03/25] add tests for parallelization --- .../testsuite/test_r_mapcalc_parallel.py | 270 ++++++++++++++++++ 1 file changed, 270 insertions(+) create mode 100644 raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py diff --git a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py new file mode 100644 index 00000000000..9b44632d41a --- /dev/null +++ b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py @@ -0,0 +1,270 @@ +from grass.gunittest.case import TestCase +from grass.gunittest.main import test +from grass.gunittest.gmodules import SimpleModule + +cell_seed_500 = """\ +north: 20 +south: 10 +east: 25 +west: 15 +rows: 10 +cols: 10 +121 12 183 55 37 96 138 117 182 40 +157 70 115 1 149 125 42 193 108 24 +83 66 82 84 186 182 179 122 67 113 +151 93 144 173 128 196 61 125 64 193 +180 175 14 41 44 27 165 27 90 60 +97 57 12 104 98 13 87 24 83 107 +174 133 146 114 115 60 78 154 49 130 +55 138 144 25 32 58 47 137 139 32 +143 193 155 190 131 124 87 81 160 154 +56 45 48 66 9 182 69 12 154 19 +""" + +dcell_seed_600 = """\ +north: 20 +south: 10 +east: 25 +west: 15 +rows: 10 +cols: 10 +130.790433856418332 101.3319248101491041 33.5781271447787759 37.4064724824657944 98.2794723130458152 73.9118866262841863 185.9530433718733775 74.5210037729812882 166.1178416001017695 90.9915650902159996 +109.2478664232956334 25.6499350759712215 150.9024447059825036 125.7544119036241312 66.7235333366722614 167.9375729129454271 123.1009291055983965 12.0922254083554606 59.389026967287819 113.2843489528100349 +40.0044184023145277 135.8273774212801186 71.6737798852435049 191.6223505280372876 4.1546013811569615 143.3082794522489962 177.043829294835632 115.0300571162354402 141.8985452774071234 127.8949967123638061 +93.2842559637482793 9.7471880052423856 118.1216452002055632 158.1474162140586088 67.2957262519499437 3.6524546350146849 147.0965842667525862 37.060628529871579 47.3408278816968959 66.2219633495724054 +175.5638637866295539 67.1399023507611901 162.2058392782793703 198.1586789345953719 36.474049475167746 49.2589028048889617 112.1169663235969836 22.0227597984432535 95.9169228571662131 86.7470895014531322 +93.5401613888204935 193.7821104138942587 193.8286564351004699 3.2623643889134684 94.6955247357847725 25.7099391122614307 155.592251526775442 25.3392337002970294 48.3979699868663005 99.6836079482556272 +104.16296861365457 190.7865884377180805 6.2841805474238619 49.3731395705159528 100.1903962703459285 116.927654961282343 19.8626348109264264 40.9693022766258466 81.6500759554420057 169.2220572316770131 +118.8112518721558217 55.8955021401724039 112.9150308215961331 62.6399760484719081 85.400498505854145 191.0144187084912062 124.2128169358724534 167.9341741649760706 170.6149695243870781 158.3034517206661462 +130.0453795775294736 64.1996403829061535 62.9317494959142465 175.1909990236256931 122.9624852869890361 79.9546265736285733 9.6594013716963367 114.0611338072915544 11.9371167643030809 186.9121199748369122 +3.2891990250261536 30.9245408751958379 46.4021422454598991 104.2378950097200203 47.424093232347019 73.4801303522840499 22.4778583078695213 132.870185207462697 48.1666164169167388 100.5504714442693057 +""" + +fcell_seed_700 = """\ +north: 20 +south: 10 +east: 25 +west: 15 +rows: 10 +cols: 10 +146.756378 192.682159 2.644822 147.270462 62.178818 192.668198 94.320778 107.710426 98.319664 114.444504 +12.995321 18.026272 151.590958 5.249451 197.266708 103.663635 115.424088 28.01062 78.555168 62.912098 +164.053619 154.652039 98.536011 44.601639 85.322289 168.383957 44.93845 128.62262 89.910591 107.242188 +111.182487 63.080284 177.791473 47.439354 42.451859 72.396568 170.597778 170.622742 141.88858 105.126854 +120.76828 148.581085 42.124866 56.432236 164.652176 98.094009 60.741329 66.286987 187.847427 160.120056 +50.530689 179.090652 138.114014 138.629211 193.147903 172.861481 133.72728 108.720459 103.508438 28.81559 +39.653179 101.948265 35.744762 25.570076 78.767021 154.600616 144.907684 82.370148 116.378654 18.218494 +35.587288 66.534409 65.744408 186.476959 137.081116 151.379272 48.261463 8.323328 130.432739 53.346546 +152.67189 15.512391 146.049072 185.276245 34.417141 127.522453 124.54998 52.08218 167.141342 87.771118 +69.0522 43.57811 63.15279 68.677063 74.202805 97.429077 167.123199 19.892767 120.593437 190.960815 +""" + +THREADS = 4 + + +class TestRandFunction(TestCase): + # TODO: replace by unified handing of maps + to_remove = [] + + @classmethod + def setUpClass(cls): + cls.use_temp_region() + cls.runModule("g.region", n=20, s=10, e=25, w=15, res=1) + + @classmethod + def tearDownClass(cls): + cls.del_temp_region() + if cls.to_remove: + cls.runModule( + "g.remove", flags="f", type="raster", name=",".join(cls.to_remove) + ) + + def rinfo_contains_number(self, raster, number): + """Test that r.info standard output for raster contains a given number + + To be used in test methods for testing presence of a given number. + """ + rinfo = SimpleModule("r.info", map=raster) + self.runModule(rinfo) + self.assertIn(str(number), rinfo.outputs.stdout) + + def test_seed_not_required(self): + """Test that seed is not required when rand() is not used""" + self.assertModule("r.mapcalc", expression="nonrand_cell = 200", nprocs=THREADS) + self.to_remove.append("nonrand_cell") + + +# TODO: add more expressions +# TODO: add tests with prepared data + + +class TestBasicOperations(TestCase): + # TODO: replace by unified handing of maps + to_remove = [] + + @classmethod + def setUpClass(cls): + cls.use_temp_region() + cls.runModule("g.region", n=20, s=10, e=25, w=15, res=1) + + @classmethod + def tearDownClass(cls): + cls.del_temp_region() + if cls.to_remove: + cls.runModule( + "g.remove", + flags="f", + type="raster", + name=",".join(cls.to_remove), + verbose=True, + ) + + def test_difference_of_the_same_map_double(self): + """Test zero difference of map with itself""" + self.runModule("r.mapcalc", flags="s", expression="a = rand(1.0, 200)") + self.to_remove.append("a") + self.assertModule("r.mapcalc", expression="diff_a_a = a - a", nprocs=THREADS) + self.to_remove.append("diff_a_a") + self.assertRasterMinMax("diff_a_a", refmin=0, refmax=0) + + def test_difference_of_the_same_map_float(self): + """Test zero difference of map with itself""" + self.runModule("r.mapcalc", flags="s", expression="af = rand(float(1), 200)") + self.to_remove.append("af") + self.assertModule( + "r.mapcalc", expression="diff_af_af = af - af", nprocs=THREADS + ) + self.to_remove.append("diff_af_af") + self.assertRasterMinMax("diff_af_af", refmin=0, refmax=0) + + def test_difference_of_the_same_map_int(self): + """Test zero difference of map with itself""" + self.runModule("r.mapcalc", flags="s", expression="ai = rand(1, 200)") + self.to_remove.append("ai") + self.assertModule( + "r.mapcalc", expression="diff_ai_ai = ai - ai", nprocs=THREADS + ) + self.to_remove.append("diff_ai_ai") + self.assertRasterMinMax("diff_ai_ai", refmin=0, refmax=0) + + def test_difference_of_the_same_expression(self): + """Test zero difference of two same expressions""" + self.assertModule( + "r.mapcalc", + expression="diff_e_e = 3 * x() * y() - 3 * x() * y()", + nprocs=THREADS, + ) + self.to_remove.append("diff_e_e") + self.assertRasterMinMax("diff_e_e", refmin=0, refmax=0) + + def test_nrows_ncols_sum(self): + """Test if sum of nrows and ncols matches one + expected from current region settings""" + self.assertModule( + "r.mapcalc", + expression="nrows_ncols_sum = nrows() + ncols()", + nprocs=THREADS, + ) + self.to_remove.append("nrows_ncols_sum") + self.assertRasterMinMax("nrows_ncols_sum", refmin=20, refmax=20) + + +class TestRegionOperations(TestCase): + # TODO: replace by unified handing of maps + to_remove = [] + + @classmethod + def setUpClass(cls): + cls.use_temp_region() + cls.runModule("g.region", n=30, s=15, e=30, w=15, res=5) + cls.runModule( + "r.mapcalc", expression="test_region_1 = 1", seed=1, nprocs=THREADS + ) + cls.runModule("g.region", n=25, s=10, e=25, w=10, res=5) + cls.runModule( + "r.mapcalc", expression="test_region_2 = 2", seed=1, nprocs=THREADS + ) + cls.runModule("g.region", n=20, s=5, e=20, w=5, res=1) + cls.runModule( + "r.mapcalc", expression="test_region_3 = 3", seed=1, nprocs=THREADS + ) + + cls.to_remove.append("test_region_1") + cls.to_remove.append("test_region_2") + cls.to_remove.append("test_region_3") + + @classmethod + def tearDownClass(cls): + cls.del_temp_region() + if cls.to_remove: + cls.runModule( + "g.remove", + flags="f", + type="raster", + name=",".join(cls.to_remove), + verbose=True, + ) + + def test_union(self): + """Test the union region option""" + self.assertModule( + "r.mapcalc", + region="union", + seed=1, + expression="test_region_4 = test_region_1 + test_region_2 + test_region_3", + nprocs=THREADS, + ) + self.to_remove.append("test_region_4") + + self.assertModuleKeyValue( + "r.info", + map="test_region_4", + flags="gr", + reference=dict( + min=6, + max=6, + cells=625, + north=30, + south=5, + west=5, + east=30, + nsres=1, + ewres=1, + ), + precision=0.01, + sep="=", + ) + + def test_intersect(self): + """Test the intersect region option""" + self.assertModule( + "r.mapcalc", + region="intersect", + seed=1, + expression="test_region_5 = test_region_1 + test_region_2 + test_region_3", + nprocs=THREADS, + ) + self.to_remove.append("test_region_5") + + self.assertModuleKeyValue( + "r.info", + map="test_region_5", + flags="gr", + reference=dict( + min=6, + max=6, + cells=25, + north=20, + south=15, + west=15, + east=20, + nsres=1, + ewres=1, + ), + precision=0.01, + sep="=", + ) + + +if __name__ == "__main__": + test() From e7fd5428edce6ebf09f2ab5bcfa67b41b72e2254 Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Tue, 20 May 2025 14:18:26 -0400 Subject: [PATCH 04/25] Update raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- .../testsuite/test_r_mapcalc_parallel.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py index 9b44632d41a..8fa0a2712ee 100644 --- a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py +++ b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py @@ -250,17 +250,17 @@ def test_intersect(self): "r.info", map="test_region_5", flags="gr", - reference=dict( - min=6, - max=6, - cells=25, - north=20, - south=15, - west=15, - east=20, - nsres=1, - ewres=1, - ), + reference={ + "min": 6, + "max": 6, + "cells": 25, + "north": 20, + "south": 15, + "west": 15, + "east": 20, + "nsres": 1, + "ewres": 1, + }, precision=0.01, sep="=", ) From 5a70b4546659b2de49fd2cdea946869e7e49c528 Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Tue, 20 May 2025 14:18:33 -0400 Subject: [PATCH 05/25] Update raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- .../testsuite/test_r_mapcalc_parallel.py | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py index 8fa0a2712ee..71b42b3b388 100644 --- a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py +++ b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py @@ -220,17 +220,17 @@ def test_union(self): "r.info", map="test_region_4", flags="gr", - reference=dict( - min=6, - max=6, - cells=625, - north=30, - south=5, - west=5, - east=30, - nsres=1, - ewres=1, - ), + reference={ + "min": 6, + "max": 6, + "cells": 625, + "north": 30, + "south": 5, + "west": 5, + "east": 30, + "nsres": 1, + "ewres": 1, + }, precision=0.01, sep="=", ) From 8c497d095105e152bbc2a9e65f0a131cc319b66d Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Tue, 20 May 2025 20:04:34 -0400 Subject: [PATCH 06/25] remove redundant parts --- raster/r.mapcalc/main.c | 18 ++++++++++-------- raster/r.mapcalc/map.c | 29 +++++------------------------ 2 files changed, 15 insertions(+), 32 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 65d598d5441..410a1f0e39a 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -121,16 +121,13 @@ int main(int argc, char **argv) nprocs = G_define_standard_option(G_OPT_M_NPROCS); + char **p = G_malloc(3 * sizeof(char *)); if (argc == 1) { - char **p = G_malloc(3 * sizeof(char *)); - p[0] = argv[0]; p[1] = G_store("file=-"); p[2] = NULL; argv = p; argc = 2; - G_free(p); - p = NULL; } if (G_parser(argc, argv)) @@ -190,17 +187,22 @@ int main(int argc, char **argv) threads = atoi(nprocs->answer); #if defined(_OPENMP) omp_set_num_threads(threads); - G_message(_("Computing in parallel, number of threads: %d"), - omp_get_max_threads()); + // G_message(_("Computing in parallel, number of threads: %d"), + // omp_get_max_threads()); #else - G_message(_("Computing in serial (no OpenMP support)")); + // G_message(_("Computing in serial (no OpenMP support)")); #endif execute(result); post_exec(); all_ok = 1; - G_free(nprocs->answer); + // Free nprocs->answer if it is not the default value, "1" + if (threads > 1) { + G_free(nprocs->answer); + } + G_free(p); + p = NULL; if (floating_point_exception_occurred) { G_warning(_("Floating point error(s) occurred in the calculation")); diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 9b7bf1bb88d..9537c1cb153 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -179,31 +179,12 @@ static void *cache_get_raw(struct row_cache *cache, int row, int data_type) return sub->buf[0]; } - tmp = G_alloca(cache->nrows * sizeof(void *)); - memcpy(tmp, sub->buf, cache->nrows * sizeof(void *)); - vtmp = G_alloca(cache->nrows); - memcpy(vtmp, sub->valid, cache->nrows); - - i = (i < 0) ? 0 : cache->nrows - 1; - newrow = row - i; - - for (j = 0; j < cache->nrows; j++) { - int r = newrow + j; - int k = r - sub->row; - int l = (k + cache->nrows) % cache->nrows; - - sub->buf[j] = tmp[l]; - sub->valid[j] = k >= 0 && k < cache->nrows && vtmp[l]; + else { + i = (i < 0) ? 0 : cache->nrows - 1; + newrow = row - i; + read_row(cache->fd, sub->buf[i], row, data_type); + return sub->buf[i]; } - - sub->row = newrow; - G_freea(tmp); - G_freea(vtmp); - - read_row(cache->fd, sub->buf[i], row, data_type); - sub->valid[i] = 1; - - return sub->buf[i]; } static void cache_get(struct row_cache *cache, void *buf, int row, int res_type) From 14ca37641c0e0f24d1d6975e4de49e3c41075d2f Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Wed, 21 May 2025 14:21:53 -0400 Subject: [PATCH 07/25] fix malloc errors related threading --- .../r.mapcalc/benchmark/benchmark_rmapcalc.py | 69 +++++++++++++ raster/r.mapcalc/evaluate.c | 7 ++ raster/r.mapcalc/main.c | 11 ++- raster/r.mapcalc/map.c | 3 +- .../testsuite/test_r_mapcalc_parallel.py | 97 +++++++++++++++++++ 5 files changed, 182 insertions(+), 5 deletions(-) create mode 100644 raster/r.mapcalc/benchmark/benchmark_rmapcalc.py diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py new file mode 100644 index 00000000000..3c993d849e4 --- /dev/null +++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py @@ -0,0 +1,69 @@ +"""Benchmarking of r.horizon +point mode, one direction + +@author Chung-Yuan Liang, 2024 +""" + +from grass.exceptions import CalledModuleError +from grass.pygrass.modules import Module + +import grass.benchmark as bm + + +def main(): + results = [] + metrics = ["time", "speedup", "efficiency"] + mapsizes = [10e6, 20e6, 40e6, 80e6] + + # run benchmarks + + for mapsize in mapsizes: + benchmark( + size=int(mapsize**0.5), + step=0, + label=f"r.mapcalc{int(mapsize / 1e6)}M", + results=results, + ) + + # plot results + + for metric in metrics: + bm.nprocs_plot( + results, + title=f"r.mapcalc {metric}", + metric=metric, + ) + + +def benchmark(size, step, label, results): + map1 = "benchmark_r_mapcalc_map1" + map2 = "benchmark_r_mapcalc_map2" + output = "benchmark_r_mapcalc" + + generate_map(rows=size, cols=size, fname=map1) + generate_map(rows=size, cols=size, fname=map2) + module = Module( + "r.mapcalc", + expression=f"{output}=({map1} + {map2})/2", + overwrite=True, + ) + + results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=8, repeat=3)) + Module( + "g.remove", quiet=True, flags="f", type="raster", pattern="benchmark_r_mapcalc*" + ) + + +def generate_map(rows, cols, fname): + Module("g.region", flags="p", rows=rows, cols=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, overwrite=True) + except CalledModuleError: + print("r.surf.fractal fails, using r.random.surface instead...") + Module("r.random.surface", output=fname, overwrite=True) + + +if __name__ == "__main__": + main() diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c index d053edd4f79..e898ea287d3 100644 --- a/raster/r.mapcalc/evaluate.c +++ b/raster/r.mapcalc/evaluate.c @@ -375,6 +375,13 @@ void execute(expr_list *ee) #if defined(_OPENMP) threads = omp_get_max_threads(); + if ((threads > rows) && (threads > 1)) { + threads = rows; + omp_set_num_threads(threads); + G_message(_("The number of rows is less than the number of threads. \ + Set the number of threads to be the same as the rows = %d"), + threads); + } #endif current_row = (int *)G_malloc(sizeof(int) * threads); diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 410a1f0e39a..1d22e7302e5 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -187,11 +187,12 @@ int main(int argc, char **argv) threads = atoi(nprocs->answer); #if defined(_OPENMP) omp_set_num_threads(threads); - // G_message(_("Computing in parallel, number of threads: %d"), - // omp_get_max_threads()); -#else - // G_message(_("Computing in serial (no OpenMP support)")); + // G_message(_("Compute in paraell. Number of threads set to %d."), + // threads); +// #else +// G_message(_("Compute in searial. OpenMP is not available.")); #endif + execute(result); post_exec(); @@ -209,6 +210,8 @@ int main(int argc, char **argv) all_ok = 0; } + printf("all_ok = %d\n", all_ok); + return all_ok ? EXIT_SUCCESS : EXIT_FAILURE; } diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 9537c1cb153..73ed9b398c4 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -460,7 +460,8 @@ static void close_map(struct map *m) if (m->use_rowio) { for (int i = 0; i < threads; i++) cache_release(&m->caches[i]); - G_free(&m->caches); + if (threads > 1) + G_free(m->caches); m->use_rowio = 0; } } diff --git a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py index 71b42b3b388..c635785daba 100644 --- a/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py +++ b/raster/r.mapcalc/testsuite/test_r_mapcalc_parallel.py @@ -93,6 +93,103 @@ def test_seed_not_required(self): self.assertModule("r.mapcalc", expression="nonrand_cell = 200", nprocs=THREADS) self.to_remove.append("nonrand_cell") + def test_seed_required(self): + """Test that seed is required when rand() is used + + This test can, and probably should, generate an error message. + """ + self.assertModuleFail( + "r.mapcalc", expression="rand_x = rand(1, 200)", nprocs=THREADS + ) + # TODO: assert map not exists but it would be handy here + # TODO: test that error message was generated + + def test_seed_cell(self): + """Test given seed with CELL against reference map""" + seed = 500 + self.runModule( + "r.in.ascii", input="-", stdin=cell_seed_500, output="rand_cell_ref" + ) + self.to_remove.append("rand_cell_ref") + self.assertModule( + "r.mapcalc", + seed=seed, + expression="rand_cell = rand(1, 200)", + nprocs=THREADS, + ) + self.to_remove.append("rand_cell") + # this assert is using r.mapcalc but we are testing different + # functionality than used by assert + self.assertRastersNoDifference( + actual="rand_cell", reference="rand_cell_ref", precision=0 + ) + self.rinfo_contains_number("rand_cell", seed) + + def test_seed_dcell(self): + """Test given seed with DCELL against reference map""" + seed = 600 + self.runModule( + "r.in.ascii", input="-", stdin=dcell_seed_600, output="rand_dcell_ref" + ) + self.to_remove.append("rand_dcell_ref") + self.assertModule( + "r.mapcalc", + seed=seed, + expression="rand_dcell = rand(1.0, 200.0)", + nprocs=THREADS, + ) + self.to_remove.append("rand_dcell") + # this assert is using r.mapcalc but we are testing different + # functionality than used by assert + self.assertRastersNoDifference( + actual="rand_dcell", reference="rand_dcell_ref", precision=0.00000000000001 + ) + self.rinfo_contains_number("rand_dcell", seed) + + def test_seed_fcell(self): + """Test given seed with FCELL against reference map""" + seed = 700 + self.runModule( + "r.in.ascii", input="-", stdin=fcell_seed_700, output="rand_fcell_ref" + ) + self.to_remove.append("rand_fcell_ref") + self.assertModule( + "r.mapcalc", + seed=seed, + expression="rand_fcell = rand(float(1), 200)", + nprocs=THREADS, + ) + self.to_remove.append("rand_fcell") + # this assert is using r.mapcalc but we are testing different + # functionality than used by assert + self.assertRastersNoDifference( + actual="rand_fcell", reference="rand_fcell_ref", precision=0.000001 + ) + self.rinfo_contains_number("rand_fcell", seed) + + def test_auto_seed(self): + """Test that two runs with -s does not give same maps""" + self.assertModule( + "r.mapcalc", + flags="s", + expression="rand_auto_1 = rand(1., 2)", + nprocs=THREADS, + ) + self.to_remove.append("rand_auto_1") + self.assertModule( + "r.mapcalc", + flags="s", + expression="rand_auto_2 = rand(1., 2)", + nprocs=THREADS, + ) + self.to_remove.append("rand_auto_2") + self.assertRastersDifference( + "rand_auto_1", + "rand_auto_2", + statistics={"min": -1, "max": 1, "mean": 0}, + precision=0.5, + ) # low precision, we have few cells + # TODO: add more expressions # TODO: add tests with prepared data From d0ca65e77c060e9424d77881dac598018d76c9c0 Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Thu, 22 May 2025 10:14:50 -0400 Subject: [PATCH 08/25] use multiple file descriptors instead of locks --- .../r.mapcalc/benchmark/benchmark_rmapcalc.py | 2 +- raster/r.mapcalc/evaluate.c | 66 ++++++++++++------- raster/r.mapcalc/expression.h | 2 +- raster/r.mapcalc/main.c | 2 - raster/r.mapcalc/map.c | 31 +++++---- 5 files changed, 60 insertions(+), 43 deletions(-) diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py index 3c993d849e4..ea0c8310daa 100644 --- a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py +++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py @@ -13,7 +13,7 @@ def main(): results = [] metrics = ["time", "speedup", "efficiency"] - mapsizes = [10e6, 20e6, 40e6, 80e6] + mapsizes = [10e6, 20e6] # , 40e6, 80e6] # run benchmarks diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c index e898ea287d3..f098958508c 100644 --- a/raster/r.mapcalc/evaluate.c +++ b/raster/r.mapcalc/evaluate.c @@ -134,8 +134,15 @@ static void initialize_map(expression *e) { allocate_buf(e); - e->data.map.idx = open_map(e->data.map.name, e->data.map.mod, - e->data.map.row, e->data.map.col); + int threads = 1; +#if defined(_OPENMP) + threads = omp_get_max_threads(); +#endif + e->data.map.idx = G_malloc(threads * sizeof(int)); + for (int t = 0; t < threads; t++) { + e->data.map.idx[t] = open_map(e->data.map.name, e->data.map.mod, + e->data.map.row, e->data.map.col); + } } static void initialize_function(expression *e) @@ -244,11 +251,21 @@ static void evaluate_map(expression *e) #if defined(_OPENMP) tid = omp_get_thread_num(); #endif - - get_map_row(e->data.map.idx, e->data.map.mod, + get_map_row(e->data.map.idx[tid], e->data.map.mod, current_depth + e->data.map.depth, current_row[tid] + e->data.map.row, e->data.map.col, e->buf[tid], e->res_type); + + // printf("reading: tid %d, map_idx %d, current_d %d, map_d %d, current_r + // %d, map_row %d, map_col %d, value %f\n", + // tid, + // e->data.map.idx, + // current_depth, + // e->data.map.depth, + // current_row[tid], + // e->data.map.row, + // e->data.map.col, + // ((double *)e->buf[tid])[0]); } static void evaluate_function(expression *e) @@ -373,18 +390,6 @@ void execute(expr_list *ee) int num_exprs = 0; int threads = 1; -#if defined(_OPENMP) - threads = omp_get_max_threads(); - if ((threads > rows) && (threads > 1)) { - threads = rows; - omp_set_num_threads(threads); - G_message(_("The number of rows is less than the number of threads. \ - Set the number of threads to be the same as the rows = %d"), - threads); - } -#endif - current_row = (int *)G_malloc(sizeof(int) * threads); - exprs = ee; G_add_error_handler(error_handler, NULL); @@ -446,11 +451,22 @@ void execute(expr_list *ee) setup_maps(); +#if defined(_OPENMP) + threads = omp_get_max_threads(); + if ((threads > rows) && (threads > 1)) { + threads = rows; + omp_set_num_threads(threads); + G_message(_("The number of rows is less than the number of threads. \ + Set the number of threads to be the same as the rows = %d"), + threads); + } +#endif + current_row = (int *)G_malloc(sizeof(int) * threads); count = rows * depths; n = 0; for (current_depth = 0; current_depth < depths; current_depth++) { -#pragma omp parallel for default(shared) schedule(static, 1) ordered +#pragma omp parallel for default(shared) schedule(static, 1) private(i) ordered for (int row = 0; row < rows; row++) { if (verbose) G_percent(n, count, 2); @@ -463,20 +479,20 @@ void execute(expr_list *ee) for (i = 0; i < num_exprs; i++) { expression *e = exp_arr[i]; int fd; - evaluate(e); #pragma omp ordered { if (e->type == expr_type_binding) { fd = e->data.bind.fd; + // printf("binding: fd %d, tid %d, row %d, value %f\n", + // fd, tid, + // row, ((double *)e->buf[tid])[0]); put_map_row(fd, e->buf[tid], e->res_type); } } } -#pragma omp critical - { - n++; - } +#pragma omp atomic update + n++; } } @@ -505,11 +521,11 @@ void execute(expr_list *ee) if (val->type == expr_type_map) { if (val->data.map.mod == 'M') { - copy_cats(var, val->data.map.idx); - copy_colors(var, val->data.map.idx); + copy_cats(var, val->data.map.idx[0]); + copy_colors(var, val->data.map.idx[0]); } - copy_history(var, val->data.map.idx); + copy_history(var, val->data.map.idx[0]); } else create_history(var, val); diff --git a/raster/r.mapcalc/expression.h b/raster/r.mapcalc/expression.h index e8ec3bea5d5..f55fffd3650 100644 --- a/raster/r.mapcalc/expression.h +++ b/raster/r.mapcalc/expression.h @@ -27,7 +27,7 @@ typedef struct expr_data_map { const char *name; int mod; int row, col, depth; - int idx; + int *idx; /* array to store fds for multi-threads*/ } expr_data_map; typedef struct expr_data_func { diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 1d22e7302e5..1a54b07b162 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -210,8 +210,6 @@ int main(int argc, char **argv) all_ok = 0; } - printf("all_ok = %d\n", all_ok); - return all_ok ? EXIT_SUCCESS : EXIT_FAILURE; } diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 73ed9b398c4..2e5c48b5c6d 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -426,6 +426,9 @@ static void read_map(struct map *m, void *buf, int res_type, int row, int col) else read_row(m->fd, buf, row, res_type); +#ifdef _OPENMP + tid = omp_get_thread_num(); +#endif if (col) column_shift(buf, res_type, col); } @@ -536,25 +539,25 @@ int open_map(const char *name, int mod, int row, int col) break; } - for (i = 0; i < num_maps; i++) { - m = &maps[i]; + // for (i = 0; i < num_maps; i++) { + // m = &maps[i]; - if (strcmp(m->name, name) != 0 || strcmp(m->mapset, mapset) != 0) - continue; + // if (strcmp(m->name, name) != 0 || strcmp(m->mapset, mapset) != 0) + // continue; - if (row < m->min_row) - m->min_row = row; - if (row > m->max_row) - m->max_row = row; + // if (row < m->min_row) + // m->min_row = row; + // if (row > m->max_row) + // m->max_row = row; - if (use_cats && !m->have_cats) - init_cats(m); + // if (use_cats && !m->have_cats) + // init_cats(m); - if (use_colors && !m->have_colors) - init_colors(m); + // if (use_colors && !m->have_colors) + // init_colors(m); - return i; - } + // return i; + // } if (num_maps >= max_maps) { max_maps += 10; From 5b3184d3c7113ebf488bfb28be5e5dc406528a89 Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Thu, 22 May 2025 14:37:51 -0400 Subject: [PATCH 09/25] disable parallelization for not implemented functions --- .../r.mapcalc/benchmark/benchmark_rmapcalc.py | 2 +- raster/r.mapcalc/evaluate.c | 32 ++++++++----------- raster/r.mapcalc/main.c | 25 ++++++++++++--- raster/r.mapcalc/map.c | 20 ------------ 4 files changed, 36 insertions(+), 43 deletions(-) diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py index ea0c8310daa..4d9e5cec60e 100644 --- a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py +++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py @@ -44,7 +44,7 @@ def benchmark(size, step, label, results): generate_map(rows=size, cols=size, fname=map2) module = Module( "r.mapcalc", - expression=f"{output}=({map1} + {map2})/2", + expression=f"{output}=({map1} + {map2} - 2*{map2} + {map1}*{map2} - {map1}/{map2})/2", overwrite=True, ) diff --git a/raster/r.mapcalc/evaluate.c b/raster/r.mapcalc/evaluate.c index f098958508c..a0bac5f9cc1 100644 --- a/raster/r.mapcalc/evaluate.c +++ b/raster/r.mapcalc/evaluate.c @@ -255,17 +255,6 @@ static void evaluate_map(expression *e) current_depth + e->data.map.depth, current_row[tid] + e->data.map.row, e->data.map.col, e->buf[tid], e->res_type); - - // printf("reading: tid %d, map_idx %d, current_d %d, map_d %d, current_r - // %d, map_row %d, map_col %d, value %f\n", - // tid, - // e->data.map.idx, - // current_depth, - // e->data.map.depth, - // current_row[tid], - // e->data.map.row, - // e->data.map.col, - // ((double *)e->buf[tid])[0]); } static void evaluate_function(expression *e) @@ -383,7 +372,7 @@ static void error_handler(void *p UNUSED) void execute(expr_list *ee) { - int verbose = isatty(2); + int verbose; expr_list *l; expression **exp_arr; int count, n, i; @@ -453,18 +442,22 @@ void execute(expr_list *ee) #if defined(_OPENMP) threads = omp_get_max_threads(); + /* Make sure the number of threads no more that the number of rows in + * rasters */ if ((threads > rows) && (threads > 1)) { threads = rows; omp_set_num_threads(threads); - G_message(_("The number of rows is less than the number of threads. \ - Set the number of threads to be the same as the rows = %d"), - threads); + G_verbose_message( + _("The number of rows is less than the number of threads. \ + Set the number of threads to be the same as the rows = %d."), + threads); } #endif current_row = (int *)G_malloc(sizeof(int) * threads); count = rows * depths; n = 0; + verbose = isatty(2); for (current_depth = 0; current_depth < depths; current_depth++) { #pragma omp parallel for default(shared) schedule(static, 1) private(i) ordered for (int row = 0; row < rows; row++) { @@ -475,6 +468,7 @@ void execute(expr_list *ee) #if defined(_OPENMP) tid = omp_get_thread_num(); #endif + /* calculate through expressions row by row */ current_row[tid] = row; for (i = 0; i < num_exprs; i++) { expression *e = exp_arr[i]; @@ -482,11 +476,9 @@ void execute(expr_list *ee) evaluate(e); #pragma omp ordered { + /* write out values to a file row by row */ if (e->type == expr_type_binding) { fd = e->data.bind.fd; - // printf("binding: fd %d, tid %d, row %d, value %f\n", - // fd, tid, - // row, ((double *)e->buf[tid])[0]); put_map_row(fd, e->buf[tid], e->res_type); } } @@ -540,6 +532,10 @@ void execute(expr_list *ee) free_buf(e); if (e->type == expr_type_function) free_argv(e); + if (e->type == expr_type_map && e->data.map.idx) { + G_free(e->data.map.idx); + e->data.map.idx = NULL; + } } G_free(exp_arr); current_row = NULL; diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 1a54b07b162..6c28d08e2ab 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -185,12 +185,29 @@ int main(int argc, char **argv) pre_exec(); threads = atoi(nprocs->answer); + +/* Determine the number of threads */ #if defined(_OPENMP) + if ((strcmp(argv[0], "r3.mapcalc") == 0) && (threads > 1)) { + threads = 1; + G_verbose_message(_("r3.mapcalc does not support parallel execution.")); + } + else if (seeded) { + threads = 1; + G_verbose_message( + _("Parallel execution is not supported for random seed.")); + } + else { + /* make the maximum num of threads be 4 to avoid I/O overheads from too + * many threads */ + threads = threads > 4 ? 4 : threads; + G_verbose_message(_("The default number of threads is %d."), threads); + } omp_set_num_threads(threads); - // G_message(_("Compute in paraell. Number of threads set to %d."), - // threads); -// #else -// G_message(_("Compute in searial. OpenMP is not available.")); + G_verbose_message(_("Compute in paralell. Number of threads set to %d."), + threads); +#else + G_verbose_message(_("Compute in searial. OpenMP is not available.")); #endif execute(result); diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 2e5c48b5c6d..664c18ec8bf 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -539,26 +539,6 @@ int open_map(const char *name, int mod, int row, int col) break; } - // for (i = 0; i < num_maps; i++) { - // m = &maps[i]; - - // if (strcmp(m->name, name) != 0 || strcmp(m->mapset, mapset) != 0) - // continue; - - // if (row < m->min_row) - // m->min_row = row; - // if (row > m->max_row) - // m->max_row = row; - - // if (use_cats && !m->have_cats) - // init_cats(m); - - // if (use_colors && !m->have_colors) - // init_colors(m); - - // return i; - // } - if (num_maps >= max_maps) { max_maps += 10; maps = G_realloc(maps, max_maps * sizeof(struct map)); From 76ae210f724ccb6d2764219741ee62c67201925f Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Thu, 22 May 2025 14:50:19 -0400 Subject: [PATCH 10/25] Update benchmark_rmapcalc.py --- raster/r.mapcalc/benchmark/benchmark_rmapcalc.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py index 4d9e5cec60e..ddccb135d24 100644 --- a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py +++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py @@ -1,7 +1,7 @@ -"""Benchmarking of r.horizon -point mode, one direction +"""Benchmarking of r.mapcalc +raster (2D) -@author Chung-Yuan Liang, 2024 +@author Chung-Yuan Liang, 2025 """ from grass.exceptions import CalledModuleError @@ -13,7 +13,7 @@ def main(): results = [] metrics = ["time", "speedup", "efficiency"] - mapsizes = [10e6, 20e6] # , 40e6, 80e6] + mapsizes = [10e6, 20e6] # run benchmarks From 224e7f3bf47f6a55db0efc08e4048e55bb9a4b6b Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Thu, 22 May 2025 18:21:18 -0400 Subject: [PATCH 11/25] remove unused variables --- raster/r.mapcalc/map.c | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 664c18ec8bf..2706be16a6a 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -152,10 +152,7 @@ static void cache_release(struct row_cache *cache) static void *cache_get_raw(struct row_cache *cache, int row, int data_type) { struct sub_cache *sub; - void **tmp; - char *vtmp; - int i, j; - int newrow; + int i; if (!cache->sub[data_type]) cache_sub_init(cache, data_type); @@ -181,7 +178,6 @@ static void *cache_get_raw(struct row_cache *cache, int row, int data_type) else { i = (i < 0) ? 0 : cache->nrows - 1; - newrow = row - i; read_row(cache->fd, sub->buf[i], row, data_type); return sub->buf[i]; } From 851d67c9e3d3b9b448647a94075371441953f81c Mon Sep 17 00:00:00 2001 From: cyliang368 Date: Thu, 22 May 2025 18:56:49 -0400 Subject: [PATCH 12/25] remove unuse variable --- raster/r.mapcalc/map.c | 1 - 1 file changed, 1 deletion(-) diff --git a/raster/r.mapcalc/map.c b/raster/r.mapcalc/map.c index 2706be16a6a..bb6475a1f6e 100644 --- a/raster/r.mapcalc/map.c +++ b/raster/r.mapcalc/map.c @@ -497,7 +497,6 @@ int map_type(const char *name, int mod) int open_map(const char *name, int mod, int row, int col) { - int i; const char *mapset; int use_cats = 0; int use_colors = 0; From 96329ebed4db908aa87d0d391f7ad605af99422b Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 6 Jun 2025 13:36:03 -0400 Subject: [PATCH 13/25] add more benchmarking --- .../r.mapcalc/benchmark/benchmark_rmapcalc.py | 43 +++++++++++++++++-- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py index ddccb135d24..5cd747f28a4 100644 --- a/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py +++ b/raster/r.mapcalc/benchmark/benchmark_rmapcalc.py @@ -13,7 +13,7 @@ def main(): results = [] metrics = ["time", "speedup", "efficiency"] - mapsizes = [10e6, 20e6] + mapsizes = [10e6, 50e6, 100e6] # run benchmarks @@ -21,7 +21,6 @@ def main(): benchmark( size=int(mapsize**0.5), step=0, - label=f"r.mapcalc{int(mapsize / 1e6)}M", results=results, ) @@ -30,25 +29,61 @@ def main(): for metric in metrics: bm.nprocs_plot( results, + filename=f"r_mapcalc_{metric}.svg", title=f"r.mapcalc {metric}", metric=metric, ) -def benchmark(size, step, label, results): +def benchmark(size, step, results): map1 = "benchmark_r_mapcalc_map1" map2 = "benchmark_r_mapcalc_map2" output = "benchmark_r_mapcalc" generate_map(rows=size, cols=size, fname=map1) generate_map(rows=size, cols=size, fname=map2) + module = Module( + "r.mapcalc", + expression=f"{output}=({map1} + {map2})", + overwrite=True, + ) + + results.append( + bm.benchmark_nprocs( + module, + label=f"r.mapcalc_simple_{int((size * size) / 1e6)}M", + max_nprocs=12, + repeat=5, + ) + ) + module = Module( "r.mapcalc", expression=f"{output}=({map1} + {map2} - 2*{map2} + {map1}*{map2} - {map1}/{map2})/2", overwrite=True, ) + results.append( + bm.benchmark_nprocs( + module, + label=f"r.mapcalc_complex_{int((size * size) / 1e6)}M", + max_nprocs=12, + repeat=5, + ) + ) - results.append(bm.benchmark_nprocs(module, label=label, max_nprocs=8, repeat=3)) + module = Module( + "r.mapcalc", + expression=f"{output}= if (({map1}[5, 5] + {map2}[-5, -5]) > 0, 1.6, 0.1)", + overwrite=True, + ) + results.append( + bm.benchmark_nprocs( + module, + label=f"r.mapcalc_neighbor_{int((size * size) / 1e6)}M", + max_nprocs=12, + repeat=5, + ) + ) Module( "g.remove", quiet=True, flags="f", type="raster", pattern="benchmark_r_mapcalc*" ) From f960c32ac3eaffadc2c49e00d7bfef68b7837e33 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 13 Jun 2025 14:29:06 -0400 Subject: [PATCH 14/25] add performance to manual --- raster/r.mapcalc/r.mapcalc.md | 20 ++++++++++++++++++ raster/r.mapcalc/r_mapcalc_benchmark_time.png | Bin 0 -> 33257 bytes 2 files changed, 20 insertions(+) create mode 100644 raster/r.mapcalc/r_mapcalc_benchmark_time.png diff --git a/raster/r.mapcalc/r.mapcalc.md b/raster/r.mapcalc/r.mapcalc.md index fd808d8c384..e2766668f56 100644 --- a/raster/r.mapcalc/r.mapcalc.md +++ b/raster/r.mapcalc/r.mapcalc.md @@ -872,6 +872,26 @@ X (map) values supplied and y (newmap) values returned: 100, 50 ``` +### Performance + +r.mapcalc is parallelized using OpenMP. The number of threads can be controlled +with the **nprocs** parameter. Note that more complex expressions can benefit +from more threads. By default (**nprocs=0**), r.mapcalc uses all available threads. +Use the **--verbose** flag to display the number of threads in use. +If you observe reduced performance when using many threads, try lowering ther number. + +Note: r.mapcalc may disable parallelization in certain cases, even when requested: + +- When a mask is active, because the current parallel implementation + does not support it. +- When the rand() function is used, to ensure reproducible results. + +![Benchmark of r.mapcalc](r_mapcalc_benchmark_time.png) +*Figure: Benchmark shows execution time for different number of cells +and different complexity of expressions. +See benchmark script in the source code. +(Intel Core i9-10940X CPU @ 3.30GHz x 28)* + ## KNOWN ISSUES The *result* variable on the left hand side of the equation should not diff --git a/raster/r.mapcalc/r_mapcalc_benchmark_time.png b/raster/r.mapcalc/r_mapcalc_benchmark_time.png new file mode 100644 index 0000000000000000000000000000000000000000..c499d491fda2b979fa91d4822b2fd1465c29c0d0 GIT binary patch literal 33257 zcmeFYcT`l((l#ldLd;dH9Gkf>$uCA)CuCDr3lVG5yK}y6(1c5+EH8HA25C|NA zKw!6#aIod3O}aby(-&c45oBZ^3it*1y108e13}?_&VX~6y9)#ow%~1`Oa^C1c79J( zgd!*&uZJv`DD504B@VUka;6}qQ)143TCYx1q{7Qv~?m@2w`( zZR@v|r3%}{ybuU9#{KG51I?>f|0NG7Mc%Dc1xzoRx!?9$Glu|22+DuZAipM>*1yJ- zTP4^8p{i*;UbupbSgha45!qdB3b^~2W4jbgpWUuk9e z#!?sjM>3J|*b_v7$HP=ju2}2aYRK*)0(%lWn14d{w;JxURhjO?)%AUwd#>S+zX!%- zL=8b6qx0HAvAZr@r%@@W%&mihclnU3@JM*D{&rtP%e+m9GT8nR{9ko5+8 zc5Ts?^9Bj)eTHh&-0vO-TCR5XTh>ll>_7SuTNu^S*%zF131bgULb0RONtr5MsN6G4 zrlw+8(2w{U$MPH_m~G8uw~?hm!^=59{)^ywheVR`10N%9!l37KX9L`*CA^&A#b1go z2U~=wF#IP?EzcJqE{d8@Tbp+%#-Ec~TenJtlRsO6-hC7_gk$DGzqRFueRMt14Y&J5 zyK^o0hUdZh33SlNe3$9{dlWsN;aSc=#|^m~-@oWoKW97a4w!plV%7oLU6?y)UlzJL zvW~vqg7!|n4$gvM-hQCrK_K! zH+M`#fU|Lgo{3|Gr=yG$hms<9tiXH@(Gj;Q{ebRR~CG}XcpoC{&WfQ zRN%1CH2|*q1~>!ag5rY00%~FIA)*|LM1XvNlZ&j8s`}p~z&8aBx1b# zK}ypnFv#A=(fNWDNLfx_jHZIt%&vxc-^rLgBK?2AT>SqJn=h{-?#j%Ra~j9H79VPEPhg_Nw+l&LB}?Q3+XL zNm&sw6A>|4acNm`aRFg*S>b>5_jPi23IBiVeX)4}`M;JN;~ogkAO5H5uN`IV?Elx> zU$0*7f3^|;{Mi(;_Ktsv5NIFb?DU6Ekk((H9Np}FT%AGh@wd7DtKI$oVhYYK&i3~9 z&Qb!B4pPnn;x5k40y1KtHHnKmIf;UD3OoF5P5-7I=<5;`Y9HXN>u+e{=PFfUZ5`& z`nNOuTW27@|9}4dy%zs}P5}V_o#cO|-~X!Xzv}v5Y2bg=`M=inUv>SjH1NOb{9o(( z|4dy(|C3HR`+!kUD3~mbpr2KMiPj4P9TPS1mz;u(o{^4?ot5hfrf&cQa>LfhOuk}0fl#>5x|10J|NANAzC6M0*ph~MOXM+r8(a6VQ9eFj>|7xc5&E%W+6SK4PGw;?` zc6S?h_jis?KK=gvd+~DvAt(dXL0dx=5|77%0bdCHFqVN32rVciZ3WcXu>TX-ND!o{ zt46SkPf8)o@Au~V3kU>&XsRlkge`o?Cw^YsGVuv-D=kftUFhH;=W@)F`Fs#lT4S7Z zVf@G6wb$L%VUa!NUcK|#wz&!|ZS$j#+b}uy0z+n{P7{RJ{2xx}2e%1-BfQQDkZvQ8 zP!;T`=9v6X^G6#e#uNZzkQv}a3?{1OopejD0}hNRe@bAj=$H!NOwxW(=A9p512HHG z&(#=B=BvCf@8q*k;6v*^n3b6xJ!WB@1vjmHW07TwRi$kk)lH-TNHcGIZQ;w{Izkw2 z7!u!9()LWn8`b9`*=vU>N-L}0riTXDqUrf;u8=ecn(2wrYW-X><#0%i;h?VA9KG*Y zS+o>0#gmja-svOub;*%13$M#t7E)qD#5!V=bNE3MxG%)@G3q31Gi(8q{)leqjaZ=l zsIlA?Cv_z5J8aFh=*xBKW4NFX+r_F)|=vMTmtsvbbd>m1ct&%!g`zn6>u2>8PXedTn(( zdxd}Y;niD9&uZqBLL%d4u$;YAV=xE3>H~ceEYV}ru&zgM$89d54+$|62%t|1_Os$i zw&$Jeg)Q2|lLzm(gHtZwKRCsTFHke%_2rq_(h(RW9wLUQ6G8PKK-}+QceZcV=*BET zY>Gs+A8;mEVsjZf)_9X1q%c2DTyLr4gBdZFq(IZipYp{K4_g=&m)v{%*K~3SP6=$f zB%tCU7uYL|PCk_MymYE6<&n$*r^oMUCB=xiJ6(qG#94CyeyRf0F(_ak&gr7b6_Lye zO(TJ5TAIL2v|(u*>b&iqh`|Mjd)B2Ok~NmV`k2>23|vxz4@2<)hIvw$(-<}=`qmO; zjhIU3s;Fe1w39uYLwz37SbKvKh`j})gl+u7raRw)pxmX=wrX_ILGIT}pF6moTU_p@ zO(l(Q68-#ZH+J>W&bx4MG2B-XGj12&{n=dNUDtjI5`nWLPDsDiV7Bo zyLi2LA6kYEg-{RfroycmZ)l3N9S>Do|JY5l(?Z=4V+W|ON!~f;IW8bT3z4is+L8k2pUm& zlnC{=J`{fYbU8R+`q$f9PyJvho*MbI?dZ_>4MW|K($@Q=w70%uDnF6^3LqxNAs7?L zOifQoGAOg?OI(zRe3A0LctotP_7t8TLV3D=F-DLR%RO#kU>36HY2nTIy-W)f72&6P z;`jNd74eD#qng9K4ux#A!Nb=5%IUPa`E4~2H1kuGGh}R@DYM`C65E;9a$`d! zGUc5QF|Mq{>!&yL{QPnG$T2HWU|I`R-F@P>n$7R6@f62=PbZO9Sz4>((hdE5&0%mK zo$g0d#Ns3Dco|rSO?9g_lSriAyWB^<?>X z@%*hb47Y%9}Ta5p(fu;?M#T`|(;A1om2f@X#VnWq2;C&P16_`?`m;SJD z{b2fevYFGpiL3XdIp7%ac^^Z3x^A-1x+Fz5Pk(im-hOI=fdP}dM9NF6bDHZSS3X(z zDl^hbdc@U2jSNr$#^|c;x25ysliuh&1$19Ztoy4re9mcqfBxXv{d2!Oj8FO*Hs)>H zicbI6jvvB6;n@HQfETFKwp6%yIn9o)gg%KIyIK{*eb;_L@OCm*wPBk<)(ozoKNbJkQi|l(? z4xO;PBB7e8a*t%)ddI+2k38P5>Av&3Z)H9mALgKF^iZj%IX#+E;nO79>Gho{%kdB_ zvO=pt<}Wg&PElo$%93uE%55g%wLW=&9gZ38262*%E%2ewFCB9Q4bCl=mF}O8af=y*Do(4-I-XwT}y?3 z5syx2uue{kOHOiM5f5%>_l8DqK3PDnBoh-r`~m%pUg;ZlN=%-k!9A*943ZyAkyUlx z-%e&BnG6ipiVr!#54%iDBIX#&Xd8Ocbl$fM{hS|RjHUiF4^@d^@w!eS0uaUdzIGnx7+koN-Moe7y`hu$}$&NKZsbEn+ zAbyUkIOP!`JS|3j4G)sgUR)}Kz+Ki%ON_#0LNS4ne5g-#2=aq?vHLK~6*hpBJrrGfT%`7Yug4DYK_#66}6 zO_Pz~U+G7FsYVgVDN)b@x@#XhE%NjCSGG^1%GNB|`steYX zsu6wum$uF$1Rm`7jeaS8{zydV1&3l)h$i=Mvhm*g%t9v&w86Mv-E7|s%pzrO(+gLi zN}!brrour|lPM+dxE%SBEc`4byc|_jmt)uQ1!>YgR@Gq-PvG{R1`E>zVqRubl;0PU44T>Fc zC_tH5BO&oYYoK8 z#v+&X0%qFo9KTc`Ohf8DgkH37wVt(OPJH)~a$QD?+m5m5AYt74y?x!_!zDXch{!bK zC26hUD3^gmf0DF<9(yiW8s=)TEJqC!uH}WQ(C!H%ltL0%Sl`}JDg7EYvT+dY;ah1o zf$%t9*n8`@%Vkljv0gd7A0lIGyH8qICYQ=tu_D8KdAyhBktjX{Q<5PsLH+E;m7&sH z^9-X#QlRlH;S$qND)TR8* z9<8;dqBItbbs9BSkV(1*Pb=ZZdlbYRiRBqdYw~zu0tFa|$(kH3NCDx^{3skLtY}fk zlo7S=PztBxQ{@_rKsU7pB_^4^hKh#N-_^PctozUR|R7HxiX$)23@`+H5r@6=*Z-9O~LT%oDVccV$C-WRidm;wk~AEDRqWlcw( z^-0`XYDHkQl?Sz;nt&NfL_U{_CA!xur6qoP-I6mhM>Vaaq9=Icd%fAP=vfzF_Qke~ zi6*!vPmtC8s}6~R(x|%LRfcSv_iu#XUWvK1C4x}TSkFQ80dqDAsZYO4SGJeV7G12 zdx$y`|2S3S1#2Xd`Zwqc-@=_d17v~N+=sy#!QLT(D1UaK?(iY#KQ}fZMaw)Xk(Mv} zOhu4SI!ZwY8ZLsR*C)O*W=)F0GhX3AvZJpv0}1#}qRAnmf4CH)#EsUR9WfA#`iuP4 zK!Kxlb#Bqj`z${G_62$pqZ!5vw0|XEHim6*K_?O2g89(j^>Icr{D+D9_$6f*lelBW z3+2Y|thTSI&=Hq+8iNN)CxCu+vtsDiSU=W}NJSA3o%kVYWNrS@Vd;$One5ms5aS=Wjv`)pUzlj+KpBY&elFx_6UcqVJ1et@Vg7<-QbyO!VL32F2C;EG39n za<~5jY<5T+EzhrQDB^364iHP4iyI+b6A!WA5F{um1)TD+CN1x(ZNUYE z{;XbMQ2>a(Qtac-l9Y;Qp|8BSwSGMv-SKlBot zUP%!w`zS`vg7hXHsLUF2)IjY`mq%}&7*}sH6{4yl$_c8{KD=jlAC~5vboCdW6>{Ki zuVDU;Bx1aqjE5oH7yr@2b(vqg_UEO!q7RO^%O)f08;>kv1w%JSkV#^y24W&DA$@Zk zcFu$N5_dTis8@S}rAbDJTIGan>2Xx>A~A?N=$Y7nlJlK#rr>b;4H_BTN*%*b zI@tPK+UD01K}j8y=QClgx%A0I3uIZhSqN@g*4fqKCR4!)Nz+xIRu5ejE|dNES_Hep z@v#}d@9kYQ`X^r@JmHC=BRqCDx~s`B;pZ)fAm(ZHCqn*&i=+@EBUFPT`Xa{u-IgN# zo`Ma(V_C0NDPHczq(y8wA8MUj@P3V(=LfQz2C5A-88*`c;!p7$DlX{pJ-k%k^z6crewmfNY48_C zQ6Kr12TW=Tz}pOxU@k9d_lapM0|W%XH=`<9#HudHBC3*276X~*BPtMVDF0=wN zRw(;xxJowV+PkDdwAvN+nBY26203~W=oK?r z{M8T=QHaxQDo@W~Fty+o#PFW50vQL;X5aXqGE3lA4*o z@DiWp&f#Tmwy4K zj>}fTp8PrV%9(BwmtPB@7pU=P^TW}jzx0-%u&W!l358xMN01b67syDEOm=kU<8$_n zP#)m_Wq~zFTtm~Z386h$dR5PN6K5W}uc<0!n74;1BhMQre`K`6s${3XozH#}-|W>W zR4kIwzviNOSFDluXNo~GewyNo#!&DOQ0bT6Ia$)t)d42EH(&e1=hTqoI)8C+l zYf5}@@Kn(c{_>%wpG|ut8aknT5lau6JmgDN)?LdhHxeQSy`iE6I!ua3UanoBk3sgC zw0?ERgOah9i_TmTQiK?p`jyqJyy;&BY0{f?@(PW&Efn8L15Jk^z1W*XM{a1pVL32N z29ZnU3-mIN7>t`KK^F_gl|;t_Rp~0B^s9= zg`*?JFy@J}agpg6VlT(v|A`E)Wozbm^Gwmjcy?n57%Qy%u$sP_}r!ztAg zu4m;w-MVzO0F$Z+=sS z<&EZaI2_V9Pg8rdkazSFP5J+BO&bI=WbH3pG>(12du*HSPl$nU_YFwO@*ygwc;5#5kn}-i0YK29f@x# zi33M9w_|ttB3;e>VkpyeVMd0i21fLUvxfGcJq3Q&3@3cU#W0Qtk-32tBFA4U3nrOd@34DX;yMO#^9aJ+oM$u?lYi?w(eSe7}2+= zsml6b5w+XMPsDNkrD^2PMEbAsxrmr39F2258N4XlP5pp;gKEK}pxD1o>;9`oyr?R~ zo)sIUZ`M%109tR@#CgbM%#Ko%p{cOBl0;Y~u6o0>b8IZtZuZF22BC{HS_OuS4JGdT& z&`1Vyf@7UNLkp4kIeXM33G7^D-IN($m8CQ!-Qg0z)0>^_P_>n@)f9PR6=S`?6Xgvn znq5vydy2rF?;RQ$_Hu^k+`*x{6Pxs@Dz9Vf>3DDSaGZy@nT856({gI=%qQ z5jK$CbdZUZzqIpAk0j2oY@Bc#I%J886+F(RgMsShgK|P>XcS3O;wqjdkZZmCspDFo zXYutbx4voj@T~rk&gK)lOQ(l|wKsQY4x(R~r4|R6wJYC54`RP7!sJx78p)(PY6 zgepe=aI|g0$e~?Y5Ir6|IT*VLU6Ht6W#p;DCGdfP=OT1R;^4{{t-r>SsE87&fU&%; zfREp(J_^N17@m#gx6GZzguSjR=M^#4N|GiC-jzLfDO1WM-}0csnq2s7DbT3Yj7>B6 z2kb49lEljRuoBU{Z(o{VGN9`nXospYAORCo-njX|RJeRz=JoInbbta6DtTbtfa zOyP-|5v`U@OonzO)vMGv7V?cn{`jj0GXqAJ>*J=^>hqTdm0zQ9lg_u9@h%<2f3u;i z{86vyGsvNlFTBPde{-wy3ygM{^k!qE(c)5MWAEn)`96H_k)E}7qSWzkHY)wJgcw4a zdz_bZp{iiKg#czVSf{6JQ7GU*q;&kH}s>Jye2Rs)Lck z)(Qudbk1Ey#X1Hbm=v6g{4mcT@L)F1NbkLKik|>geP%%O#Lr~)>q079pwrOpt*S}% zLX@;yh_tx9u_@cF={G&6iw74m8YTAT&&}PVCi3`F{))(>-Mo9z3F_v4y;0Lp+AbFJ z)z!MkPip9Otm|j#9H+*pJls@@n96n=d>J~4>PC55QebGERpU-ULXAyn$4$b#Cn@=N z^;^$Pu-yAaS}%R^!P9;Q5?=?UsIYAR31y5lsm;MFP7SefifazvbW}(#=9V|gZLJRx zsfG9K#SG>=`Es`l0~=$bzD@jD$+DJH;)m=Fu8n%fX#A3NMYe7B7LYX3{kjT6mKOm? z`ypgw;7-2acPz;Q+Xr*Am>$HNFXU#SbfUt33hx>h*>ra7+=E|;*xj`t?rw+|$4*5= zTl>Bb?XQ9wuYSk~t;jDpfVv;&@9jT~tVYH*U&ZHQ2Lrm=4jX0MZF*c(_;}B*-5ngY znErX{FCe2o{ct)Ei8Ja5j9mfjcxI|!pU!xFupgqof?FHD^XWwG72bjXM77iy>HXv(P*EjXn+o#&iAFBElqm#8t*%D#Ms0@@%cG_eGAqkDdwW>hy= zYT|l6GfD#c$YLr6AC*wL=tG3E2n7*{b!g?jB|2!|3V~2By|$yydS-t&B<=?*4i@2; z!u7#aUU?=!`cNF0drw=V7oWh*_lZt!H$u&;A$=Un8TnI#0>?kxf`_@PJje^letP(t z{Ww>q?58rf&I9Q54cbn6`V^e+fZ2C=<)kv`412g}YZgl*ksW0U>PnOd0w#$ac!?~(F*Pqnl z-|rP>CC@nNzt#D#hB%5U(b@9RX%SPpaq2H_$DEIisv4fd6poCz zmcDf^kN~{MLJ4+q#D;v%p+RT+eNWc)@TVlGDzp|BgtZKh#5b3n%`9_Fwc#rLcdp9< z!P7(DkX?O)cXT)+4MCQJ4redA!?U;@;ldR zg4P6i-Hk7JCGkgPG_m&))ch2pw2Qd~Ap!swf<1j_83E9#WrtsBRaEx3MaQs1Mm>Y%g21&0jnOCgi)Nw6()-TPLY|L^!4_zxotbXofZ5l%fp%IH@2{&^ zFD4w~t|=*fxDt{F{k8coPtAmnFi&XCest4v;q92g69gg=%Cg0q7bW+#=0Pe<&L#B$(>z%+y7E9C+nlBfjmdDFY9LZKo-@i(B{XG2wAO`MYUO2 zTCR+pMKzb>^$h_f%lgtaB}9CXL8mh`;od?OQMu-ekb>(^ENv@Xr18zK;U+j6b9F*5 z3If>=t0AlWDrnxakwSGO(uB>__xmh7py*DJ!UYoL{^JtFa8u`~rq1>SynggqNQBxW zXNv;!h+`GSD~*Ru5^tot-8R(wQK*UM(;D!i+QwfZ*Mc>koQebnTZI)uD6H9(p&a@% zl~&o#-|A~~8foT{nydJDu%e`LB*E4~JHLbg#iRgKW#5O5Aby>;YkgI=lvjO6s!ARf zh?a5GXq$h+qs|zdLArSP!!Yw&5&825$dI*#uK2UQ7-#?IcSvuvg|kJQu=1miOz{W_ zF=cuzNCe`t&wECSXRs#`vI+* zq}RJt;WxU+B*ZRKaykd!8kiA;rd&8t$Njxy(Glw4S9bZ&Qhw)2PL`4r!r`xbwVkvUej{*i9 zUm4Lkz%jB6rdQfOh%=A>eVcJ%_|!ZkWJ+{WdS8168+SwwS#|XYgHoPQT{Wo^Zq$ zMvAFh)B<|G{w>KpbbvlJcDdakuVw`03N-k%C&scOq&FTS?}Qy6RU?yiI%Bh-{l>?|sVXU}v_1u3)#=IX9B1&PN=BA? zrTE`)&if>z%I1k{Zud|h=0)D{?6+@?4vf;1GY@RnCn#?OYjiLATp*pSS^cskn6l#% zMw~lw(N%k&Yq%36Fb#x3(DF9B1KAL#m8hv3k`wW#)3E%0Bfae6?*5hw>=Uy+OAxQO zu*Hbt!=FE`+AVguq5uSnV8}81bFoh$D4QQ$PE^=u-yq*b2ES@xUic@K5X`>$m{!;= zjCI1f&(U9?+T?huoI@mltd}rL*?TH4{;U(yM|H7Q>nB+86k}Fs$w-t0K5F@aSRM?~ zXz^sS#gNwXn+CTRoD24Y781{0aIn(?=~+{vj*%z2$P;22-@K13hko9u_KKV6coW6H z?T=W$WSJUnMl5-r7m_Bg|1*FVwc7K`eY`aMTG)phN7COx!!`0FY2Hk@DxtG(as@(p z*ZtuK@e!YQ$F7dp7Z-Sr84^R>JAb=$k*!3SGNoG^UcA@dSj)K8VpXRZ2x)t2qd*)# ze$+bl(c#G4=WTMS`*}(x_6iLUUDK9UUk!-7p*SkAYP5kbk7i$ z_no*6j!BS7`SgLsWYANR0Iox~syfQIUwhk{m=Pgq|6us|oGDHQlp^Jq(U+r9d;MoH zsHWqkp$YH$shm5vDfWuwG_9ZO`~`;XQ5aImtuk<1O@52Ex*d+c3!2vfQi2MZjqtZagGaQ|^ z%u`DQ))4N$j0g>L@O`lT)BZ?&?8^12@a<>Y;ewJRV1FSP?roy^8-XB3=*ySzD&{KL z$}Oj6aw-<4?|69}J!xuwO;lE|pV0pm29-mR`u)vVzb%AMsOV!gQ-L8) zg;T4>nmmw3ikbBC6GbDZ>*9|u6T-V*cGPaZV;8rsp5#V?w_M`oz37uIcpRJK76gq6 zDEwz!>9GqQGE% zG_$!ZtIsSLL<-T51P;vaE|5S2<@Kd77|Rc`-PPJ6esOI_u|_JdrWT*&5kV;;gMN0U zzIig#_I~X@p|}l3Kuba>toh~FIcC0ZLGOg35(F8YA^d*PXah|jrSF( zW^m+6#2@>EpwU#PzK1HUDJBv;xGh+)pn&jB|0DeI!b|uxLRrNt)i07GA4*xgW`q)E z?*W3M-zo7P-h3KPjy$1GuqY7Zh_-9*doLSpui~L& z3~xV92uFt~P_YQenHp=`5M|OaW1h8=Vbqidm!b7J^g6UDBMLQ$Qo7p@!qiYRD z7knkhWB)=rNe(xX)3ch<6@y~lwPiQjm`UjRz+aX>##VLNMbP66sym@2c5=oICE4)4 z7#7&v3)$3fx)4qm1EsKitO4HSG_)wmeD;EL2yKIZ(GSO(iIkmP{zVlwADU3tYgO@R z;hS3S-2%>VD8S~ESprq_Osg@`1hat zZ9jwV_!YgP$VCN9wY=twK9M!jh>qSNoyuMU7%HtQ@CpsakZ`eg3LkAh)xDwXv+wf1 zxXU00x?6>+j?B2ND=W%oWJsFd7ncS9Mpf-@il7qsafvYOK1mW$*C!i+Qf3-%MwqA-$-5dD zN*M%&0@#7u>33NGxzg_$NmhrImEXSMxyp4J4u`}+G^^@8Bt0Qu!15i&;Zgs)=+GDK zcxPp1#`3$q>X#6L_=@tw!r(fWy-r;Pe|r|-?Eo{T72ifm`32+0n_oFj>OSOuRF8vL zz8~m<{;4+OzC+f&>&{K2P{0*=sq!a;;@_sA6xB)XrbVV+QtN=tbgQaH0ce738mOP?9<~6X6+ac^Upzy z#%g`!w=-mvi$@PGNhibI?w|&f<-t3GLoBpqnj7ht%-yeqe!gn(+$L>!edbfhPTx_C z{i`C*nD#^O7ci;RUzBIh(L2oAYPE|ib8Gsoz==qnrIK3$WmDv)I|I= zgrfqsD}$Zt-)-0~kRJVb$1M|YJn)5~2Vb3qy|=qfB84{@Z{D6luFs{4!vrjHnMs%A z(1Bd-x}nXF$v{5BZ7Bgk&C^?F+G zLh`0U`NX|=KWdz;@oZ|xLPISS_b>_<2w5H4H9CYloR&;~_Oc_=^5?rng^jwAgHSa6 zh}&Zu|cy6 z*+7fh+5OIc7h|2E1`H*I#!vA;(>p{C7W|Cd5s2iOzzrCp5Lor9+@Gbz3!O)dYr@DM zmzu+wst2ExKi(Ptes%V2y2|E-eULHgeDJ9_KzWINR8Pm=@YT^-0iOFAl#m1%j^O4C zpvG}R-T2&_zFI*ld5|)2T-0O1)7vBtU0EJ?fayN+m!Ir{Ls!KH$a@P2IP->2ACA6* zi}ls=C1K$4o2@Ve}T^SqbbdvnoLmrRIAAqSt3k{;{Ad}`*zw3KVXZy zT@ofSFtZNE0u7ukG8l^(8&UEv1eN?{ORqPFwjn&Hm!6L^4aVI`o-DZf>T;6RHRXtv z$Gkh&rx(1@(k_mA(j`2y0b<;>5+RSf9R%vcm)6-a`V6;V2}E zqF!s*fwo#*?I{7z41{KBA6h9)AdyygU&n|FcjJH=fcm`8B@K{UHLFfH6Ww$Q#%{RP zIoO#>>8kos;N(}!Aw~_-$EIhxVbuxB0CJT&W{LD7#6IG&Y0`CAV3WQ*J;UzicTiX2 z&h&`8nb_mJ#nCkMe)i+mqI(2V^)GE2sE=MGxV2SZThoIU;j;$bh{oZyx$W~~@#V9K zS;GV5b+6KvRkBtPda^%ltGAtE0q~{{?axdaMb@?Vk}OkPKwA zj&X;eMB-2-Bdwv?c*q!}*VGl8;6%xvL|RL7-mcKq6xqz?$Z@Li(lpOSL%HHq%(4MR ztv=koE4n)m`HgNH>hKD=pG18miG6>8PWX@`w`Y{Tsr9U8KApI`t(iCTHFe6aXMHY= zaqlC-*2zDQu!qKi>u}l3y(Wa&rg4GmnGN(`vDJbrMxYhH+;OWFx-m(W z2OFHgGaf4j6KSVp?8d~oEPc@WyN+YKdwMtJ&)^?cjg{abP*u=^D0gUHDhXU%kmu zzkO&Cv&%*oLkV$t$&1GC$Z&}l{id-~;jb;@x!xs}ROE;{>IF+(WFxx}D41>ph(&Ol zd?VWam3goX0W%623CBf7z%nROopa(*kxG3`CEd1`_w4VA`!sMFII1vna5hXT@#fx~ z6hak`-C587+bc+0MOgF2U7v%oj~!3j@|Z6=y=0 z&4OqXRGN0LT1@9n^x0+rTJ-uX3ou8*-kWIs6}FGAkgd*26R`9Hk>T26f)_Wt1m0z_RTQpV*10i>9W%g?T+ zmT^;IKa9l32U=9VN7|VdQ(8E$rkK<9Ymh$LPB(7K+?|~e$^5b!EY{vUlv-zjzuFx7 zDT97yGvxV6g=Q&W`u5AuNmsF|9*z*s1ao|_#!sY%L;B%tgyj?kcZ00|vJ<@)Vl?@P_c5Gvw~A1stSqBy?s#M4x(ugX?z!pH6r3>&GQW00s$wkJL5t`WwO8tw$Dal zwN8j9u5_CVs4J2W%8AC7M26Hex_XNsrv%>($-~l1^|vNF)r2U*lx38=F^$?(qDyOzgB3gz}NPvAM7M@uGt3*F|XQ< z*?gjY4KAVhSHb=Wi;3!*AV8YWujRZRC-EtBZwWj5FF zX)WsXy(7$jeVMlEUCtAvlx?R~KOpZGy$87HH4ZBBAMa zM!&cCG4f5oQIg;J^Ij~sG1r~MK1WX`f>d<@eG+!Ov73+K8+Tme=PvH^pYib}vJk($ z7f6b!GIg2zRS=ZNy(zCS`KzE&bFTz1@TLGhN`rj#oC*sdI`K_ny_FAy!UfSmD~$2ifTR4Qw(K{<*rqn0R)S9sY^D1pWR#Vrjzdy|%ist)31!d99_ zEN6stsQ@W}cXo7A4uOU+VtAb9@Yz|#i*le&(rK*oLgOf9o4o^A5psx9>uiPb;#kc% zQEr$zPR8?JdNn=T(`Io--1_Dv3$)}tk5;`RdZDZdS!gn| zhW#L~lO098UiiGUBiA)0!lOFUgLbW;QT1~0D18VYC4GoFQ%{HgDg>BwD8Txk*?rQ{ zK7I79gm7a3cKYbPn@jh8zXQQmU&3iWH##qCB$-7=Gr=NXx-b3Y_Mt>z6Sm6S=muNQyrZ`HS4hptp(}2efXgl$Uoos4nLzGSG``|V(7inleV+O6Xf!{u z-}>I6ON--=7q-E7sx*V9Kc~O~UIuQs@e9k1QWCr{r9MO+fW@dR+ z9vL@;R%1B}Q`bY-V+%Ote|v5&YrgbwM}PGez#N}u{%#MMuF*Ry1IbkOO2Ph|3t)Tx<7mWFH(tv*dlh9>@(Ybt#GBVYmh>MD#C^lZHrd|c z(H(5o|Cy84g%~V`4CUpPg=Y|+f8N_Lwh4+cXvl=LT>)g0!VII9oZ36S8s^}(w<1&$ zeWYWIED?i0Aqlur*F3xjiN%;zlBNi8i$5zKX#4mQVlL?pX3o;hmpXln zjAI!heI)BLf`|-PWp) zO$+pfuBY9AOn$0n+bxH24Bz3QYEWI%q%>kwaa-XSB8Boi74irUu;CYG_ocHWuZpS+L8clKpr)9!0)w=N(cX2OAZzB9ZVxPyNtvvdpJ8*{LOI8266K3Sv4tU6;OgqNcWT+m zY^XfHhIG+t6Eo$D&u}DYE2Rz|8N+)QRu29T5pZzI8+`df4cEOsUBlZ;|G?G80338d z-906&zBZZUJwaBIu#5YqI(^6P$Rs8GKXrX|RGUrkZg4AJ+=9EiQ>?*4afjmWR-6LG zDZ$-c3dIT(EgrOZaVQR@&;Z55O~2o{=iGblxqs!AeRn6%%xq@!&OXCj;6xF8^q+(- zbdc$&*kMd|4us|qIzde&*x@x$d9a4r=z5TUAXW&AjtMQ_=mOUeFM|D9v(#K>taP46>Ie z#i`cgm|%!P_D$75nR}a3k2Z*=r!oCM_>LgGB8Ii@m6@ZR*oj3QaHqOZ6Y$Yn@KRLi zAdCSG>EC!1n#-x8f;a!rTtin(=QHkbXf4(@+VDVFE$GV&1b)rcyr?x?(*FVKiejZc zmdb7<%vq=TSGyS8Y$>-&nP9_B$%Tn!s`MY?2uMFfe&AWD!(<KGRwNr2wM0B5F`hDaW9SsxJ!lYUUC zn!>W?a-7B}&&DW@z+%N2J%)@64AR>ug51z!nFg$zX1{)i#!uqrdLz|@J%o^kS9ZLr zCDSqDqfR|2AS`*%K$wBjiVzKuePYpsOPUz+{B$V$oKqxWXlXrDF?zK?WD5GrV&s^J z*=fWyF$$ys4&d%=eJ4CbVqgOHP36X-zs4aeDnzLIP3!Wl(Z6xAK~RBPQL0Oig<)-$ zy&EH9B$mC=id`SWw`?^4##GXOZE!&ilEDYN;G)sUt1Ucnqjmv$$arn-I;>?Z7P5{< zfJcvr1)K438<@n=$L5`}#=`6i>Q9XFARmY|5Dn2oNfAKiA*D81Z9BC=_K?M2jE%bY z-km_wm{k!z<5NqfZR8*pD4Rgk%jDk*MBb}TxW3PG56y}S&h{=nu@~McH-4w8k5j`H zE5jkgK%pT`-^mQj8UfBlb<3nf<%wZ`&~fr}RT+^tN@r&n80Ec)4Mi2|`vJ z^|9slYDEU8$k=`R6hit=B;WwnHTb8goa(X76BAZ|r~)7tQguR#MyO;)N*;gppSb%- zjyila%xbGGNO!mQdCPq?V>XLo_MUfQA>oLAv<9+lu;j*tGz=k?HG$2d-#%uoU&eF% zS@_?%05gRN%;S&8pVOE<2B4o*uR)4la^*KIY&&Yr0AT0zZmEIZD6Fx=Z3azI;3FN9ti_se5}{@KS7NKhAAb~dR3bcl;ug|(lpImEu<@ zt5bDBA|`YfsVSdJCP_}GTK7Cj?IhPSWqPJGiWkZO(25EU#sx1~2%76=3`V7`l?QSw z@*S>ez$&y9QiTDfg9K+baZud115n&b8x4jz4Hz?GD0Dd=a<6_UGBNC$#tD3CNhr;h z`p-*HfaK8*taUj;^lx_dXs1f9KF9v~XE#!`1k%N!{m?mY(z>9l>S0U5UHgpm3uZ(@ zofFdhD^@iJVy|7X9tB0_*2E&EvS=oJ48;%htI(fmV0~HRcVsab-n}FLl#NI~wjuZW z{?82ly7iQJ7?xlD*FPyi#%YN67*v2|&rKtr|Iqodz7&B83?#R}@yEggyi|@bzkCm` zEW{v8Gi`)W?A$BEWe;;K-DiCTz(ZRoL=w zS!_dD@!_Al(k`MKN7Ii$bH@%u(7z=#@g^B})|Rb9+zO4aks zg~XPB9ANN}o^GRt2{}c2Z+ecuZF45BJ&;H`E*B!K<6O~F0ZakxwLp%zvYMG~R1vKJ zj`(5-U^ewRPgU2%O&CEy6c_8ciU?>r0gw!Z;cO^B%{B+pb{%n&&n@hHO2A}mIY1=9 zj5_wga;!e>10sELp{-5N_cZ^L(+d@?6n1(!QWQn2i-WwaZ%|77YBge589#mYlyO$X z^>pT*zgKpYIs#W=Ajfh7n3OIrwsI_pxOrNHwx>lnN>L-7RYyLY3jpFAuUZxalV(GB zyI#ynBpez~XGWCYe0j-|nwI+Sd0d1g23$_vAUw;OCErllCdtWG(!?6XGKrOgV2qgL z8$wJ^Mzz8X&GZ)VQ>X|-9nXm! z%h+yyJ*`U|G+mtiA84;D0t-Exo$g$o2nmj#>bUd1q59`#i;E~qyWA>2`dAWW=P*a z!bjnSfD_!_z7VwPcvUmXG(2EeZ^xOE?weclF}mt}Z4H zhWzT3qfz3~<%HVBuh|E`itS5DRy49lV{#_YjEWOvvGMVIkb zuDjqNek#ydIXUgzb3QH=am${57(SVZc}lx9dci`&S|AjoI0;?d&gP@8=xN@*MQbF0 zd=UXJy^maqjYUlKdA{+&IhS8qyvNhDr*T+txCVZIO!}OVxV9yzFql%*CE}VC9(#f- zVe}ZjY(r$6S)})-J~S1syYJ&6P^cW-pq{w-V|zWHir#5{+*a3k#H)6i8|kyMoGy6; zSU>cr?kGzC&Nl`T(j8hNe}kr!+_{OUvE~rEB@7Qc{IDy^qi|84UhDMGt(?)4F#c?0 z-eT&V5#D+|x(n7X`re$m*AT-0`-I5s!lnT8Cm94uQyL^iH_`fFmwcfR3HeXvVA^zNvb!0t#muyn- z{p67=ev#~FMIMSDrRd>!6iu`GYz}n=qe!N9_#*P0q!Fu_diM8|CoZoRx2dGhEMTy@ zRG~08#-1E(0+4Yo>A&p-josB>ywS%2ajMqe7Lru`DJL>6C{l7;D(mkIerJ~N-prLv ziX|q323iYV053U?jyRZQWe2f{o+{jkkaTlNxQ0NK4@go-{(VdM| zxe3?~u@ESoLujg$zu#el3tB->n3{zPT1!8&v3-1@Xw%w|wI`!h4~CT@Ii9la+tgTU zhdNyo<%gC&)!Pd&Fx{XIsYi)Z5u{jT>dLcn%Hj}ahr-IAxVSRG5Xd?jB&tN{5P3j0 z%vVcWr8uqPbw+#|Vajaap*2JuR* z1I}t_X?6@!La>^qHQM!3h4EUKcov)Pz-%%O(4!K zT2Xp^SVbqP48Oiko}nbF9aiz8l9L}9=pon>15V}ThvB&`-Mt5;*`E`nZ1vvnAFCtS zw%RRjAdk0yb6#LUKLO(AxUP*BPW36-v67Ksi}n`#;vGhh!J>6MUYqo=`>9TS)Pn$- zr`zkS%yuONF31G$9uuN$HzW=H*d#KC?|-$-^uCLXGEF^q<-kdC4-JTP=ATGw?lr_v zeg*S8ay_Ob8^VNIg*#+;`~8$6JKZHqC2&-1R@>mB@FELH9y_)0X!8mduNsRtVGq0 z-{~2~247CBT@p~pSp}HbfIzgrtZMY9J4G0>&$)T>>^P;lsR${W>InN!s24>bQYO(+ z-JdjdKQ{VH;Xo*(2dE@LvV;Vnj$T9n{x{OIh38S#_xhW>qVlisu>hXb_cnX3gp5Tm zsooH{%>Knd!$}ibX)?^m#l&+dBivc;FjDP`iS35(@7~4!fI!->YnqUzbUquOmxc=& z{SxGb-pGhKO2-N?$QBd6lIO7?T4RI_ioR7YT(DZZL~+H1t~IqQ&3`ffvclHoH*hy< zf_j{V6rqf2?OSl((+wLVPg(AhGK&H&t6MoiqFz8JtOtFj?{2TICU_ey!Y-%g0|8Bm zV`L7B&72&`AVqjHp%&r82OLbjYwK&>J5M2y_IhpR3uw=5wPed|O38kfFp?Ejit`2+ zPRjG6P4BOHBcrXRcIF~pch<^s2zvmb-R|*m64P#8O1C)~&_R(;Kg{=2o-HHn!TBo=H6&KnXoy1z&j?l$&vXy6#?^S*hF4c)wS>+UFE$bTgB!&?C!X) ze6f#Ky?1%Acq^9Mv4Ah#pHkQZ;?jRCDlX4EUl_OidBI|xnpnkZ*5;XE=Jr~7Z_yWR zlS^^)ajw&NQ8+mWak*QX~(e`6nqlh{L^4#uUp_e2nsxq zGO=G1#1F-sGoGK+};YEtA{LXNzjQT)GcGix+($*b_|wgGlnP>VeLH> zoqlWzT}|#)>=aBLf3e5#sA3l9Y0?|R_9!b_gx&kzoiX8za(veO9gLX+k9e~yD^t=P zP!<7DGrxgN^=2odca13=l}Xny8f0iBOYKNexEIztiyCZFjrtXY8!GR4Q|tFNqn0LN zhRrMeqH7GBzB&HH&&;<^8=d!tOhWiw4{KlIev)Kf{`eT8kP+%$_{xo?9<~N$%`6iE2u_u(g z@IOs)CSg98dbc#6a%v=eZ6>(8J-9B=lRxc`()6DE)S27>mciLHg!#M=rXIJsJdk#5 z$ow4z!*u3lU*X|i#q_4LI@>2&J&T_t3)A1W}Vc0 z<7F;85owNjd^LfgNc~;MQSsjJX;ENGluDAc^SYY0p)Y{;W4D@W>!SnIgmDIAm65E) zZG}l2$B-rUHOgGNw&vkNt003MsCQDe0x#tme*K%)SrU?3Rl+W=1V@8ild~$#&R<^c zU5upMqF;!+k|){-hi^+&Qg90^@42hJsJmHzy54#w2Zhd@H+?~g+2lp+8hjN^jWsnt zdN}K_v&Uy?S&CrdXCy@(Sv??6QI7I&y05E{z)|J2BH)LtTauFC^{AnI7uPQcxTTyL z{!~oy;*0<1m!zD+P0i+VwW@*wHKil5N0!$g+)Aru4HH4Cq0b&gHPNLj6Ff;pw==ch z%Lr=@W%1sB1u*4NVrg`as={H06w}I*7DKfF@uAr2@3CAbb~~5dw-c3neSLkLzx%d> z^NR`w?SFAN1`CPU>!3`M+Z4t(OJ4Pw5gvx7y~d00$2VGze%9H)yMsO!9vlW6Qvjn= zu+UsNVi%vkROc66JO@tJ zM=4H*IWCbYCHy`~B-ln@-^1=3W1AmC=t*N|>-9)XcU^rV)aj<@9w&d+N`flIpi>X? zcZZKgi&jCBd({d@6CchXBKS`i)9Q~cqe>e>REAn?@K1X_%HLzx=cXjENEX(R`#l+y!$xB%xG_qG^D$x*3vv!mf^9sYB>vFR@;HII zZq&b8mxcJTbE$X~aDG7qp7Bo}~qo2!$j#sVH!Z9z3YHzlRCS!vcgvMjHjn`)^ zdd6?9{Eu&+CZ0tjhUE7LV^$oSt#Eu^+Rrut8AC?zZMVtIO5Ebhnjtt&(4NfocxBly z63V>~ySx0#dzU*HXLb}IA?BZdlMTc7`gF3dlq(FKCU3z$)cBGgF_?aSwYVnrK#(P@ zW^)Ob`;H194`lznF&l%Arcjmko|6@Y%+g1|NwvLvk0P|GePol&`huKcTVJoC6+VfJ@vwRpkEzZH=my_2x;oL;RWQXKa$e}}j3MOKuZ1`ddE zUz?EFL1LdT1`WcfUX>=|ctp&^o9yvRBy_)y7vZ?9=Hj>*dN^w@1mf<{``3HBk>bZC zqn~{hKat*%b2ps#IFUvq2zjiv*7+J6t07l#)~F4q{~f|pAxOmnyO9>FmwElDUdIJ` zT=ArC1&O=HY{|hTSzpVT8%{r!TwpQM@sGCdd7j~`3&2}jNo5lIOLW6v291?K{(KK< zfA7A4VWS1=vlXky#(9Q}yiMO!Bg)oAJ@odx!FOJTrzHpq57$RcHRqL<$aE~qTq3Ge zB-Y^)q8LG8$C8XJ?R!|O>Z>$mTd)0%BZ26uFBD%!i%!2LyuF&w!_0l@*xdU7?|QD< zYr(p+QVlia2I>!mIvG2QNrJRf?Z2?!t!@1cOKq&&k5M_U&99TSP|->Y^$}$+8;QY9 z`7n+-Ah;S=AevFnD)HoCFkzoHVZe!Qv+Y!n zy+k>Kr%lpkuJM{VAwFO`5i{r#9sBKa(6;*Mq}5l0ob;RXElwUXf&kAGdCf=dTfNN} ztSMo^fhePQCquj2DZRHp_og1Ube}WbEBnR@)w9Dn-G7ZE9izYQUzCU>?DB)OD&7H7 z4IpIjUvPmDhx3o2e+`>rZ{G#_9qo>a9=`XFMHiEt?lV(8IBn?-!U%v(^JZ_K?@Azc z(>2a`wSm>0&t`%lS_*C$%yY6oAE)qN7s}MlZJH3`sNUIqZ`rD3_h$J@Tv=2|+hM9k<)DR^*-TANEiV2^m|qeU`j8wsMq?*7O{OL~L+6Zq zT|VtIQy8+~W*y#;hT{?i7w@tt=Hu*07K$Hmo0@wSMTaqS-xLs=Pq5)zCjazf(iq@4 z9^By<8pdrSE9T{8Jb}G>HuzEj^pe?54|X&f)rMs^y;yT#DsOFPaD$ue3v4=?$oN6T6cexF=1yyAslh|nZ3b>)Noa)k9X~a7}IDk(l!d!o1V85%* zN3(;U=O12uC|@xa9x%=h9|<1$bMhj-;15lC%(Am;R7Q&A*HV>yla8I2W!k1MhS>8L zs0R&<5a);Q4Qvwx&HV}2Y71#G%S78F|jH)2;Jd>*P|JG6urxIr+|#Vd_(DU*mYyA z-m9Rlyg0&_$d0#WTutdcR@^Ob{I@nk1CJZrI+xGn=#+hf-RSh3B5rc&>g$ug1PbMq zqVV?mxLQmveP$epg zk2ylExN)?2T;~{$?-Kk04(eQrl}VC z#g>+tb5QKo);_qQ#>9q|^+u`rkLD@SwdU;&X#B7)-GFR7xM1y-__}pW=hKI2)8@W zq1I1^gYgXuXk}C}H|_lG^_u4Srp9}il=ANLf}=}Uvt}|2C`WSMi=}Jto2r*ZT;h=S zh21A>=`q>X;D@uKzEC;3lYni56})A%;zZC2{`TvR#}bEt(HPnc{reTIO-C>sztOJI zU?j9U)aB=sdUZlc6cb;DfLFXS-l*p1@7GmD)86K#Mt2W{;_x>Q4@N{(ovXxZVep4? zh8>+xeXN{NtGKWiN~>uW3B&{Yc;zS-ejMOexCMMHL9LIn72r7RAA7U6=3zTVV7T}o zRIW3tZd{6F%$n)?&(nh}9FiP){+AAwDC6dQM7*unT8rNU$X(1S#8fo{< z+Rb%hWdHu7#3+=aa}?I!$p0yCBVW9&taVqQAAbDs?62F^JNQHC-a3qJdB5i+rd2e2 z)+HRK+5LVWIa-0YG4y$!?sYNX{7-VHYdGA0W$e@A;mqqFY}lMG6Yh-pRBLT7j90Cw zyG^en$44h0^Z*Pep<-yvO-o@%4*k&GPqXV0dwvM{a!NH`-7fnckH#`~@ zMSkBk!Vk>?#<88tQXM#YB~?G{ktNb*OOj*HOHOB}7nu*kjQ90ygF|B*x^bz_TgNb_ zHQB-*OF{>3&Nq9gw7q47B!;rPHsPT+vLkeL07YJc&nE=q=e+zmdwYXM$*YHh#p;Ld zF-;ML8MPrzI1V`jqddZkI8LaK{C9CFWtGTJmHs-rC4;LtPLQKwjk{x-Evgbpnk4XA zknOU?QJ4922uzN2zlxtvEXen06}KQJbm>FLd)boxg3kDBfPi7K7FC+SxnIF3)KuUv zm;Zp*C;?eW`%K?dJArQgO^}nIv}-Yr=wu^orf;`9XS4G5?+kR#2qVC(Xwz`I+O5s8uHh?zE9kCjBtdYu$xQlpmi^K0J_a&99)okKWx~Tn{Bo4IL{) zQBauk@Zx17Rv+xJCu#Gk^(DK4ZA`~6{L5irCz6+OIBc^+u?ofek>WCBCaQot>v;_s zXnk0hEyPMZ`BF6~T$fsi^6u}7m>93N(p8#>ZT(e9e$@T+IZ!_p zPMz9!DoEyS-S!S;sk!t5IZbRT1 z^0NYrPn%!*_?b*I?~At|Q*;>5i1ua*KXYBS1+mg0@nS=a+|VZyb(M$XypeKjRxpr^mg!IehvjM zDkD;YX~TYU3uz!mS=hU%qM9A@`!g%yw7TuoXJNQH+BTqf&n<8<>yI^oVo}#3SW3((7{^IA_u~R zgy1Kff|ywU9X=6AFH2pNsJ;5f@9vs+eAI1|6ToT7_K)ga011RtHqDOUp$|utJ$M<3 zXd#X9vOJ9d*kU!;&0hD9wcbWk-XJLui`N0`-^swJsv_rHQ;gET-c~!Ye}-2ecn1)V ztK_gat`PmJ9bwc!_U{UU)&AesGQ`pUFx(@?`#%Zq|2v&i8>c2p&0=$*z#*zwIb*_X zi2;q5zAg|8n2fW|lsOnmxy-(#UH-&mRoCw22`eT(UpDMXUI?n-JZBPjruF{ma>Isy z(&2`vTk(&d2*L37|9(Iq^%nm>+Wb!alV&!$MK*mQd3?X#P7=`^F+jw(0XoA$ z1t@OxSUAKepq)x3v&u<#1aBA5zU$-b$@kTUe*ki-f3tIhNfLId?BGu>ncX!sK#F4e z$89`weeSPzZ$A08HLqX{vP4|jtk0j&;_0Wi;scmMzX!9;=Lp;Q<0?qV1v^YLkzC8{ zKN%fsW)krt>wCal}Hhx_7R!#t{tde11d~_-T=;-?e3KMC%KPX-j?5BIZ_U?ur z%WoI&M>60&39kWj4igz$?CY%n$I{3+EahPJV|Hw7|Dq(49Wl(ln~2-`7T`W3fPd@3 zd{9mr_>+&06N4%21>Hun_p?P6V15>Ly5MzaxqFpvi%<9HivQF9~2owYj#2|CF4Rx-eYrabVkq~;pP(%00 za348A%K!ixypotF&YFl#AjK8x%_})NX(tEvTBQr zBI@dD(Z48}<;#3TQY;j!OQe{y&~nSnVNO3uy&z|3Trl|X$SuuhP$}>?Y6o#$sP0rGgfIwWS`a5g^uPy+T70hT{!Lsk;#dKAj6o z{e|O(2g!CL)ojk+@2oaU79Y9wUo06_sR#dI01(m#6QU}TU2VUmC6RhRGf;We19=sK z;O06QBE0-#pqfo}*Wo8^A25VG8W=0CH5Ng6qQTDJ7Xp3vYZPf|UWYgPT*S4>T5G6V z-h}#Y;Jx#hKUeyp{&_++n&)0EiXJV_E&*+WLxEvx3-GO>bq(HLm7Bz8lMgcyvEl*v zAjUqW)k@p6nB8UQ8cS5d+XsBNTQm9k@l*lc zIg2lLMgzf(+ClfQgz{T6f35EJ$@ZZ%n-iA3N&K%CMS?Rk&tX9ZbBvwXku8j#uYL}s zIDZsZBE*C}KQ9i7IH#TjXlvYDC4L2GRP?W5d`WAj`j$ta3}eq^eWA|JufoJM=8PU- zCQUM-CR4yhV8g&YWvjJXd`!8gV8mpOu7wi4PPr=1#jB*2C zKP&VnMb^4uVc6LSTx?g5k4EbnqPNAF@d;-+fO#?pIHY^;Cwdb<#M@vzTU>uuKK+R82N|m>ZU#uC~iMwev)P!Pf2%B_&DX(YuAU7UZLMVb-9n|>q%R{5GZR`*(=C1U7{S}{^9NE zz!S9StuVq!2$61jFMTe@x2knv5G&~mPn7i=4wxWKn`A$Nbq(0d7LVO^n_{bi!?5!0 z&Z)mfQ8^#IqxQ%@UEqla@G2n?c})AT{USg45E4Km9Z!?Ko=KxNJO;8R3$GT;v(&y; zXLfk}EM$Azl2~#vyt@AtIar1VV1&IPAwi`ohPvz)vCynR3-B2~Vjj$HVtpqs)rIk8*k2k6wo|k;v+NG97QXT5t2h7wrx*2`|S)(scw|ov}MWwpfZe<+>zAU zf_@l1lJlMnb!>t4r$(5+lF@!!r~)4z%p@@JF%9LmIqe^VIEJBwP*v zA}?yRqc}cVz;+Z)65d}54?`j{NG~j@>zNxhRY36l7WSczEj3(J7nQAZ{(NzpXU0vV zbtEG0Gp~ZA(w*zYIWZpAQ_~(IO9hXit5KP745l4=rkr8>1I`e4lV^WKqyRFwJb%Uo~xwoDe zW047qUk5QEm0K$)a=D>uGk{k(KS{P+;{9ou4Vp+RiRDryb)qfKzg!W@=@#VVzFQHR zTh$h_-sjfk8LnN9$PADYbX|)~HsJ)CUt{3-taDPV{m>cdG26X;05C;m*nkBa<$t{n zhq(E{YL^TE;#oCL%pORXoT}=8xTH-aR5gNm^5JM z`Vmboxa_u!NPf?o7t_m3V5O|IIl_~VX#*wsH#(i8q5MxXy*FT{!ddAUQ zRJ%rFl0+83Bx>RN8yVuwv?e_*!K^pTXGt;Q%yVhQzfk}n`UqRd?46REb?svJ=8AW`e%J)Y>ufiLj)& zbkJb!;QntH%N8b3L$AI)z$zMK9Alwf=)bD7KkPA-)wc=TB8 zL8lPgQ~RL@0OBKUFyQXDF%FMvijo3<0Z4}|%1z?iaJ~=cYem)MP=pB(BuyS*oIUWN zT6vt{GF4No?yvqHYXS62A)&H@A}l0HZJHnkw$Un$sdx*RH|i~tYv#kzFZp(K79T&$ z5Ax+FZlocztYL$}{TJyZ1XP8uNxsQN=qd(bhduPX3fG0C)N<&WYlH%qln663zAnKEU!JFdZo>1&R7oK@nq*$j?Qa&Ak(? zznzEgw=Ry{3j;g{ADrMwWA-_b5kKBPzlY0)dY>ad=#0S*4N!*Q=-khMy);tZsVp^* ji_1GNvQyDuIErY{$Cbo6)<1|!hX55NuwsM!+lc=KxkfOZ literal 0 HcmV?d00001 From 45f7763b0fd81e833a7c3b9d50074413095d5e4c Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 13 Jun 2025 14:30:01 -0400 Subject: [PATCH 15/25] add openmp to cmake --- raster/r.mapcalc/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/raster/r.mapcalc/CMakeLists.txt b/raster/r.mapcalc/CMakeLists.txt index ab9155ce0bb..e851137e2b5 100644 --- a/raster/r.mapcalc/CMakeLists.txt +++ b/raster/r.mapcalc/CMakeLists.txt @@ -30,7 +30,8 @@ build_program( OPTIONAL_DEPENDS Readline::Readline Readline::History - Threads::Threads) + Threads::Threads + OPENMP) build_program( NAME From b2ba22fa14b9aa689109cd4aa9a51a8a43ab2e05 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Fri, 13 Jun 2025 14:50:26 -0400 Subject: [PATCH 16/25] use new functions for setting openmp threads --- raster/r.mapcalc/main.c | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 6c28d08e2ab..9b581dc11b2 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -184,31 +184,24 @@ int main(int argc, char **argv) } pre_exec(); - threads = atoi(nprocs->answer); -/* Determine the number of threads */ -#if defined(_OPENMP) - if ((strcmp(argv[0], "r3.mapcalc") == 0) && (threads > 1)) { + /* Determine the number of threads */ + threads = atoi(nprocs->answer); + if ((strcmp(argv[0], "r3.mapcalc") == 0) && (threads != 1)) { threads = 1; G_verbose_message(_("r3.mapcalc does not support parallel execution.")); } - else if (seeded) { + else if ((seeded) && (threads != 1)) { threads = 1; G_verbose_message( _("Parallel execution is not supported for random seed.")); } else { - /* make the maximum num of threads be 4 to avoid I/O overheads from too - * many threads */ - threads = threads > 4 ? 4 : threads; - G_verbose_message(_("The default number of threads is %d."), threads); + threads = G_set_omp_num_threads(nprocs); + threads = Rast_disable_omp_on_mask(threads); + if (threads < 1) + G_fatal_error(_("<%d> is not valid number of nprocs."), threads); } - omp_set_num_threads(threads); - G_verbose_message(_("Compute in paralell. Number of threads set to %d."), - threads); -#else - G_verbose_message(_("Compute in searial. OpenMP is not available.")); -#endif execute(result); post_exec(); From 42b3f97f45468bb5f57d8bfa7b786b0cfb091c55 Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sat, 14 Jun 2025 11:34:43 -0400 Subject: [PATCH 17/25] correct the location and the condition to free nprocs->answer --- raster/r.mapcalc/main.c | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 9b581dc11b2..add6d18a3aa 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -199,6 +199,8 @@ int main(int argc, char **argv) else { threads = G_set_omp_num_threads(nprocs); threads = Rast_disable_omp_on_mask(threads); + if (strcmp(nprocs->answer, "0")) + G_free(nprocs->answer); if (threads < 1) G_fatal_error(_("<%d> is not valid number of nprocs."), threads); } @@ -208,10 +210,6 @@ int main(int argc, char **argv) all_ok = 1; - // Free nprocs->answer if it is not the default value, "1" - if (threads > 1) { - G_free(nprocs->answer); - } G_free(p); p = NULL; From 8746bbbf994863f7f17a49757898e3f18177d95e Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sun, 15 Jun 2025 18:08:11 -0400 Subject: [PATCH 18/25] Ensure r3 is detected and ensure correct #threads is assigned --- raster/r.mapcalc/main.c | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index add6d18a3aa..0a409c89994 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -187,24 +187,31 @@ int main(int argc, char **argv) /* Determine the number of threads */ threads = atoi(nprocs->answer); - if ((strcmp(argv[0], "r3.mapcalc") == 0) && (threads != 1)) { + + /* Check if the program name is r3.mapcalc */ + const char *progname = strrchr(argv[0], '/'); + progname = progname ? progname + 1 : argv[0]; + + if ((strcmp(progname, "r3.mapcalc") == 0) && (threads != 1)) { threads = 1; + nprocs->answer = "1"; G_verbose_message(_("r3.mapcalc does not support parallel execution.")); } else if ((seeded) && (threads != 1)) { threads = 1; + nprocs->answer = "1"; G_verbose_message( _("Parallel execution is not supported for random seed.")); } - else { - threads = G_set_omp_num_threads(nprocs); + + /* Ensure the proper number of threads is assigned */ + threads = G_set_omp_num_threads(nprocs); + if (threads > 1) threads = Rast_disable_omp_on_mask(threads); - if (strcmp(nprocs->answer, "0")) - G_free(nprocs->answer); - if (threads < 1) - G_fatal_error(_("<%d> is not valid number of nprocs."), threads); - } + if (threads < 1) + G_fatal_error(_("<%d> is not valid number of nprocs."), threads); + /* Execute calculations */ execute(result); post_exec(); From bed8011c9a2971e8e740579433fda14b099425bb Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sun, 15 Jun 2025 19:51:06 -0400 Subject: [PATCH 19/25] Handle Windows path separators --- raster/r.mapcalc/main.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 0a409c89994..453fa60a615 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -189,9 +189,14 @@ int main(int argc, char **argv) threads = atoi(nprocs->answer); /* Check if the program name is r3.mapcalc */ + /* Handle both Unix and Windows path separators */ const char *progname = strrchr(argv[0], '/'); + const char *win_progname = strrchr(argv[0], '\\'); + if (!progname || (win_progname && win_progname > progname)) + progname = win_progname; progname = progname ? progname + 1 : argv[0]; + if ((strcmp(progname, "r3.mapcalc") == 0) && (threads != 1)) { threads = 1; nprocs->answer = "1"; From dfe3a1bdc2323aac871e6ab57672cde23028f16d Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sun, 15 Jun 2025 19:53:20 -0400 Subject: [PATCH 20/25] remove redundant empty line Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- raster/r.mapcalc/main.c | 1 - 1 file changed, 1 deletion(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 453fa60a615..fe993c21a2d 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -196,7 +196,6 @@ int main(int argc, char **argv) progname = win_progname; progname = progname ? progname + 1 : argv[0]; - if ((strcmp(progname, "r3.mapcalc") == 0) && (threads != 1)) { threads = 1; nprocs->answer = "1"; From 895a7d9e1ed68fe7c1722c7b1931fdd26e111fd6 Mon Sep 17 00:00:00 2001 From: Chung-Yuan Liang <77927944+cyliang368@users.noreply.github.com> Date: Sun, 15 Jun 2025 19:53:42 -0400 Subject: [PATCH 21/25] remove space Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- raster/r.mapcalc/main.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index fe993c21a2d..b83adc32659 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -192,7 +192,7 @@ int main(int argc, char **argv) /* Handle both Unix and Windows path separators */ const char *progname = strrchr(argv[0], '/'); const char *win_progname = strrchr(argv[0], '\\'); - if (!progname || (win_progname && win_progname > progname)) + if (!progname || (win_progname && win_progname > progname)) progname = win_progname; progname = progname ? progname + 1 : argv[0]; From 8ab77d2219eaeb4fa88605efcc96565c0f6fc57c Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Tue, 24 Jun 2025 12:05:07 -0400 Subject: [PATCH 22/25] test1 --- raster/r.mapcalc/main.c | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index b83adc32659..d3b849decee 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -191,15 +191,19 @@ int main(int argc, char **argv) /* Check if the program name is r3.mapcalc */ /* Handle both Unix and Windows path separators */ const char *progname = strrchr(argv[0], '/'); - const char *win_progname = strrchr(argv[0], '\\'); - if (!progname || (win_progname && win_progname > progname)) - progname = win_progname; + if (!progname) + progname = strrchr(argv[0], '\\'); progname = progname ? progname + 1 : argv[0]; + // const char *progname = strrchr(argv[0], '/'); + // const char *win_progname = strrchr(argv[0], '\\'); + // if (!progname || (win_progname && win_progname > progname)) + // progname = win_progname; + // progname = progname ? progname + 1 : argv[0]; - if ((strcmp(progname, "r3.mapcalc") == 0) && (threads != 1)) { + if ((strcmp(progname, "r3.mapcalc", 10) == 0) && (threads != 1)) { threads = 1; nprocs->answer = "1"; - G_verbose_message(_("r3.mapcalc does not support parallel execution.")); + G_warning(_("r3.mapcalc does not support parallel execution.")); } else if ((seeded) && (threads != 1)) { threads = 1; From a4862b25cae98c0ffff58de9779b059e422845bf Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Tue, 24 Jun 2025 12:43:46 -0400 Subject: [PATCH 23/25] fix --- raster/r.mapcalc/main.c | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index d3b849decee..051a2591162 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -194,13 +194,8 @@ int main(int argc, char **argv) if (!progname) progname = strrchr(argv[0], '\\'); progname = progname ? progname + 1 : argv[0]; - // const char *progname = strrchr(argv[0], '/'); - // const char *win_progname = strrchr(argv[0], '\\'); - // if (!progname || (win_progname && win_progname > progname)) - // progname = win_progname; - // progname = progname ? progname + 1 : argv[0]; - if ((strcmp(progname, "r3.mapcalc", 10) == 0) && (threads != 1)) { + if ((strncmp(progname, "r3.mapcalc", 10) == 0) && (threads != 1)) { threads = 1; nprocs->answer = "1"; G_warning(_("r3.mapcalc does not support parallel execution.")); From cbb864ea4e883fd0f58627ab8203c10c1a7c08df Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Tue, 24 Jun 2025 16:35:07 -0400 Subject: [PATCH 24/25] test nmedian with nprocs=1 --- raster/r.mapcalc/main.c | 2 +- .../testsuite/test_nmedian_bug_3296.py | 35 +++++++++++++++++++ 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/raster/r.mapcalc/main.c b/raster/r.mapcalc/main.c index 051a2591162..faa2fed12c1 100644 --- a/raster/r.mapcalc/main.c +++ b/raster/r.mapcalc/main.c @@ -198,7 +198,7 @@ int main(int argc, char **argv) if ((strncmp(progname, "r3.mapcalc", 10) == 0) && (threads != 1)) { threads = 1; nprocs->answer = "1"; - G_warning(_("r3.mapcalc does not support parallel execution.")); + G_verbose_message(_("r3.mapcalc does not support parallel execution.")); } else if ((seeded) && (threads != 1)) { threads = 1; diff --git a/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py b/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py index 3ab7d7c255f..c5a8e8f4690 100644 --- a/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py +++ b/raster/r.mapcalc/testsuite/test_nmedian_bug_3296.py @@ -125,6 +125,41 @@ def test_dcell(self): actual=self.output, reference=self.output_ref, precision=0 ) + def test_cell_nprocs1(self): + expression = "{o}=nmedian(({i}[0,-1] - {i})^2,({i}[0,1] - {i})^2)".format( + o=self.output, i=self.input + ) + self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True) + self.assertRasterExists(self.output) + self.to_remove.append(self.output) + self.assertRastersNoDifference( + actual=self.output, reference=self.output_cell, precision=0 + ) + + def test_fcell_nprocs1(self): + expression = ( + "{o}=nmedian(float(({i}[0,-1] - {i})^2), float(({i}[0,1] - {i})^2))".format( + o=self.output, i=self.input + ) + ) + self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True) + self.assertRasterExists(self.output) + self.to_remove.append(self.output) + self.assertRastersNoDifference( + actual=self.output, reference=self.output_ref, precision=0 + ) + + def test_dcell_nprocs1(self): + expression = "{o}=nmedian(double(({i}[0,-1] - {i})^2), double(({i}[0,1] - {i})^2))".format( + o=self.output, i=self.input + ) + self.assertModule("r.mapcalc", expression=expression, nprocs=1, overwrite=True) + self.assertRasterExists(self.output) + self.to_remove.append(self.output) + self.assertRastersNoDifference( + actual=self.output, reference=self.output_ref, precision=0 + ) + if __name__ == "__main__": test() From dd0468defb383685c372ee8905bdcfa5a4b96105 Mon Sep 17 00:00:00 2001 From: Anna Petrasova Date: Wed, 25 Jun 2025 10:23:51 -0400 Subject: [PATCH 25/25] remove static variable --- lib/calc/xnmedian.c | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/lib/calc/xnmedian.c b/lib/calc/xnmedian.c index 3dab302e90e..e10efbce1ac 100644 --- a/lib/calc/xnmedian.c +++ b/lib/calc/xnmedian.c @@ -44,8 +44,7 @@ static int dcmp(const void *aa, const void *bb) int f_nmedian(int argc, const int *argt, void **args) { - static void *array; - static int alloc; + void *array; int size = argc * Rast_cell_size(argt[0]); int i, j; @@ -56,10 +55,7 @@ int f_nmedian(int argc, const int *argt, void **args) if (argt[i] != argt[0]) return E_ARG_TYPE; - if (size > alloc) { - alloc = size; - array = G_realloc(array, size); - } + array = G_malloc(size); switch (argt[0]) { case CELL_TYPE: {