diff --git a/include/evalhyd/detail/determinist/elements.hpp b/include/evalhyd/detail/determinist/elements.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8c8fe1b11f0da79a4771c38050aebe40f6e64715
--- /dev/null
+++ b/include/evalhyd/detail/determinist/elements.hpp
@@ -0,0 +1,460 @@
+#ifndef EVALHYD_DETERMINIST_ELEMENTS_HPP
+#define EVALHYD_DETERMINIST_ELEMENTS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+
+#include "../maths.hpp"
+
+
+namespace evalhyd
+{
+ namespace determinist
+ {
+ namespace elements
+ {
+ // Compute the mean of the observations.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (1, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Mean observed streamflow.
+ // shape: (subsets, samples, series, 1)
+ template <class D2>
+ inline xt::xtensor<double, 4> calc_mean_obs(
+ const D2& q_obs,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 4> mean_obs =
+ xt::zeros<double>({n_msk, n_exp, n_srs, std::size_t {1}});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(mean_obs, m, e) =
+ xt::nanmean(obs, -1, xt::keep_dims);
+ }
+ }
+
+ return mean_obs;
+ }
+
+ // Compute the mean of the predictions.
+ //
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Mean predicted streamflow.
+ // shape: (subsets, samples, series, 1)
+ template <class D2>
+ inline xt::xtensor<double, 4> calc_mean_prd(
+ const D2& q_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 4> mean_prd =
+ xt::zeros<double>({n_msk, n_exp, n_srs, std::size_t {1}});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // apply the bootstrap sampling
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ xt::view(mean_prd, m, e) =
+ xt::nanmean(prd, -1, xt::keep_dims);
+ }
+ }
+
+ return mean_prd;
+ }
+
+ // Compute the quadratic error between observations and predictions.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (1, time)
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \return
+ // Quadratic errors between observations and predictions.
+ // shape: (series, time)
+ template <class D2>
+ inline xt::xtensor<double, 2> calc_quad_err(
+ const D2& q_obs,
+ const D2& q_prd
+ )
+ {
+ return xt::square(q_obs - q_prd);
+ }
+
+ // Compute the quadratic error between observations and mean observation.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (series, time)
+ // \param mean_obs
+ // Mean observed streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_tim
+ // Number of time steps.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Quadratic errors between observations and mean observation.
+ // shape: (subsets, samples, series, time)
+ template <class D2>
+ inline xt::xtensor<double, 4> calc_quad_obs(
+ const D2& q_obs,
+ const xt::xtensor<double, 4>& mean_obs,
+ const xt::xtensor<bool, 3>& t_msk,
+ std::size_t n_srs,
+ std::size_t n_tim,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 4> quad_obs =
+ xt::zeros<double>({n_msk, n_exp, n_srs, n_tim});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ xt::view(quad_obs, m, e) = xt::square(
+ obs_masked - xt::view(mean_obs, m, e)
+ );
+ }
+ }
+
+ return quad_obs;
+ }
+
+ // Compute the quadratic error between predictions and mean prediction.
+ //
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \param mean_prd
+ // Mean predicted streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_tim
+ // Number of time steps.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Quadratic errors between predictions and mean prediction.
+ // shape: (subsets, samples, series, time)
+ template <class D2>
+ inline xt::xtensor<double, 4> calc_quad_prd(
+ const D2& q_prd,
+ const xt::xtensor<double, 4>& mean_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ std::size_t n_srs,
+ std::size_t n_tim,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 4> quad_prd =
+ xt::zeros<double>({n_msk, n_exp, n_srs, n_tim});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ xt::view(quad_prd, m, e) = xt::square(
+ prd_masked - xt::view(mean_prd, m, e)
+ );
+ }
+ }
+
+ return quad_prd;
+ }
+
+ // Compute the Pearson correlation coefficient.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (series, time)
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \param mean_obs
+ // Mean observed streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param mean_prd
+ // Mean predicted streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param quad_obs
+ // Quadratic errors between observations and mean observation.
+ // shape: (subsets, samples, series, time)
+ // \param quad_prd
+ // Quadratic errors between predictions and mean prediction.
+ // shape: (subsets, samples, series, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Pearson correlation coefficients.
+ // shape: (subsets, samples, series)
+ template <class D2>
+ inline xt::xtensor<double, 3> calc_r_pearson(
+ const D2& q_obs,
+ const D2& q_prd,
+ const xt::xtensor<double, 4>& mean_obs,
+ const xt::xtensor<double, 4>& mean_prd,
+ const xt::xtensor<double, 4>& quad_obs,
+ const xt::xtensor<double, 4>& quad_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // calculate error in timing and dynamics $r_{pearson}$
+ // (Pearson's correlation coefficient)
+ xt::xtensor<double, 3> r_pearson =
+ xt::zeros<double>({n_msk, n_exp, n_srs});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ auto r_num = xt::nansum(
+ (prd - xt::view(mean_prd, m, e))
+ * (obs - xt::view(mean_obs, m, e)),
+ -1
+ );
+
+ auto prd2 = xt::view(quad_prd, m, e, xt::all(), b_exp[e]);
+ auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
+ auto r_den = xt::sqrt(
+ xt::nansum(prd2, -1) * xt::nansum(obs2, -1)
+ );
+
+ xt::view(r_pearson, m, e) = r_num / r_den;
+ }
+ }
+
+ return r_pearson;
+ }
+
+ // Compute alpha.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (series, time)
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \param mean_obs
+ // Mean observed streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param mean_prd
+ // Mean predicted streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Alphas, ratios of standard deviations.
+ // shape: (subsets, samples, series)
+ template <class D2>
+ inline xt::xtensor<double, 3> calc_alpha(
+ const D2& q_obs,
+ const D2& q_prd,
+ const xt::xtensor<double, 4>& mean_obs,
+ const xt::xtensor<double, 4>& mean_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // calculate error in spread of flow $alpha$
+ xt::xtensor<double, 3> alpha =
+ xt::zeros<double>({n_msk, n_exp, n_srs});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(alpha, m, e) =
+ maths::nanstd(prd, xt::view(mean_prd, m, e))
+ / maths::nanstd(obs, xt::view(mean_obs, m, e));
+ }
+ }
+
+ return alpha;
+ }
+
+ // Compute the bias.
+ //
+ // \param q_obs
+ // Streamflow observations.
+ // shape: (series, time)
+ // \param q_prd
+ // Streamflow predictions.
+ // shape: (series, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Biases.
+ // shape: (subsets, samples, series)
+ template <class D2>
+ inline xt::xtensor<double, 3> calc_bias(
+ const D2& q_obs,
+ const D2& q_prd,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // calculate $bias$
+ xt::xtensor<double, 3> bias =
+ xt::zeros<double>({n_msk, n_exp, n_srs});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
+ auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
+ auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
+ xt::view(bias, m, e) =
+ xt::nansum(prd, -1) / xt::nansum(obs, -1);
+ }
+ }
+
+ return bias;
+ }
+ }
+ }
+}
+
+#endif //EVALHYD_DETERMINIST_ELEMENTS_HPP
diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp
index 74b92d0348d3ca49ae0668c4ea6365ec7c904a59..8e494af48fcd524b2e43e7bdfeb4798e2f6d47f5 100644
--- a/include/evalhyd/detail/determinist/evaluator.hpp
+++ b/include/evalhyd/detail/determinist/evaluator.hpp
@@ -3,10 +3,13 @@
#include <vector>
+#include <xtl/xoptional.hpp>
#include <xtensor/xexpression.hpp>
#include <xtensor/xtensor.hpp>
-#include "../maths.hpp"
+#include "elements.hpp"
+#include "metrics.hpp"
+
namespace evalhyd
{
@@ -27,492 +30,195 @@ namespace evalhyd
std::size_t n_msk;
std::size_t n_srs;
std::size_t n_exp;
- std::vector<std::size_t> inner_dims;
- std::vector<std::size_t> inter_dims;
- std::vector<std::size_t> mean_dims;
- std::vector<std::size_t> final_dims;
// members for computational elements
- xt::xtensor<double, 4> mean_obs;
- xt::xtensor<double, 4> mean_prd;
- xt::xtensor<double, 2> quad_err;
- xt::xtensor<double, 4> quad_obs;
- xt::xtensor<double, 4> quad_prd;
- xt::xtensor<double, 3> r_pearson;
- xt::xtensor<double, 3> alpha;
- xt::xtensor<double, 3> bias;
-
- public:
- // constructor method
- Evaluator(const D2& obs,
- const D2& prd,
- const B2& msk,
- const std::vector<xt::xkeep_slice<int>>& exp) :
- q_obs{obs}, q_prd{prd}, b_exp{exp}
- {
- // determine size for recurring dimensions
- n_tim = q_prd.shape(1);
- n_srs = q_prd.shape(0);
- n_msk = msk.shape(0);
- n_exp = b_exp.size();
+ xtl::xoptional<xt::xtensor<double, 4>, bool> mean_obs;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> mean_prd;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> quad_err;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> quad_obs;
+ xtl::xoptional<xt::xtensor<double, 4>, bool> quad_prd;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> r_pearson;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> alpha;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> bias;
- // determine dimensions for inner components, intermediate
- // elements, for time average elements, and for final metrics:
- // -> inner components shape: (samples, subsets, series)
- inner_dims = {n_msk, n_exp, n_srs};
- // -> intermediate elements shape: (samples, subsets, series, time)
- inter_dims = {n_msk, n_exp, n_srs, n_tim};
- // -> time average elements shape: (samples, subsets, series, 1)
- mean_dims = {n_msk, n_exp, n_srs, 1};
- // -> final metrics shape: (series, subsets, samples)
- final_dims = {n_srs, n_msk, n_exp};
-
- // drop time steps where observations or predictions are NaN
- auto obs_nan = xt::isnan(q_obs);
- auto prd_nan = xt::isnan(q_prd);
+ // members for evaluation metrics
+ xtl::xoptional<xt::xtensor<double, 3>, bool> RMSE;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> NSE;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> KGE;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> KGEPRIME;
- t_msk = xt::ones<bool>({n_msk, n_srs, n_tim});
- xt::view(t_msk, xt::all()) =
- xt::view(msk, xt::all(), xt::newaxis(), xt::all());
- for (std::size_t m = 0; m < n_msk; m++)
+ // methods to compute elements
+ xt::xtensor<double, 4> get_mean_obs()
+ {
+ if (!mean_obs.has_value())
{
- xt::view(t_msk, m) =
- xt::where(obs_nan || prd_nan,
- false, xt::view(t_msk, m));
+ mean_obs = elements::calc_mean_obs<D2>(
+ q_obs, t_msk, b_exp, n_srs, n_msk, n_exp
+ );
}
+ return mean_obs.value();
};
- // members for evaluation metrics
- xt::xtensor<double, 3> RMSE;
- xt::xtensor<double, 3> NSE;
- xt::xtensor<double, 3> KGE;
- xt::xtensor<double, 3> KGEPRIME;
-
- // methods to compute elements
- void calc_mean_obs();
- void calc_mean_prd();
- void calc_quad_err();
- void calc_quad_obs();
- void calc_quad_prd();
- void calc_r_pearson();
- void calc_alpha();
- void calc_bias();
-
- // methods to compute metrics
- void calc_RMSE();
- void calc_NSE();
- void calc_KGE();
- void calc_KGEPRIME();
- };
-
- // Compute the mean of the observations.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \assign mean_obs:
- // Mean observed streamflow.
- // shape: (subsets, samples, series, 1)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_mean_obs()
- {
- mean_obs = xt::zeros<double>(mean_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 4> get_mean_prd()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!mean_prd.has_value())
{
- // apply the bootstrap sampling
- auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
- xt::view(mean_obs, m, e) =
- xt::nanmean(obs, -1, xt::keep_dims);
+ mean_prd = elements::calc_mean_prd<D2>(
+ q_prd, t_msk, b_exp, n_srs, n_msk, n_exp
+ );
}
- }
- }
+ return mean_prd.value();
+ };
- // Compute the mean of the predictions.
- //
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \assign mean_prd:
- // Mean predicted streamflow.
- // shape: (subsets, samples, series, 1)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_mean_prd()
- {
- mean_prd = xt::zeros<double>(mean_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 2> get_quad_err()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!quad_err.has_value())
{
- // apply the bootstrap sampling
- auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
- xt::view(mean_prd, m, e) =
- xt::nanmean(prd, -1, xt::keep_dims);
+ quad_err = elements::calc_quad_err<D2>(
+ q_obs, q_prd
+ );
}
- }
- }
-
- // Compute the quadratic error between observations and predictions.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \assign quad_err:
- // Quadratic errors between observations and predictions.
- // shape: (series, time)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_quad_err()
- {
- quad_err = xt::square(q_obs - q_prd);
- }
+ return quad_err.value();
+ };
- // Compute the quadratic error between observations and mean observation.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \require mean_obs:
- // Mean observed streamflow.
- // shape: (samples, series, 1)
- // \assign quad_obs:
- // Quadratic errors between observations and mean observation.
- // shape: (subsets, samples, series, time)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_quad_obs()
- {
- quad_obs = xt::zeros<double>(inter_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 4> get_quad_obs()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
-
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!quad_obs.has_value())
{
- xt::view(quad_obs, m, e) = xt::square(
- obs_masked - xt::view(mean_obs, m, e)
+ quad_obs = elements::calc_quad_obs<D2>(
+ q_obs, get_mean_obs(), t_msk,
+ n_srs, n_tim, n_msk, n_exp
);
}
- }
- }
+ return quad_obs.value();
+ };
- // Compute the quadratic error between predictions and mean prediction.
- //
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \require mean_prd:
- // Mean predicted streamflow.
- // shape: (samples, series, 1)
- // \assign quad_prd:
- // Quadratic errors between predictions and mean prediction.
- // shape: (subsets, samples, series, time)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_quad_prd()
- {
- quad_prd = xt::zeros<double>(inter_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 4> get_quad_prd()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
-
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!quad_prd.has_value())
{
- xt::view(quad_prd, m, e) = xt::square(
- prd_masked - xt::view(mean_prd, m, e)
+ quad_prd = elements::calc_quad_prd<D2>(
+ q_prd, get_mean_prd(), t_msk,
+ n_srs, n_tim, n_msk, n_exp
);
}
- }
- }
+ return quad_prd.value();
+ };
- // Compute the Pearson correlation coefficient.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \require mean_obs:
- // Mean observed streamflow.
- // shape: (samples, series, 1)
- // \require mean_prd:
- // Mean predicted streamflow.
- // shape: (samples, series, 1)
- // \require quad_obs:
- // Quadratic errors between observations and mean observation.
- // shape: (samples, series, time)
- // \require quad_prd:
- // Quadratic errors between predictions and mean prediction.
- // shape: (samples, series, time)
- // \assign r_pearson:
- // Pearson correlation coefficients.
- // shape: (subsets, samples, series)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_r_pearson()
- {
- // calculate error in timing and dynamics $r_{pearson}$
- // (Pearson's correlation coefficient)
- r_pearson = xt::zeros<double>(inner_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_r_pearson()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!r_pearson.has_value())
{
- auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
- auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
- auto r_num = xt::nansum(
- (prd - xt::view(mean_prd, m, e))
- * (obs - xt::view(mean_obs, m, e)),
- -1
+ r_pearson = elements::calc_r_pearson<D2>(
+ q_obs, q_prd, get_mean_obs(), get_mean_prd(),
+ get_quad_obs(), get_quad_prd(), t_msk, b_exp,
+ n_srs, n_msk, n_exp
);
+ }
+ return r_pearson.value();
+ };
- auto prd2 = xt::view(quad_prd, m, e, xt::all(), b_exp[e]);
- auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
- auto r_den = xt::sqrt(
- xt::nansum(prd2, -1) * xt::nansum(obs2, -1)
+ xt::xtensor<double, 3> get_alpha()
+ {
+ if (!alpha.has_value())
+ {
+ alpha = elements::calc_alpha<D2>(
+ q_obs, q_prd, get_mean_obs(), get_mean_prd(),
+ t_msk, b_exp, n_srs, n_msk, n_exp
);
-
- xt::view(r_pearson, m, e) = r_num / r_den;
}
- }
- }
+ return alpha.value();
+ };
- // Compute alpha.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \require mean_obs:
- // Mean observed streamflow.
- // shape: (samples, series, 1)
- // \require mean_prd:
- // Mean predicted streamflow.
- // shape: (samples, series, 1)
- // \assign alpha:
- // Alphas, ratios of standard deviations.
- // shape: (subsets, samples, series)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_alpha()
- {
- // calculate error in spread of flow $alpha$
- alpha = xt::zeros<double>(inner_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_bias()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!bias.has_value())
{
- auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
- auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
- xt::view(alpha, m, e) =
- evalhyd::maths::nanstd(prd, xt::view(mean_prd, m, e))
- / evalhyd::maths::nanstd(obs, xt::view(mean_obs, m, e));
+ bias = elements::calc_bias<D2>(
+ q_obs, q_prd, t_msk, b_exp, n_srs, n_msk, n_exp
+ );
}
- }
- }
+ return bias.value();
+ };
- // Compute the bias.
- //
- // \require q_obs:
- // Streamflow observations.
- // shape: (series, time)
- // \require q_prd:
- // Streamflow predictions.
- // shape: (series, time)
- // \assign bias:
- // Biases.
- // shape: (subsets, samples, series)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_bias()
- {
- // calculate $bias$
- bias = xt::zeros<double>(inner_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ public:
+ // constructor method
+ Evaluator(const D2& obs,
+ const D2& prd,
+ const B2& msk,
+ const std::vector<xt::xkeep_slice<int>>& exp) :
+ q_obs{obs}, q_prd{prd}, b_exp{exp}
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto prd_masked = xt::where(xt::view(t_msk, m), q_prd, NAN);
- auto obs_masked = xt::where(xt::view(t_msk, m), q_obs, NAN);
+ // determine size for recurring dimensions
+ n_tim = q_prd.shape(1);
+ n_srs = q_prd.shape(0);
+ n_msk = msk.shape(0);
+ n_exp = b_exp.size();
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ // drop time steps where observations or predictions are NaN
+ auto obs_nan = xt::isnan(q_obs);
+ auto prd_nan = xt::isnan(q_prd);
+
+ t_msk = xt::ones<bool>({n_msk, n_srs, n_tim});
+ xt::view(t_msk, xt::all()) =
+ xt::view(msk, xt::all(), xt::newaxis(), xt::all());
+ for (std::size_t m = 0; m < n_msk; m++)
{
- auto prd = xt::view(prd_masked, xt::all(), b_exp[e]);
- auto obs = xt::view(obs_masked, xt::all(), b_exp[e]);
- xt::view(bias, m, e) =
- xt::nansum(prd, -1) / xt::nansum(obs, -1);
+ xt::view(t_msk, m) =
+ xt::where(obs_nan || prd_nan,
+ false, xt::view(t_msk, m));
}
- }
- }
+ };
- // Compute the root-mean-square error (RMSE).
- //
- // \require quad_err:
- // Quadratic errors between observations and predictions.
- // shape: (series, time)
- // \assign RMSE:
- // Root-mean-square errors.
- // shape: (series, subsets, samples)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_RMSE()
- {
- // compute RMSE
- RMSE = xt::zeros<double>(final_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ // methods to compute metrics
+ xt::xtensor<double, 3> get_RMSE()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!RMSE.has_value())
{
- auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
- xt::view(RMSE, xt::all(), m, e) = xt::sqrt(xt::nanmean(err2, -1));
+ RMSE = metrics::calc_RMSE(
+ get_quad_err(), t_msk, b_exp, n_srs, n_msk, n_exp
+ );
}
- }
- }
+ return RMSE.value();
+ };
- // Compute the Nash-Sutcliffe Efficiency (NSE).
- //
- // \require quad_err:
- // Quadratic errors between observations and predictions.
- // shape: (series, time)
- // \require quad_obs:
- // Quadratic errors between observations and mean observation.
- // shape: (samples, series, time)
- // \assign NSE:
- // Nash-Sutcliffe efficiencies.
- // shape: (series, subsets, samples)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_NSE()
- {
- NSE = xt::zeros<double>(final_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_NSE()
{
- // apply the mask
- // (using NaN workaround until reducers work on masked_view)
- auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
-
- // compute variable one sample at a time
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!NSE.has_value())
{
- // compute squared errors operands
- auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
- xt::xtensor<double, 1> f_num =
- xt::nansum(err2, -1);
- auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
- xt::xtensor<double, 1> f_den =
- xt::nansum(obs2, -1);
-
- // compute NSE
- xt::view(NSE, xt::all(), m, e) = 1 - (f_num / f_den);
+ NSE = metrics::calc_NSE(
+ get_quad_err(), get_quad_obs(), t_msk, b_exp,
+ n_srs, n_msk, n_exp
+ );
}
- }
- }
+ return NSE.value();
+ };
- // Compute the Kling-Gupta Efficiency (KGE).
- //
- // \require r_pearson:
- // Pearson correlation coefficients.
- // shape: (samples, series)
- // \require alpha:
- // Alphas, ratios of standard deviations.
- // shape: (samples, series)
- // \require bias:
- // Biases.
- // shape: (samples, series)
- // \assign KGE:
- // Kling-Gupta efficiencies.
- // shape: (series, subsets, samples)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_KGE()
- {
- KGE = xt::zeros<double>(final_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_KGE()
{
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!KGE.has_value())
{
- // compute KGE
- xt::view(KGE, xt::all(), m, e) = 1 - xt::sqrt(
- xt::square(xt::view(r_pearson, m, e) - 1)
- + xt::square(xt::view(alpha, m, e) - 1)
- + xt::square(xt::view(bias, m, e) - 1)
+ KGE = metrics::calc_KGE(
+ get_r_pearson(), get_alpha(), get_bias(),
+ n_srs, n_msk, n_exp
);
}
- }
- }
+ return KGE.value();
+ };
- // Compute the modified Kling-Gupta Efficiency (KGEPRIME).
- //
- // \require mean_obs:
- // Mean observed streamflow.
- // shape: (samples, series, 1)
- // \require mean_prd:
- // Mean predicted streamflow.
- // shape: (samples, series, 1)
- // \require r_pearson:
- // Pearson correlation coefficients.
- // shape: (samples, series)
- // \require alpha:
- // Alphas, ratios of standard deviations.
- // shape: (samples, series)
- // \require bias:
- // Biases.
- // shape: (samples, series)
- // \assign KGEPRIME:
- // Modified Kling-Gupta efficiencies.
- // shape: (series, subsets, samples)
- template <class D2, class B2>
- void Evaluator<D2, B2>::calc_KGEPRIME()
- {
- KGEPRIME = xt::zeros<double>(final_dims);
- for (std::size_t m = 0; m < n_msk; m++)
+ xt::xtensor<double, 3> get_KGEPRIME()
{
- for (std::size_t e = 0; e < n_exp; e++)
+ if (!KGEPRIME.has_value())
{
- // calculate error in spread of flow $gamma$
- auto gamma = xt::view(alpha, m, e)
- * (xt::view(mean_obs, m, e, xt::all(), 0)
- / xt::view(mean_prd, m, e, xt::all(), 0));
-
- // compute KGEPRIME
- xt::view(KGEPRIME, xt::all(), m, e) = 1 - xt::sqrt(
- xt::square(xt::view(r_pearson, m, e) - 1)
- + xt::square(gamma - 1)
- + xt::square(xt::view(bias, m, e) - 1)
+ KGEPRIME = metrics::calc_KGEPRIME(
+ get_mean_obs(), get_mean_prd(),
+ get_r_pearson(), get_alpha(), get_bias(),
+ n_srs, n_msk, n_exp
);
}
- }
- }
+ return KGEPRIME.value();
+ };
+ };
}
}
diff --git a/include/evalhyd/detail/determinist/metrics.hpp b/include/evalhyd/detail/determinist/metrics.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..23f93429194442bc5e5b7a17a9e03b5a8c49da58
--- /dev/null
+++ b/include/evalhyd/detail/determinist/metrics.hpp
@@ -0,0 +1,238 @@
+#ifndef EVALHYD_DETERMINIST_METRICS_HPP
+#define EVALHYD_DETERMINIST_METRICS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xoperation.hpp>
+
+
+namespace evalhyd
+{
+ namespace determinist
+ {
+ namespace metrics
+ {
+ // Compute the root-mean-square error (RMSE).
+ //
+ // \param quad_err
+ // Quadratic errors between observations and predictions.
+ // shape: (series, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Root-mean-square errors.
+ // shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_RMSE(
+ const xt::xtensor<double, 2>& quad_err,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // compute RMSE
+ xt::xtensor<double, 3> RMSE =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
+ xt::view(RMSE, xt::all(), m, e) = xt::sqrt(xt::nanmean(err2, -1));
+ }
+ }
+
+ return RMSE;
+ }
+
+ // Compute the Nash-Sutcliffe Efficiency (NSE).
+ //
+ // \param quad_err
+ // Quadratic errors between observations and predictions.
+ // shape: (series, time)
+ // \param quad_obs
+ // Quadratic errors between observations and mean observation.
+ // shape: (subsets, samples, series, time)
+ // \param t_msk
+ // Temporal subsets of the whole record.
+ // shape: (subsets, series, time)
+ // \param b_exp
+ // Boostrap samples.
+ // shape: (samples, time slice)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Nash-Sutcliffe efficiencies.
+ // shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_NSE(
+ const xt::xtensor<double, 2>& quad_err,
+ const xt::xtensor<double, 4>& quad_obs,
+ const xt::xtensor<bool, 3>& t_msk,
+ const std::vector<xt::xkeep_slice<int>>& b_exp,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 3> NSE =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // apply the mask
+ // (using NaN workaround until reducers work on masked_view)
+ auto quad_err_masked = xt::where(xt::view(t_msk, m), quad_err, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // compute squared errors operands
+ auto err2 = xt::view(quad_err_masked, xt::all(), b_exp[e]);
+ xt::xtensor<double, 1> f_num =
+ xt::nansum(err2, -1);
+ auto obs2 = xt::view(quad_obs, m, e, xt::all(), b_exp[e]);
+ xt::xtensor<double, 1> f_den =
+ xt::nansum(obs2, -1);
+
+ // compute NSE
+ xt::view(NSE, xt::all(), m, e) = 1 - (f_num / f_den);
+ }
+ }
+
+ return NSE;
+ }
+
+ // Compute the Kling-Gupta Efficiency (KGE).
+ //
+ // \param r_pearson
+ // Pearson correlation coefficients.
+ // shape: (subsets, samples, series)
+ // \param alpha
+ // Alphas, ratios of standard deviations.
+ // shape: (subsets, samples, series)
+ // \param bias
+ // Biases.
+ // shape: (subsets, samples, series)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Kling-Gupta efficiencies.
+ // shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_KGE(
+ const xt::xtensor<double, 3>& r_pearson,
+ const xt::xtensor<double, 3>& alpha,
+ const xt::xtensor<double, 3>& bias,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 3> KGE =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // compute KGE
+ xt::view(KGE, xt::all(), m, e) = 1 - xt::sqrt(
+ xt::square(xt::view(r_pearson, m, e) - 1)
+ + xt::square(xt::view(alpha, m, e) - 1)
+ + xt::square(xt::view(bias, m, e) - 1)
+ );
+ }
+ }
+
+ return KGE;
+ }
+
+ // Compute the modified Kling-Gupta Efficiency (KGEPRIME).
+ //
+ // \param mean_obs
+ // Mean observed streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param mean_prd
+ // Mean predicted streamflow.
+ // shape: (subsets, samples, series, 1)
+ // \param r_pearson
+ // Pearson correlation coefficients.
+ // shape: (subsets, samples, series)
+ // \param alpha
+ // Alphas, ratios of standard deviations.
+ // shape: (subsets, samples, series)
+ // \param bias
+ // Biases.
+ // shape: (subsets, samples, series)
+ // \param n_srs
+ // Number of prediction series.
+ // \param n_msk
+ // Number of temporal subsets.
+ // \param n_exp
+ // Number of bootstrap samples.
+ // \return
+ // Modified Kling-Gupta efficiencies.
+ // shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_KGEPRIME(
+ const xt::xtensor<double, 4>& mean_obs,
+ const xt::xtensor<double, 4>& mean_prd,
+ const xt::xtensor<double, 3>& r_pearson,
+ const xt::xtensor<double, 3>& alpha,
+ const xt::xtensor<double, 3>& bias,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ xt::xtensor<double, 3> KGEPRIME =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ // calculate error in spread of flow $gamma$
+ auto gamma = xt::view(alpha, m, e)
+ * (xt::view(mean_obs, m, e, xt::all(), 0)
+ / xt::view(mean_prd, m, e, xt::all(), 0));
+
+ // compute KGEPRIME
+ xt::view(KGEPRIME, xt::all(), m, e) = 1 - xt::sqrt(
+ xt::square(xt::view(r_pearson, m, e) - 1)
+ + xt::square(gamma - 1)
+ + xt::square(xt::view(bias, m, e) - 1)
+ );
+ }
+ }
+
+ return KGEPRIME;
+ }
+ }
+ }
+}
+
+#endif //EVALHYD_DETERMINIST_METRICS_HPP
diff --git a/include/evalhyd/detail/utils.hpp b/include/evalhyd/detail/utils.hpp
index 1d11b4bd1274a4b281b80d5a472b56f85454d6e2..64934cab69a06e81fe19bbb923578244f2d06b51 100644
--- a/include/evalhyd/detail/utils.hpp
+++ b/include/evalhyd/detail/utils.hpp
@@ -12,67 +12,6 @@ namespace evalhyd
{
namespace utils
{
- /// Procedure to determine based on a list of metrics which elements
- /// and which metrics (and their associated elements) require to be
- /// pre-computed for memoisation purposes.
- ///
- /// \param [in] metrics:
- /// Vector of strings for the metric(s) to be computed.
- /// \param [in] elements:
- /// Map between metrics and their required computation elements.
- /// \param [in] dependencies:
- /// Map between metrics and their dependencies.
- /// \param [out] required_elements:
- /// Set of elements identified as required to be pre-computed.
- /// \param [out] required_dependencies:
- /// Set of metrics identified as required to be pre-computed.
- inline void find_requirements (
- const std::vector<std::string>& metrics,
- std::unordered_map<std::string, std::vector<std::string>>& elements,
- std::unordered_map<std::string, std::vector<std::string>>& dependencies,
- std::vector<std::string>& required_elements,
- std::vector<std::string>& required_dependencies
- )
- {
- std::unordered_set<std::string> found_elements;
- std::unordered_set<std::string> found_dependencies;
-
- for (const auto& metric : metrics)
- {
- // add elements to pre-computation set
- for (const auto& element : elements[metric])
- {
- if (found_elements.find(element) == found_elements.end())
- {
- found_elements.insert(element);
- required_elements.push_back(element);
- }
- }
-
- // add metric dependencies to pre-computation set
- if (dependencies.find(metric) != dependencies.end())
- {
- for (const auto& dependency : dependencies[metric])
- {
- if (found_dependencies.find(dependency) == found_dependencies.end())
- {
- found_dependencies.insert(dependency);
- required_dependencies.push_back(dependency);
- }
- // add dependency elements to pre-computation set
- for (const auto& element : elements[dependency])
- {
- if (found_elements.find(element) == found_elements.end())
- {
- found_elements.insert(element);
- required_elements.push_back(element);
- }
- }
- }
- }
- }
- }
-
/// Procedure to check that all elements in the list of metrics are
/// valid metrics.
///
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index f24254e522055a9b20615261cd5f3edebfcc4425..c7e7c149e2c2cd8568b514b89d330439a41859fc 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -340,67 +340,6 @@ namespace evalhyd
// instantiate determinist evaluator
determinist::Evaluator<D2, B2> evaluator(obs, prd, msk, exp);
- // declare maps for memoisation purposes
- std::unordered_map<std::string, std::vector<std::string>> elt;
- std::unordered_map<std::string, std::vector<std::string>> dep;
-
- // register potentially recurring computation elt across metrics
- elt["RMSE"] = {"quad_err"};
- elt["NSE"] = {"mean_obs", "quad_obs", "quad_err"};
- elt["KGE"] = {"mean_obs", "mean_prd", "quad_obs", "quad_prd",
- "r_pearson", "alpha", "bias"};
- elt["KGEPRIME"] = {"mean_obs", "mean_prd", "quad_obs", "quad_prd",
- "r_pearson", "alpha", "bias"};
-
- // TODO: register nested metrics (i.e. metric dependent on another metric)
- // dep[...] = {...};
-
- // determine required elt/dep to be pre-computed
- std::vector<std::string> req_elt;
- std::vector<std::string> req_dep;
-
- utils::find_requirements(metrics, elt, dep, req_elt, req_dep);
-
- // pre-compute required elt
- for ( const auto& element : req_elt )
- {
- if ( element == "mean_obs" )
- {
- evaluator.calc_mean_obs();
- }
- else if ( element == "mean_prd" )
- {
- evaluator.calc_mean_prd();
- }
- else if ( element == "quad_err" )
- {
- evaluator.calc_quad_err();
- }
- else if ( element == "quad_obs" )
- {
- evaluator.calc_quad_obs();
- }
- else if ( element == "quad_prd" )
- {
- evaluator.calc_quad_prd();
- }
- else if ( element == "r_pearson" )
- {
- evaluator.calc_r_pearson();
- }
- else if ( element == "alpha" )
- {
- evaluator.calc_alpha();
- }
- else if ( element == "bias" )
- {
- evaluator.calc_bias();
- }
- }
-
- // TODO: pre-compute required dep
- // for ( const auto& dependency : req_dep ) {...}
-
// retrieve or compute requested metrics
std::vector<xt::xarray<double>> r;
@@ -410,39 +349,27 @@ namespace evalhyd
{
if ( metric == "RMSE" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_RMSE();
- }
- r.emplace_back(uncertainty::summarise(evaluator.RMSE, summary));
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_RMSE(), summary)
+ );
}
else if ( metric == "NSE" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_NSE();
- }
- r.emplace_back(uncertainty::summarise(evaluator.NSE, summary));
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_NSE(), summary)
+ );
}
else if ( metric == "KGE" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_KGE();
- }
- r.emplace_back(uncertainty::summarise(evaluator.KGE, summary));
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_KGE(), summary)
+ );
}
else if ( metric == "KGEPRIME" )
{
- if (std::find(req_dep.begin(), req_dep.end(), metric)
- == req_dep.end())
- {
- evaluator.calc_KGEPRIME();
- }
- r.emplace_back(uncertainty::summarise(evaluator.KGEPRIME, summary));
+ r.emplace_back(
+ uncertainty::summarise(evaluator.get_KGEPRIME(), summary)
+ );
}
}