Commit 3c9dd882 authored by Thibault Hallouin's avatar Thibault Hallouin
Browse files

centralise computation of dimension sizes

These calculations were redundantly calculated for each metric
separately when there was no reason to because the dimensions do
not change overtime.
Showing with 24 additions and 46 deletions
+24 -46
...@@ -16,6 +16,12 @@ namespace evalhyd ...@@ -16,6 +16,12 @@ namespace evalhyd
const xt::xtensor<double, 1>& q_thr; const xt::xtensor<double, 1>& q_thr;
xt::xtensor<bool, 2> t_msk; xt::xtensor<bool, 2> t_msk;
// members for dimensions
std::size_t n;
std::size_t n_msk;
std::size_t n_mbr;
std::size_t n_thr;
// members for computational elements // members for computational elements
xt::xtensor<double, 2> o_k; xt::xtensor<double, 2> o_k;
xt::xtensor<double, 2> bar_o; xt::xtensor<double, 2> bar_o;
...@@ -30,8 +36,16 @@ namespace evalhyd ...@@ -30,8 +36,16 @@ namespace evalhyd
const xt::xtensor<bool, 2>& msk) : const xt::xtensor<bool, 2>& msk) :
q_obs{obs}, q_prd{prd}, q_thr{thr}, t_msk(msk) q_obs{obs}, q_prd{prd}, q_thr{thr}, t_msk(msk)
{ {
if (msk.size() < 1) // initialise a mask if none provided
t_msk = xt::ones<bool>(obs.shape()); // (corresponding to no temporal subset)
if (msk.size() < 1)
t_msk = xt::ones<bool>(obs.shape());
// determine size for recurring dimensions
n = obs.shape(1);
n_msk = t_msk.shape(0);
n_mbr = prd.shape(0);
n_thr = thr.size();
}; };
// members for intermediate evaluation metrics // members for intermediate evaluation metrics
......
...@@ -48,10 +48,6 @@ namespace evalhyd ...@@ -48,10 +48,6 @@ namespace evalhyd
// shape: (subsets, thresholds, 1) // shape: (subsets, thresholds, 1)
void Evaluator::calc_BS() void Evaluator::calc_BS()
{ {
// define some dimensions
std::size_t n_msk = t_msk.shape(0);
std::size_t n_thr = bs.shape(0);
// initialise output variable // initialise output variable
// shape: (subsets, thresholds, 1) // shape: (subsets, thresholds, 1)
BS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}}); BS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}});
...@@ -104,18 +100,12 @@ namespace evalhyd ...@@ -104,18 +100,12 @@ namespace evalhyd
// shape: (bins,) // shape: (bins,)
xt::xtensor<double, 1> y_i; xt::xtensor<double, 1> y_i;
// define some dimensions
std::size_t n_thr = o_k.shape(0);
std::size_t n_mbr = q_prd.shape(0);
std::size_t n_msk = t_msk.shape(0);
std::size_t n_cmp = 3;
// compute range of forecast values $y_i$ // compute range of forecast values $y_i$
y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr; y_i = xt::arange<double>(double(n_mbr + 1)) / n_mbr;
// initialise output variable // initialise output variable
// shape: (subsets, thresholds, components) // shape: (subsets, thresholds, components)
BS_CRD = xt::zeros<double>({n_msk, n_thr, n_cmp}); BS_CRD = xt::zeros<double>({n_msk, n_thr, std::size_t {3}});
// compute variable one mask at a time to minimise memory imprint // compute variable one mask at a time to minimise memory imprint
for (int m = 0; m < n_msk; m++) for (int m = 0; m < n_msk; m++)
...@@ -211,17 +201,12 @@ namespace evalhyd ...@@ -211,17 +201,12 @@ namespace evalhyd
// shape: (bins,) // shape: (bins,)
xt::xtensor<double, 1> o_j; xt::xtensor<double, 1> o_j;
// define some dimensions
std::size_t n_thr = o_k.shape(0);
std::size_t n_msk = t_msk.shape(0);
std::size_t n_cmp = 3;
// set the range of observed values $o_j$ // set the range of observed values $o_j$
o_j = {0., 1.}; o_j = {0., 1.};
// declare and initialise output variable // declare and initialise output variable
// shape: (subsets, thresholds, components) // shape: (subsets, thresholds, components)
BS_LBD = xt::zeros<double>({n_msk, n_thr, n_cmp}); BS_LBD = xt::zeros<double>({n_msk, n_thr, std::size_t {3}});
// compute variable one mask at a time to minimise memory imprint // compute variable one mask at a time to minimise memory imprint
for (int m = 0; m < n_msk; m++) for (int m = 0; m < n_msk; m++)
...@@ -319,10 +304,6 @@ namespace evalhyd ...@@ -319,10 +304,6 @@ namespace evalhyd
// shape: (subsets, thresholds, 1) // shape: (subsets, thresholds, 1)
void Evaluator::calc_BSS() void Evaluator::calc_BSS()
{ {
// define some dimensions
std::size_t n_thr = o_k.shape(0);
std::size_t n_msk = t_msk.shape(0);
// declare and initialise output variable // declare and initialise output variable
// shape: (subsets, thresholds, components) // shape: (subsets, thresholds, components)
BSS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}}); BSS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}});
......
...@@ -70,13 +70,9 @@ namespace evalhyd ...@@ -70,13 +70,9 @@ namespace evalhyd
// calculate how many members have exceeded threshold(s) // calculate how many members have exceeded threshold(s)
xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1); xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
// determine correspondence between number of exceeding members // determine probability of threshold(s) exceedance
// and forecast probabilities of exceedance
// /!\ probability calculation dividing by n (the number of // /!\ probability calculation dividing by n (the number of
// members), not n+1 (the number of ranks) like in other metrics // members), not n+1 (the number of ranks) like in other metrics
std::size_t n_mbr = q_prd.shape(0);
// determine probability of threshold(s) exceedance
y_k = n_frc / n_mbr; y_k = n_frc / n_mbr;
} }
......
...@@ -24,13 +24,10 @@ namespace evalhyd ...@@ -24,13 +24,10 @@ namespace evalhyd
// shape: (quantiles, time) // shape: (quantiles, time)
void Evaluator::calc_qs() void Evaluator::calc_qs()
{ {
// define some dimensions
std::size_t n_qnt = q_qnt.shape(0);
// compute the quantile order $alpha$ // compute the quantile order $alpha$
xt::xtensor<double, 1> alpha = xt::xtensor<double, 1> alpha =
xt::arange<double>(1., double(n_qnt + 1)) xt::arange<double>(1., double(n_mbr + 1))
/ double(n_qnt + 1); / double(n_mbr + 1);
// calculate the difference // calculate the difference
xt::xtensor<double, 2> diff = q_qnt - q_obs; xt::xtensor<double, 2> diff = q_qnt - q_obs;
...@@ -56,13 +53,9 @@ namespace evalhyd ...@@ -56,13 +53,9 @@ namespace evalhyd
// shape: (subsets, quantiles, 1) // shape: (subsets, quantiles, 1)
void Evaluator::calc_QS() void Evaluator::calc_QS()
{ {
// define some dimensions
std::size_t n_msk = t_msk.shape(0);
std::size_t n_qnt = qs.shape(0);
// initialise output variable // initialise output variable
// shape: (subsets, thresholds, 1) // shape: (subsets, quantiles, 1)
QS = xt::zeros<double>({n_msk, n_qnt, std::size_t {1}}); QS = xt::zeros<double>({n_msk, n_mbr, std::size_t {1}});
// compute variable one mask at a time to minimise memory imprint // compute variable one mask at a time to minimise memory imprint
for (int m = 0; m < n_msk; m++) for (int m = 0; m < n_msk; m++)
...@@ -94,13 +87,10 @@ namespace evalhyd ...@@ -94,13 +87,10 @@ namespace evalhyd
// shape: (1, time) // shape: (1, time)
void Evaluator::calc_crps() void Evaluator::calc_crps()
{ {
// define some dimensions
std::size_t n_qnt = q_qnt.shape(0);
// integrate with trapezoidal rule // integrate with trapezoidal rule
crps = xt::view( crps = xt::view(
// xt::trapz(y, dx=1/(n+1), axis=0) // xt::trapz(y, dx=1/(n+1), axis=0)
xt::trapz(qs, 1./(double(n_qnt) + .1), 0), xt::trapz(qs, 1./(double(n_mbr) + .1), 0),
xt::newaxis(), xt::all() xt::newaxis(), xt::all()
); );
} }
...@@ -119,9 +109,6 @@ namespace evalhyd ...@@ -119,9 +109,6 @@ namespace evalhyd
// shape: (subsets, 1, 1) // shape: (subsets, 1, 1)
void Evaluator::calc_CRPS() void Evaluator::calc_CRPS()
{ {
// define some dimensions
std::size_t n_msk = t_msk.shape(0);
// initialise output variable // initialise output variable
// shape: (subsets, thresholds, 1) // shape: (subsets, thresholds, 1)
CRPS = xt::zeros<double>({n_msk, std::size_t {1}, std::size_t {1}}); CRPS = xt::zeros<double>({n_msk, std::size_t {1}, std::size_t {1}});
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment