From 46cf0db38860fb392ecb5d87364acfd3ad411253 Mon Sep 17 00:00:00 2001 From: John Halley Gotway Date: Tue, 7 Nov 2023 12:42:23 -0700 Subject: [PATCH] Feature #2724 mode_openmp (#2726) --- docs/Users_Guide/config_options.rst | 12 +- src/basic/vx_util/data_plane.cc | 77 +- src/basic/vx_util/vx_util.h | 1 + src/libcode/vx_shapedata/engine.cc | 731 +++++++----------- src/libcode/vx_shapedata/interest.cc | 67 +- src/libcode/vx_shapedata/interest.h | 8 +- src/libcode/vx_shapedata/shapedata.cc | 242 +++--- src/tools/core/ensemble_stat/ensemble_stat.cc | 5 - src/tools/core/mode/mode.cc | 5 + 9 files changed, 490 insertions(+), 658 deletions(-) diff --git a/docs/Users_Guide/config_options.rst b/docs/Users_Guide/config_options.rst index 35805c71d3..92470cfda5 100644 --- a/docs/Users_Guide/config_options.rst +++ b/docs/Users_Guide/config_options.rst @@ -459,12 +459,20 @@ Regions of parallelized code are: * :code:`fractional_coverage (data_plane_util.cc)` + * Called by `gen_ens_prod` to compute NMEP outputs. + * Called by `grid_stat` when applying neighborhood verification methods. + + * :code:`ShapeData::conv_filter_circ() (shapedata.cc)` + + * Called by `mode` to apply a convolution smoothing operation when + defining objects. + Only the following top-level executables can presently benefit from OpenMP parallelization: * :code:`grid_stat` - * :code:`ensemble_stat` * :code:`grid_ens_prod` + * :code:`mode` **Thread Binding** @@ -474,7 +482,7 @@ guarantees that threads remain evenly distributed across the available cores. Otherwise, the operating system may migrate threads between cores during a run. OpenMP provides some environment variables to handle this: :code:`OMP_PLACES` -and :code:`OMP_PROC_BIND`. We anticipate that the effect of setting only +and :code:`OMP_PROC_BIND`. We anticipate that the effect of setting only :code:`OMP_PROC_BIND=true` would be neutral-to-positive. However, there are sometimes compiler-specific environment variables. Instead, diff --git a/src/basic/vx_util/data_plane.cc b/src/basic/vx_util/data_plane.cc index 2dd2dcaf38..52d2b72f9e 100644 --- a/src/basic/vx_util/data_plane.cc +++ b/src/basic/vx_util/data_plane.cc @@ -23,6 +23,7 @@ using namespace std; #include +#include #include "data_plane.h" @@ -230,57 +231,57 @@ void DataPlane::dump(ostream & out, int depth) const { void DataPlane::debug_examine(bool show_all_values) const { - vector values; - vector count; + // Nothing to print if verbosity level is below 4 + if(mlog.verbosity_level() < 4) return; + + map value_count_map; int total_count = 0; - - for (int x=0; x::iterator vi; - vi = find(values.begin(), values.end(), v); - if (vi == values.end()) { - values.push_back(v); - count.push_back(1); - } else { - int ii = vi - values.begin(); - count[ii] = count[ii] + 1; - } + + for(int i=0; i 0) total_count++; + + if (show_all_values) { + + // Add a new map entry + if(value_count_map.count(Data[i]) == 0) { + value_count_map[Data[i]] = 0; } + + // Increment count + value_count_map[Data[i]] += 1; } } - if (show_all_values) { - for (size_t i=0; i 0 = " << total_count + << " of " << Data.size() << "\n"; + + return; } /////////////////////////////////////////////////////////////////////////////// string DataPlane::sdebug_examine() const{ + ConcatString cs; + int n = 0; - int total_count = 0; - - for (int x=0; x 0) n++; } - char buf[1000]; - sprintf(buf, "Total count = %d", total_count); - string retval = buf; - return retval; + + cs << "Total count > 0 = " << n << " of " << (int) Data.size(); + + return cs; } /////////////////////////////////////////////////////////////////////////////// diff --git a/src/basic/vx_util/vx_util.h b/src/basic/vx_util/vx_util.h index 84638028ff..c50e7d16fc 100644 --- a/src/basic/vx_util/vx_util.h +++ b/src/basic/vx_util/vx_util.h @@ -37,6 +37,7 @@ #include "fix_float.h" #include "get_filenames.h" #include "grib_constants.h" +#include "GridTemplate.h" #include "int_array.h" #include "interp_mthd.h" #include "interp_util.h" diff --git a/src/libcode/vx_shapedata/engine.cc b/src/libcode/vx_shapedata/engine.cc index d16ba7de6c..0bfab73929 100644 --- a/src/libcode/vx_shapedata/engine.cc +++ b/src/libcode/vx_shapedata/engine.cc @@ -6,7 +6,6 @@ // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* - /////////////////////////////////////////////////////////////////////// using namespace std; @@ -17,6 +16,8 @@ using namespace std; #include #include #include +#include +#include #include "engine.h" #include "mode_columns.h" @@ -591,12 +592,17 @@ void ModeFuzzyEngine::do_fcst_convolution() { // Apply a circular convolution to the raw field // - mlog << Debug(4) << " Before fcst convolution: " << fcst_conv->sdebug_examine() << "\n"; - + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " Before fcst convolution: " + << fcst_conv->sdebug_examine() << "\n"; + } if(r > 0) fcst_conv->conv_filter_circ(2*r + 1, conf_info.Fcst->vld_thresh); - mlog << Debug(4) << " After fcst convolution: " << fcst_conv->sdebug_examine() << "\n"; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After fcst convolution: " + << fcst_conv->sdebug_examine() << "\n"; + } need_fcst_conv = false; need_fcst_thresh = true; @@ -625,14 +631,20 @@ void ModeFuzzyEngine::do_obs_convolution() { mlog << Debug(3) << "Applying circular convolution of radius " << r << " to the observation field.\n"; - mlog << Debug(4) << " Before obs convolution: " << obs_conv->sdebug_examine() << "\n"; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " Before obs convolution: " + << obs_conv->sdebug_examine() << "\n"; + } // // Apply a circular convolution to the raw field // if(r > 0) obs_conv->conv_filter_circ(2*r + 1, conf_info.Obs->vld_thresh); - mlog << Debug(4) << " After obs convolution: " << obs_conv->sdebug_examine() << "\n"; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After obs convolution: " + << obs_conv->sdebug_examine() << "\n"; + } need_obs_conv = false; need_obs_thresh = true; @@ -657,13 +669,21 @@ void ModeFuzzyEngine::do_fcst_thresholding() { *fcst_thresh = *fcst_raw; fcst_thresh->threshold(conf_info.Fcst->conv_thresh); - mlog << Debug(4) << " After thresholding raw fcst field: " << fcst_thresh->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After thresholding raw fcst field: " + << fcst_thresh->sdebug_examine() << "\n"; + } // // Threshold the convolved field // fcst_mask->threshold(conf_info.Fcst->conv_thresh); - mlog << Debug(4) << " After thresholding convolved fcst field: " << fcst_mask->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After thresholding convolved fcst field: " + << fcst_mask->sdebug_examine() << "\n"; + } mlog << Debug(3) << "Applying convolution threshold " << conf_info.Fcst->conv_thresh.get_str() @@ -692,13 +712,21 @@ void ModeFuzzyEngine::do_obs_thresholding() { *obs_thresh = *obs_raw; obs_thresh->threshold(conf_info.Obs->conv_thresh); - mlog << Debug(4) << " After thresholding raw obs field: " << obs_thresh->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After thresholding raw obs field: " + << obs_thresh->sdebug_examine() << "\n"; + } // // Threshold the convolved field // obs_mask->threshold(conf_info.Obs->conv_thresh); - mlog << Debug(4) << " After thresholding convolved obs field: " << obs_mask->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After thresholding convolved obs field: " + << obs_mask->sdebug_examine() << "\n"; + } mlog << Debug(3) << "Applying convolution threshold " << conf_info.Obs->conv_thresh.get_str() @@ -733,7 +761,11 @@ void ModeFuzzyEngine::do_fcst_filtering() { fcst_raw, conf_info.Fcst->conv_thresh, grid, conf_info.Fcst->var_info->is_precipitation()); - mlog << Debug(4) << " After attribute filtering of fcst field: " << fcst_mask->sdebug_examine() << "\n"; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After attribute filtering of fcst field: " + << fcst_mask->sdebug_examine() << "\n"; + } + mlog << Debug(3) << "Applying object attribute filtering" << " resulted in " << fcst_mask->n_objects() << " simple forecast objects.\n"; @@ -769,7 +801,11 @@ void ModeFuzzyEngine::do_obs_filtering() { obs_raw, conf_info.Obs->conv_thresh, grid, conf_info.Obs->var_info->is_precipitation()); - mlog << Debug(4) << " After attribute filtering of obs field: " << obs_mask->sdebug_examine() << "\n"; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After attribute filtering of obs field: " + << obs_mask->sdebug_examine() << "\n"; + } + mlog << Debug(3) << "Applying object attribute filtering" << " resulted in " << obs_mask->n_objects() << " simple observation objects.\n"; @@ -797,7 +833,11 @@ void ModeFuzzyEngine::do_fcst_splitting() { if(!need_fcst_split) return; *fcst_split = split(*fcst_mask, n_fcst); - mlog << Debug(4) << " After splitting of fcst field: " << fcst_split->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After splitting of fcst field: " + << fcst_split->sdebug_examine() << "\n"; + } need_fcst_split = false; need_fcst_merge = true; @@ -816,7 +856,11 @@ void ModeFuzzyEngine::do_obs_splitting() { if(!need_obs_split) return; *obs_split = split(*obs_mask, n_obs); - mlog << Debug(4) << " After splitting of obs field: " << obs_split->sdebug_examine() << "\n"; + + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After splitting of obs field: " + << obs_split->sdebug_examine() << "\n"; + } need_obs_split = false; need_obs_merge = true; @@ -1033,32 +1077,22 @@ void ModeFuzzyEngine::do_matching() { void ModeFuzzyEngine::do_no_match() { int j, k, n; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * obs_shape = (ShapeData *) 0; + ShapeData cur_shape; do_fcst_splitting(); do_obs_splitting(); clear_colors(); - fcst_shape = new ShapeData [n_fcst]; - obs_shape = new ShapeData [n_obs]; - - if(!fcst_shape || !obs_shape) { - - mlog << Error << "\nModeFuzzyEngine::do_no_match() -> " - << "memory allocation error\n\n"; - exit(1); - } - // // Do the single features // fcst_single.set_size(n_fcst); for(j=0; jvar_info->is_precipitation()); fcst_single[j].object_number = j+1; @@ -1067,8 +1101,9 @@ void ModeFuzzyEngine::do_no_match() { obs_single.set_size(n_obs); for(j=0; jvar_info->is_precipitation()); obs_single[j].object_number = j+1; @@ -1112,9 +1147,6 @@ void ModeFuzzyEngine::do_no_match() { // Done // - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] obs_shape; obs_shape = (ShapeData *) 0; - return; } @@ -1123,32 +1155,22 @@ void ModeFuzzyEngine::do_no_match() { void ModeFuzzyEngine::do_match_merge() { int j, k, n; InterestInfo junkinfo; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * obs_shape = (ShapeData *) 0; + ShapeData cur_shape; do_fcst_splitting(); do_obs_splitting(); clear_colors(); - fcst_shape = new ShapeData [n_fcst]; - obs_shape = new ShapeData [n_obs]; - - if(!fcst_shape || !obs_shape) { - - mlog << Error << "\nModeFuzzyEngine::do_match_merge() -> " - << "memory allocation error\n\n"; - exit(1); - } - // // Do the single features // fcst_single.set_size(n_fcst); for(j=0; jvar_info->is_precipitation()); fcst_single[j].object_number = j+1; @@ -1157,8 +1179,9 @@ void ModeFuzzyEngine::do_match_merge() { obs_single.set_size(n_obs); for(j=0; jvar_info->is_precipitation()); obs_single[j].object_number = j+1; @@ -1277,9 +1300,6 @@ void ModeFuzzyEngine::do_match_merge() { // Done // - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] obs_shape; obs_shape = (ShapeData *) 0; - return; } @@ -1292,12 +1312,10 @@ void ModeFuzzyEngine::do_match_merge() { /////////////////////////////////////////////////////////////////////// void ModeFuzzyEngine::do_fcst_merge_thresh() { - int j, k, x, y; - int n_fcst_merge, intersection; - int count, first_k; + const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> "; + int j, mid, oid; + int n_fcst_merge; ShapeData fcst_merge_mask, fcst_merge_split; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * fcst_merge_shape = (ShapeData *) 0; do_fcst_splitting(); @@ -1318,90 +1336,60 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() { fcst_merge_split = split(fcst_merge_mask, n_fcst_merge); // - // Allocate space for all of the simple forecast shapes and - // forecast merge shapes + // Mapping of simple object id to merge object id values and vice-versa // - fcst_shape = new ShapeData [n_fcst]; - fcst_merge_shape = new ShapeData [n_fcst_merge]; + map > obj_to_mrg_map; + map > mrg_to_obj_map; + std::set empty_set; - if(!fcst_shape || !fcst_merge_shape) { + // Add map entries for each object + for(j=1; j<=n_fcst; j++) obj_to_mrg_map[j] = empty_set; + for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set; - mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> " - << "memory allocation error\n\n"; - exit(1); - } + // Loop over all the data points + for(j=0; jdata.nxy(); j++) { + + oid = fcst_split->data.buf()[j]; + mid = fcst_merge_split.data.buf()[j]; + + // Update object id maps + if(oid > 0) obj_to_mrg_map[oid].insert(mid); + if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid); + + } // end for j // - // Select all of the simple forecast shapes and - // forecast merge shapes + // Each simple object should be fully contained in a merge object // - for(j=0; j= fcst area - // - - if(intersection >= fcst_shape[k].area()) { - - collection.make_room(); - - // - // Keep track of the first embedded shape. You only want to - // create a composite if there are more than 1 shapes in it. - // - if(count == 0) first_k = k+1; + // Ignore merge objects containing a single simple object + if(mrg_to_obj_map[mid].size() <= 1) continue; - else if(count == 1) { - collection.set[collection.n_sets].add_pair(first_k, -1); - collection.set[collection.n_sets].add_pair(k+1, -1); - } - else { - collection.set[collection.n_sets].add_pair(k+1, -1); - } + // Add all simple objects into one new set + collection.make_room(); + for(auto it = mrg_to_obj_map[mid].begin(); + it != mrg_to_obj_map[mid].end(); it++) { + collection.set[collection.n_sets].add_pair(*it, -1); + } - count++; - } - } // end for k + // Increment set counter + collection.n_sets++; - if(count > 0) collection.n_sets++; - } // end for j - - // - // Done - // - - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0; + } // end for mid return; } @@ -1415,12 +1403,10 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() { /////////////////////////////////////////////////////////////////////// void ModeFuzzyEngine::do_obs_merge_thresh() { - int j, k, x, y; - int n_obs_merge, intersection; - int count, first_k; + const char *method_name = "ModeFuzzyEngine::do_obs_merge_thresh() -> "; + int j, mid, oid; + int n_obs_merge; ShapeData obs_merge_mask, obs_merge_split; - ShapeData * obs_shape = (ShapeData *) 0; - ShapeData * obs_merge_shape = (ShapeData *) 0; do_obs_splitting(); @@ -1441,328 +1427,243 @@ void ModeFuzzyEngine::do_obs_merge_thresh() { obs_merge_split = split(obs_merge_mask, n_obs_merge); // - // Allocate space for all of the simple observation shapes and - // observation merge shapes + // Mapping of simple object id to merge object id values and vice-versa // - obs_shape = new ShapeData [n_obs]; - obs_merge_shape = new ShapeData [n_obs_merge]; + map > obj_to_mrg_map; + map > mrg_to_obj_map; + std::set empty_set; - if(!obs_shape || !obs_merge_shape) { + // Add map entries for each object + for(j=1; j<=n_obs; j++) obj_to_mrg_map[j] = empty_set; + for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set; - mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> " - << "memory allocation error\n\n"; - exit(1); - } + // Loop over all the data points + for(j=0; jdata.nxy(); j++) { + + oid = obs_split->data.buf()[j]; + mid = obs_merge_split.data.buf()[j]; + + // Update object id maps + if(oid > 0) obj_to_mrg_map[oid].insert(mid); + if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid); + + } // end for j // - // Select all of the simple observation shapes and - // simple merge shapes + // Each simple object should be fully contained in a merge object // - for(j=0; j= obs area - // - if(intersection >= obs_shape[k].area()) { - - collection.make_room(); + for(mid=1; mid<=n_obs_merge; mid++) { - // - // Keep track of the first embedded shape. You only want to - // create a composite if there are more than 1 shapes in it. - // - if(count == 0) first_k = k+1; + // Ignore merge objects containing a single simple object + if(mrg_to_obj_map[mid].size() <= 1) continue; - else if(count == 1) { - collection.set[collection.n_sets].add_pair(-1, first_k); - collection.set[collection.n_sets].add_pair(-1, k+1); - } - else { - collection.set[collection.n_sets].add_pair(-1, k+1); - } - - count++; - } - } // end for k - - if(count > 0) collection.n_sets++; - } // end for j + // Add all simple objects into one new set + collection.make_room(); + for(auto it = mrg_to_obj_map[mid].begin(); + it != mrg_to_obj_map[mid].end(); it++) { + collection.set[collection.n_sets].add_pair(-1, *it); + } - // - // Done - // + // Increment set counter + collection.n_sets++; - delete [] obs_shape; obs_shape = (ShapeData *) 0; - delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0; + } // end for mid return; } +/////////////////////////////////////////////////////////////////////// void ModeFuzzyEngine::do_fcst_merge_thresh(const ShapeData &merge_data) { - int j, k, x, y; - int n_fcst_merge, intersection; - int count, first_k; - //ShapeData fcst_merge_mask; + const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> "; + int j, mid, oid; + int n_fcst_merge; ShapeData fcst_merge_split; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * fcst_merge_shape = (ShapeData *) 0; do_fcst_splitting(); // - // Define the forecast merge field by applying the specified threshold - // to the convolved field - // - fcst_merge_split = merge_data; - n_fcst_merge = fcst_merge_split.n_objects(); - - // fcst_merge_mask = *fcst_conv; - - // // - // // Threshold the forecast merge field - // // - // fcst_merge_mask.threshold(conf_info.Fcst->merge_thresh); - - // // - // // Split up the forecast merge field - // // - // fcst_merge_split = split(fcst_merge_mask, n_fcst_merge); - - // - // Allocate space for all of the simple forecast shapes and - // forecast merge shapes + // Check dimensions // - fcst_shape = new ShapeData [n_fcst]; - fcst_merge_shape = new ShapeData [n_fcst_merge]; - - if(!fcst_shape || !fcst_merge_shape) { - - mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> " - << "memory allocation error\n\n"; + if(fcst_split->data.nx() != merge_data.data.nx() || + fcst_split->data.ny() != merge_data.data.ny()) { + mlog << Error << "\n" << method_name + << "grid dimensions do not match (" << fcst_split->data.nx() << ", " + << fcst_split->data.ny() << ") != (" << merge_data.data.nx() << ", " + << merge_data.data.ny() << ")!\n\n"; exit(1); } // - // Select all of the simple forecast shapes and - // forecast merge shapes + // Assume that the input merge data has already been split // - for(j=0; j > obj_to_mrg_map; + map > mrg_to_obj_map; + std::set empty_set; - // - // Calculate intersection area - // - intersection = 0; + // Add map entries for each object + for(j=1; j<=n_fcst; j++) obj_to_mrg_map[j] = empty_set; + for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set; - for(x=0; xdata.nxy(); j++) { - if(fcst_merge_shape[j].s_is_on(x, y) && - fcst_shape[k].s_is_on(x, y)) intersection++; - } - } + oid = fcst_split->data.buf()[j]; + mid = fcst_merge_split.data.buf()[j]; - // - // Add to set collection only if the fcst object is - // completely contained in the merge object. Meaning, - // intersection area >= fcst area - // + // Update object id maps + if(oid > 0) obj_to_mrg_map[oid].insert(mid); + if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid); - if(intersection >= fcst_shape[k].area()) { + } // end for j - collection.make_room(); + // + // Each simple object should be fully contained in a merge object + // + for(oid=1; oid<=n_fcst; oid++) { + if(obj_to_mrg_map[oid].size() != 1) { + mlog << Warning << "\n" << method_name + << "No forecast threshold merging done because simple object " + << oid << " is not fully contained within one merge object!\n\n"; + return; + } + } - // - // Keep track of the first embedded shape. You only want to - // create a composite if there are more than 1 shapes in it. - // - if(count == 0) first_k = k+1; + // + // Calculate the composite object sets + // + for(mid=1; mid<=n_fcst_merge; mid++) { - else if(count == 1) { - collection.set[collection.n_sets].add_pair(first_k, -1); - collection.set[collection.n_sets].add_pair(k+1, -1); - } - else { - collection.set[collection.n_sets].add_pair(k+1, -1); - } + // Ignore merge objects containing a single simple object + if(mrg_to_obj_map[mid].size() <= 1) continue; - count++; - } - } // end for k + // Add all simple objects into one new set + collection.make_room(); + for(auto it = mrg_to_obj_map[mid].begin(); + it != mrg_to_obj_map[mid].end(); it++) { + collection.set[collection.n_sets].add_pair(*it, -1); + } - if(count > 0) collection.n_sets++; - } // end for j + // Increment set counter + collection.n_sets++; - // - // Done - // + } // end for mid - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After merging of fcst field: " + << fcst_split->sdebug_examine() << "\n"; + } - mlog << Debug(4) << " After merging of fcst field: " << fcst_split->sdebug_examine() << "\n"; return; } +/////////////////////////////////////////////////////////////////////// + void ModeFuzzyEngine::do_obs_merge_thresh(const ShapeData &merge_data) { - int j, k, x, y; - int n_obs_merge, intersection; - int count, first_k; - //ShapeData obs_merge_mask; + const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> "; + int j, mid, oid; + int n_obs_merge; ShapeData obs_merge_split; - ShapeData * obs_shape = (ShapeData *) 0; - ShapeData * obs_merge_shape = (ShapeData *) 0; do_obs_splitting(); // - // Define the forecast merge field by applying the specified threshold - // to the convolved field - // - // obs_merge_mask = *obs_conv; - - // - // Threshold the forecast merge field + // Check dimensions // - // obs_merge_mask.threshold(conf_info.Obs->merge_thresh); + if(obs_split->data.nx() != merge_data.data.nx() || + obs_split->data.ny() != merge_data.data.ny()) { + mlog << Error << "\n" << method_name + << "grid dimensions do not match (" << obs_split->data.nx() << ", " + << obs_split->data.ny() << ") != (" << merge_data.data.nx() << ", " + << merge_data.data.ny() << ")!\n\n"; + exit(1); + } // - // Split up the forecast merge field + // Assume that the input merge data has already been split // - obs_merge_split = merge_data; //split(obs_merge_mask, n_obs_merge); + obs_merge_split = merge_data; n_obs_merge = obs_merge_split.n_objects(); // - // Allocate space for all of the simple observation shapes and - // observation merge shapes + // Mapping of simple object id to merge object id values and vice-versa // - obs_shape = new ShapeData [n_obs]; - obs_merge_shape = new ShapeData [n_obs_merge]; + map > obj_to_mrg_map; + map > mrg_to_obj_map; + std::set empty_set; - if(!obs_shape || !obs_merge_shape) { + // Add map entries for each object + for(j=1; j<=n_obs; j++) obj_to_mrg_map[j] = empty_set; + for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set; - mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> " - << "memory allocation error\n\n"; - exit(1); - } + // Loop over all the data points + for(j=0; jdata.nxy(); j++) { + + oid = obs_split->data.buf()[j]; + mid = obs_merge_split.data.buf()[j]; + + // Update object id maps + if(oid > 0) obj_to_mrg_map[oid].insert(mid); + if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid); + + } // end for j // - // Select all of the simple observation shapes and - // simple merge shapes + // Each simple object should be fully contained in a merge object // - for(j=0; j= obs area - // - if(intersection >= obs_shape[k].area()) { - - collection.make_room(); - - // - // Keep track of the first embedded shape. You only want to - // create a composite if there are more than 1 shapes in it. - // - if(count == 0) first_k = k+1; - - else if(count == 1) { - collection.set[collection.n_sets].add_pair(-1, first_k); - collection.set[collection.n_sets].add_pair(-1, k+1); - } - else { - collection.set[collection.n_sets].add_pair(-1, k+1); - } + // Add all simple objects into one new set + collection.make_room(); + for(auto it = mrg_to_obj_map[mid].begin(); + it != mrg_to_obj_map[mid].end(); it++) { + collection.set[collection.n_sets].add_pair(-1, *it); + } - count++; - } - } // end for k + // Increment set counter + collection.n_sets++; - if(count > 0) collection.n_sets++; - } // end for j + } // end for mid - // - // Done - // - - delete [] obs_shape; obs_shape = (ShapeData *) 0; - delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0; + if(mlog.verbosity_level() >= 4) { + mlog << Debug(4) << " After merging of obs field: " + << obs_split->sdebug_examine() << "\n"; + } - mlog << Debug(4) << " After merging of obs field: " << obs_split->sdebug_examine() << "\n"; return; } @@ -2114,31 +2015,22 @@ void ModeFuzzyEngine::do_obs_merge_engine(const char *default_config, void ModeFuzzyEngine::do_match_fcst_merge() { int j, k, n; InterestInfo junkinfo; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * obs_shape = (ShapeData *) 0; + ShapeData cur_shape; do_fcst_splitting(); do_obs_splitting(); clear_colors(); - fcst_shape = new ShapeData [n_fcst]; - obs_shape = new ShapeData [n_obs]; - - if(!fcst_shape || !obs_shape) { - mlog << Error << "\nModeFuzzyEngine::do_match_fcst_merge() -> " - << "memory allocation error\n\n"; - exit(1); - } - // // Do the single features // fcst_single.set_size(n_fcst); for(j=0; jvar_info->is_precipitation()); fcst_single[j].object_number = j+1; @@ -2147,8 +2039,9 @@ void ModeFuzzyEngine::do_match_fcst_merge() { obs_single.set_size(n_obs); for(j=0; jvar_info->is_precipitation()); obs_single[j].object_number = j+1; @@ -2275,9 +2168,6 @@ void ModeFuzzyEngine::do_match_fcst_merge() { // Done // - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] obs_shape; obs_shape = (ShapeData *) 0; - return; } @@ -2291,31 +2181,22 @@ void ModeFuzzyEngine::do_match_fcst_merge() { void ModeFuzzyEngine::do_match_only() { int j, k, n; InterestInfo junkinfo; - ShapeData * fcst_shape = (ShapeData *) 0; - ShapeData * obs_shape = (ShapeData *) 0; + ShapeData cur_shape; do_fcst_splitting(); do_obs_splitting(); clear_colors(); - fcst_shape = new ShapeData [n_fcst]; - obs_shape = new ShapeData [n_obs]; - - if(!fcst_shape || !obs_shape) { - mlog << Error << "\nModeFuzzyEngine::do_match_only() -> " - << "memory allocation error\n\n"; - exit(1); - } - // // Do the single features // fcst_single.set_size(n_fcst); for(j=0; jvar_info->is_precipitation()); fcst_single[j].object_number = j+1; @@ -2324,8 +2205,9 @@ void ModeFuzzyEngine::do_match_only() { obs_single.set_size(n_obs); for(j=0; jvar_info->is_precipitation()); obs_single[j].object_number = j+1; @@ -2446,9 +2328,6 @@ void ModeFuzzyEngine::do_match_only() { // Done // - delete [] fcst_shape; fcst_shape = (ShapeData *) 0; - delete [] obs_shape; obs_shape = (ShapeData *) 0; - return; } @@ -2537,8 +2416,7 @@ void ModeFuzzyEngine::do_obs_clus_splitting() { void ModeFuzzyEngine::do_cluster_features() { int j; - ShapeData * fcst_clus_shape = (ShapeData *) 0; - ShapeData * obs_clus_shape = (ShapeData *) 0; + ShapeData cur_shape; if(need_fcst_clus_split) do_fcst_clus_splitting(); if(need_obs_clus_split) do_obs_clus_splitting(); @@ -2548,16 +2426,6 @@ void ModeFuzzyEngine::do_cluster_features() { // n_clus = collection.n_sets; - fcst_clus_shape = new ShapeData [n_clus]; - obs_clus_shape = new ShapeData [n_clus]; - - if(!fcst_clus_shape || !obs_clus_shape) { - - mlog << Error << "\nModeFuzzyEngine::do_cluster_features() -> " - << "memory allocation error\n\n"; - exit(1); - } - // // Do the single features for clusters // @@ -2565,16 +2433,18 @@ void ModeFuzzyEngine::do_cluster_features() { obs_cluster.set_size(n_clus); for(j=0; jvar_info->is_precipitation()); + cur_shape = select(*fcst_clus_split, j+1); + fcst_cluster[j].set(*fcst_raw, *fcst_thresh, + *fcst_clus_split, cur_shape, + conf_info.inten_perc_value, + conf_info.Fcst->var_info->is_precipitation()); fcst_cluster[j].object_number = j+1; - obs_clus_shape[j] = select(*obs_clus_split, j+1); - obs_cluster[j].set(*obs_raw, *obs_thresh, obs_clus_shape[j], - conf_info.inten_perc_value, - conf_info.Obs->var_info->is_precipitation()); + cur_shape = select(*obs_clus_split, j+1); + obs_cluster[j].set(*obs_raw, *obs_thresh, + *obs_clus_split, cur_shape, + conf_info.inten_perc_value, + conf_info.Obs->var_info->is_precipitation()); obs_cluster[j].object_number = j+1; } @@ -2585,7 +2455,7 @@ void ModeFuzzyEngine::do_cluster_features() { for(j=0; jcentroid(centroid_x, centroid_y); + mask_sd.centroid(centroid_x, centroid_y); // // Axis angle // - axis_ang = Mask->angle_degrees(); + axis_ang = mask_sd.angle_degrees(); // // Length and Width // - Mask->calc_length_width(length, width); + mask_sd.calc_length_width(length, width); // // Aspect ratio @@ -228,28 +228,28 @@ void SingleFeature::set(const ShapeData &raw, const ShapeData &thresh, // // Object area // - area = Mask->area(); + area = mask_sd.area(); // // Object threshold area: the area of the raw field inside the mask // area that meets the threshold criteria // - area_thresh = (double) ShapeData_intersection(*Thresh, *Mask); + area_thresh = (double) ShapeData_intersection(*Thresh, mask_sd); // // Curvature, Curvature_x, Curvature_y // - curvature = Mask->curvature(curvature_x, curvature_y); + curvature = mask_sd.curvature(curvature_x, curvature_y); // // Complexity // - complexity = Mask->complexity(); + complexity = mask_sd.complexity(); // // Compute the Intensity Percentiles // - get_percentiles(intensity_ptile, raw, mask, perc, precip_flag); + get_percentiles(intensity_ptile, raw_sd, mask_sd, perc, precip_flag); // // User Percentile @@ -259,17 +259,17 @@ void SingleFeature::set(const ShapeData &raw, const ShapeData &thresh, // // Convex hull // - convex_hull = Mask->convex_hull(); + convex_hull = mask_sd.convex_hull(); // // Boundary: // Split the mask field and store the boundary for each object. // - split_wd = split(mask, n_bdy); + cur_split_sd = split(mask_sd, n_bdy); boundary = new Polyline [n_bdy]; for(i=0; iarea)/(Fcst->area), (Fcst->area)/(Obs->area) ); - // // Intersection, union, and symmetric diff areas // intersection_area = union_area = 0.0; symmetric_diff = 0.0; - for(x=0; x<(Fcst->Mask->data.nx()); ++x) { - for(y=0; y<(Fcst->Mask->data.ny()); ++y) { + int nxy = Fcst->Split->data.nxy(); + for(i=0; iMask->s_is_on(x, y); - obs_on = Obs->Mask->s_is_on(x, y); + fcst_on = (nint(Fcst->Split->data.data()[i]) == Fcst->object_number ? 1 : 0); + obs_on = (nint( Obs->Split->data.data()[i]) == Obs->object_number ? 1 : 0); - if(fcst_on && obs_on) intersection_area++; - if(fcst_on || obs_on) union_area++; - if((fcst_on && !obs_on) || - (!fcst_on && obs_on)) symmetric_diff++; - } + if( fcst_on && obs_on) intersection_area++; + if( fcst_on || obs_on) union_area++; + if( ( fcst_on && !obs_on) || + (!fcst_on && obs_on) ) symmetric_diff++; } // diff --git a/src/libcode/vx_shapedata/interest.h b/src/libcode/vx_shapedata/interest.h index de8832ed2b..6fd6f9cc2d 100644 --- a/src/libcode/vx_shapedata/interest.h +++ b/src/libcode/vx_shapedata/interest.h @@ -49,13 +49,13 @@ class SingleFeature { void clear(); - void set(const ShapeData &raw, const ShapeData &thresh, - const ShapeData &mask, const int perc, - const bool precip_flag); + void set(const ShapeData &raw, const ShapeData &thresh, + const ShapeData &split, const ShapeData &mask, + const int perc, const bool precip_flag); const ShapeData * Raw; // NOT allocated, so DON'T delete! const ShapeData * Thresh; // NOT allocated, so DON'T delete! - const ShapeData * Mask; // NOT allocated, so DON'T delete! + const ShapeData * Split; // NOT allocated, so DON'T delete! int object_number; diff --git a/src/libcode/vx_shapedata/shapedata.cc b/src/libcode/vx_shapedata/shapedata.cc index 9051d00819..88b63e195b 100644 --- a/src/libcode/vx_shapedata/shapedata.cc +++ b/src/libcode/vx_shapedata/shapedata.cc @@ -15,8 +15,9 @@ // // Mod# Date Name Description // ---- ---- ---- ----------- -// 000 11-05-31 Halley Gotway Adapated from wrfdata.cc. -// 001 14-05-29 Halley Gotway Add ShapeData::n_objects() +// 000 05/31/11 Halley Gotway Adapated from wrfdata.cc +// 001 05/29/14 Halley Gotway Add ShapeData::n_objects() +// 002 11/02/23 Halley Gotway MET #2724 add OpenMP to convolution // /////////////////////////////////////////////////////////////////////////////// @@ -513,160 +514,117 @@ double ShapeData::get_attr(const ConcatString &attr_name, /////////////////////////////////////////////////////////////////////////////// +void ShapeData::conv_filter_circ(int diameter, double vld_thresh) { + const char *method_name = "ShapeData::conv_filter_circ() -> "; + GridPoint *gp = nullptr; + int x, y, count, n_vld; + double v, v_sum; + DataPlane conv_dp; -void ShapeData::conv_filter_circ(int diameter, double vld_thresh) - -{ - -int x, y, xx, yy, u, v; -int dn, fn, nn; -int vpr, upr; -int count, bd_count; -double center, cur, sum; -double dx, dy, dist; -double vld_ratio; -const int nx = data.nx(); -const int ny = data.ny(); -bool * f = (bool *) 0; -bool center_bad = false; -DataPlane in_data = data; -const bool vld_thresh_one = is_eq(vld_thresh, 1.0); - - -if ( (diameter%2 == 0) || (diameter < 3) ) { - - mlog << Error << "\nShapeData::conv_filter_circ() -> " - << "diameter must be odd and >= 3 ... diameter = " - << diameter << "\n\n"; - - exit(1); - -} - -const int radius = (diameter - 1)/2; - -const vector * in = &(in_data.Data); - vector * out = &(data.Data); - -f = new bool [diameter*diameter]; - - // - // set up the filter - // - -for (y=0; y= 3 ... diameter = " + << diameter << "\n\n"; + exit(1); } -} - - // - // do the convolution - // - -dn = -1; - -for(y=0; y= ny) ) continue; - - vpr = v + radius; - - for(u=-radius; u<=radius; ++u) { - - xx = x + u; - - if ( (xx < 0) || (xx >= nx) ) continue; - - upr = u + radius; - - fn = STANDARD_XY_YO_N(diameter, upr, vpr); - - if ( !(f[fn]) ) continue; - - nn = STANDARD_XY_YO_N(nx, xx, yy) ; - - cur = (*in)[nn]; - - if( ::is_bad_data(cur) ) { bd_count++; continue; } - - sum += cur; - - count++; - - } // for v - - } // for u - - // - // If the center of the convolution contains bad data and the ratio - // of bad data in the convolution area is too high, set the convoled - // value to bad data. - // - - if ( count == 0 ) sum = bad_data_double; - else { +#pragma omp parallel default(none) \ + shared(mlog, data, conv_dp, diameter, vld_thresh) \ + private(x, y, count, n_vld, v, v_sum, gp) + { + + // Build the grid template with shape circle and wrap_lon false + GridTemplateFactory gtf; + GridTemplate* gt = gtf.buildGT(GridTemplateFactory::GridTemplate_Circle, + diameter, false); + +#pragma omp single + { + // Initialize the convolved field to bad data + conv_dp = data; + conv_dp.set_constant(bad_data_double); + } - vld_ratio = ((double) count)/(bd_count + count); + // Compute the convolved values +#pragma omp for schedule (static) + for(x=0; xgetFirstInGrid(x, y, data.nx(), data.ny()); + gp != nullptr; + gp = gt->getNextInGrid()) { + v = data.get(gp->x, gp->y); + count += 1; + if(::is_bad_data(v)) continue; + n_vld += 1; + v_sum += v; + } + } + // Subtract off the bottom edge, shift up, and add the top + else { + + // Subtract points from the the bottom edge + for(gp = gt->getFirstInBotEdge(); + gp != nullptr; + gp = gt->getNextInBotEdge()) { + v = data.get(gp->x, gp->y); + count -= 1; + if(::is_bad_data(v)) continue; + n_vld -= 1; + v_sum -= v; + } + + // Increment Y + gt->incBaseY(1); + + // Add points from the the top edge + for(gp = gt->getFirstInTopEdge(); + gp != nullptr; + gp = gt->getNextInTopEdge()) { + v = data.get(gp->x, gp->y); + count += 1; + if(::is_bad_data(v)) continue; + n_vld += 1; + v_sum += v; + } + } - if ( center_bad && (vld_ratio < vld_thresh) ) sum = bad_data_double; - else sum /= count; + // If the center of the convolution contains bad data and the ratio + // of bad data in the convolution area is too high, set the convoled + // value to bad data. + if(count == 0 || n_vld == 0) v = bad_data_double; + else if(::is_bad_data(data.get(x, y)) && + (double)(n_vld)/count < vld_thresh) v = bad_data_double; + else v = (double) v_sum/n_vld; + conv_dp.set(v, x, y); - } + } // end for y - (*out)[dn] = sum; + // Increment X + if(x < (data.nx() - 1)) gt->incBaseX(1); - } // for y + } // end for x -} // for x + delete gt; - // - // done - // + } // End of omp parallel -if ( f ) { delete [] f; f = (bool *) 0; } - -return; + // Save the result + data = conv_dp; + return; } - //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/ensemble_stat/ensemble_stat.cc b/src/tools/core/ensemble_stat/ensemble_stat.cc index 62885a5f5e..b778e97f2f 100644 --- a/src/tools/core/ensemble_stat/ensemble_stat.cc +++ b/src/tools/core/ensemble_stat/ensemble_stat.cc @@ -99,8 +99,6 @@ using namespace netCDF; #include "nc_obs_util.h" #include "nc_point_obs_in.h" -#include "handle_openmp.h" - #ifdef WITH_PYTHON #include "data2d_nc_met.h" #include "pointdata_python.h" @@ -177,9 +175,6 @@ static void set_compress(const StringArray &); int met_main(int argc, char *argv[]) { - // Set up OpenMP (if enabled) - init_openmp(); - // Process the command line arguments process_command_line(argc, argv); diff --git a/src/tools/core/mode/mode.cc b/src/tools/core/mode/mode.cc index 88e18d0b40..a9277a0884 100644 --- a/src/tools/core/mode/mode.cc +++ b/src/tools/core/mode/mode.cc @@ -56,6 +56,7 @@ // 019 04/01/19 Fillmore Add FCST and OBS units. // 020 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main // 021 06/09/23 Albo Major changes for multivariate mode +// 022 11/02/23 Halley Gotway MET #2724 add OpenMP to convolution // //////////////////////////////////////////////////////////////////////// @@ -73,6 +74,7 @@ using namespace std; #include #include "main.h" +#include "handle_openmp.h" #include "string_array.h" #include "mode_usage.h" #include "mode_frontend.h" @@ -131,6 +133,9 @@ int met_main(int argc, char * argv []) string s; const char * user_config_filename = 0; + // Set up OpenMP (if enabled) + init_openmp(); + for (j=0,n=0; j