Skip to content

Commit

Permalink
Per #1451, instead of computing the climo crps on the fly, compute an…
Browse files Browse the repository at this point in the history
…d store it separately for each point.
  • Loading branch information
JohnHalleyGotway committed Feb 15, 2021
1 parent e5a08ea commit e39201d
Showing 1 changed file with 20 additions and 41 deletions.
61 changes: 20 additions & 41 deletions met/src/libcode/vx_statistics/ens_stats.cc
Expand Up @@ -228,10 +228,26 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
// Store the number of ensemble members
n_ens = pd.n_ens;

// Get the average empirical and Gaussian CRPS and
// components values
crps_emp = pd.crps_emp_na.wmean(pd.wgt_na);
crps_gaus = pd.crps_gaus_na.wmean(pd.wgt_na);
// Compute empirical CRPS scores
crps_emp = pd.crps_emp_na.wmean(pd.wgt_na);
crpscl_emp = pd.crpscl_emp_na.wmean(pd.wgt_na);
crpss_emp = (is_bad_data(crps_emp) ||
is_bad_data(crpscl_emp) ||
is_eq(crpscl_emp, 0.0) ?
bad_data_double :
(crpscl_emp - crps_emp) / crpscl_emp);

// Compute Gaussian CRPS scores
crps_gaus = pd.crps_gaus_na.wmean(pd.wgt_na);
crpscl_gaus = pd.crpscl_gaus_na.wmean(pd.wgt_na);
crpss_gaus = (is_bad_data(crps_gaus) ||
is_bad_data(crpscl_gaus) ||
is_eq(crpscl_gaus, 0.0) ?
bad_data_double :
(crpscl_gaus - crps_gaus) / crpscl_gaus);

// Compute the average IGN value
ign = pd.ign_na.wmean(pd.wgt_na);

// Get the sum of the weights
for(i=0, n_pair=0, w_sum=0.0; i<pd.wgt_na.n(); i++) {
Expand All @@ -241,43 +257,6 @@ void ECNTInfo::set(const PairDataEnsemble &pd) {
}
}

// Check for bad data
if(is_bad_data(crps_emp) ||
is_bad_data(crps_gaus) ||
pd.cmn_na.n() != pd.o_na.n() ||
pd.cmn_na.n() == 0 ||
pd.cmn_na.has(bad_data_double)) {
crpss_emp = crpss_gaus = bad_data_double;
}
else {

// Compute the climatological CRPS
ffbar = oobar = fobar = 0.0;
for(i=0; i<pd.n_obs; i++) {

if(pd.skip_ba[i]) continue;

// Track running sums
w = pd.wgt_na[i]/w_sum;
ffbar += w * pd.cmn_na[i] * pd.cmn_na[i];
oobar += w * pd.o_na[i] * pd.o_na[i];
fobar += w * pd.cmn_na[i] * pd.o_na[i];
}

// TODO: MET #1451 correct this computation
crpscl_emp = ffbar + oobar - 2.0*fobar;
crpscl_gaus = ffbar + oobar - 2.0*fobar;

// Compute CRPS skill scores
crpss_emp = (is_eq(crpscl_emp, 0.0) ? bad_data_double :
(crpscl_emp - crps_emp) /crpscl_emp);
crpss_gaus = (is_eq(crpscl_gaus, 0.0) ? bad_data_double :
(crpscl_gaus - crps_gaus) /crpscl_gaus);
}

// Compute the average IGN value
ign = pd.ign_na.wmean(pd.wgt_na);

// Compute ME and RMSE values
fbar = obar = ffbar = oobar = fobar = 0.0;
for(i=0; i<pd.n_obs; i++) {
Expand Down

0 comments on commit e39201d

Please sign in to comment.