Skip to content

Commit

Permalink
Per #2724, reimplement ShapeData::conv_filter_circ() using OpenMP to …
Browse files Browse the repository at this point in the history 
…make it more efficient.
  • Loading branch information
@JohnHalleyGotway
JohnHalleyGotway committed Nov 2, 2023
1 parent 3ea3a8a commit ba5cc6e
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 142 deletions.
1 change: 1 addition & 0 deletions src/basic/vx_util/vx_util.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
235 changes: 93 additions & 142 deletions src/libcode/vx_shapedata/shapedata.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
///////////////////////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -513,160 +514,110 @@ 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, 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<double> * in = &(in_data.Data);
vector<double> * out = &(data.Data);

f = new bool [diameter*diameter];

//
// set up the filter
//

for (y=0; y<diameter; ++y) {

dy = y - radius;

for (x=0; x<diameter; ++x) {

dx = x - radius;

dist = sqrt( dx*dx + dy*dy );

fn = STANDARD_XY_YO_N(diameter, x, y) ;

f[fn] = (dist <= radius);

// Check the diameter
if(diameter%2 == 0 || diameter < 3) {
mlog << Error << "\n" << method_name
<< "diameter must be odd and >= 3 ... diameter = "
<< diameter << "\n\n";
exit(1);
}

}

//
// do the convolution
//

dn = -1;

for(y=0; y<ny; y++) {

for(x=0; x<nx; x++) {

++dn;

//
// If the bad data threshold is set to zero and the center of the
// convolution radius contains bad data, set the convolved value to
// bad data and continue.
//

center = (*in)[dn];

center_bad = ::is_bad_data(center);

if ( center_bad && vld_thresh_one ) { (*out)[dn] = bad_data_double; continue; }

sum = 0.0;
count = 0;
bd_count = 0;

for (v=-radius; v<=radius; ++v) {

yy = y + v;

if ( (yy < 0) || (yy >= 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) \
shared(vld_thresh, bad_data_double) \
private(x, y, 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; x<data.nx(); x++) {
for(y=0; y<data.ny(); y++) {

// For a new column, reset the grid template and counts
if(y == 0) {

// Initialize count and sum
n_vld = 0;
v_sum = 0.0;

// Sum all the points
for(gp = gt->getFirstInGrid(x, y, data.nx(), data.ny());
gp != nullptr;
gp = gt->getNextInGrid()) {
v = data.get(gp->x, gp->y);
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);
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);
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;
// Check for enough valid data and compute the mean value
if((double)(n_vld)/gt->size() >= vld_thresh && n_vld != 0) {
conv_dp.set((double) v_sum/n_vld, 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;
}


////////////////////////////////////////////////////////////////////////


Expand Down
5 changes: 5 additions & 0 deletions src/tools/core/mode/mode.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
//
////////////////////////////////////////////////////////////////////////

Expand All @@ -73,6 +74,7 @@ using namespace std;
#include <sys/types.h>

#include "main.h"
#include "handle_openmp.h"
#include "string_array.h"
#include "mode_usage.h"
#include "mode_frontend.h"
Expand Down Expand Up @@ -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<argc; ++j) {

//
Expand Down

0 comments on commit ba5cc6e

Please sign in to comment.