diff --git a/include/evalhyd/detail/determinist/errors.hpp b/include/evalhyd/detail/determinist/errors.hpp
index 4fb1ff894f36d73ba69291e5199e1cc7ae1bfc4a..f6e12354f9d45d05ebf34a2bef08793750fd50e1 100644
--- a/include/evalhyd/detail/determinist/errors.hpp
+++ b/include/evalhyd/detail/determinist/errors.hpp
@@ -121,7 +121,7 @@ namespace evalhyd
return mean_prd;
}
- /// Compute the quadratic error between observations and predictions.
+ /// Compute the error between observations and predictions.
///
/// \param q_obs
/// Streamflow observations.
@@ -130,15 +130,45 @@ namespace evalhyd
/// Streamflow predictions.
/// shape: (series, time)
/// \return
- /// Quadratic errors between observations and predictions.
+ /// Errors between observations and predictions.
/// shape: (series, time)
template <class XD2>
- inline xt::xtensor<double, 2> calc_quad_err(
+ inline xt::xtensor<double, 2> calc_err(
const XD2& q_obs,
const XD2& q_prd
)
{
- return xt::square(q_obs - q_prd);
+ return q_obs - q_prd;
+ }
+
+ /// Compute the absolute error between observations and predictions.
+ ///
+ /// \param err
+ /// Errors between observations and predictions.
+ /// shape: (series, time)
+ /// \return
+ /// Quadratic errors between observations and predictions.
+ /// shape: (series, time)
+ inline xt::xtensor<double, 2> calc_abs_err(
+ const xt::xtensor<double, 2>& err
+ )
+ {
+ return xt::abs(err);
+ }
+
+ /// Compute the quadratic error between observations and predictions.
+ ///
+ /// \param err
+ /// Errors between observations and predictions.
+ /// shape: (series, time)
+ /// \return
+ /// Quadratic errors between observations and predictions.
+ /// shape: (series, time)
+ inline xt::xtensor<double, 2> calc_quad_err(
+ const xt::xtensor<double, 2>& err
+ )
+ {
+ return xt::square(err);
}
/// Compute the quadratic error between observations and mean observation.
@@ -252,7 +282,101 @@ namespace evalhyd
namespace metrics
{
- /// Compute the root-mean-square error (RMSE).
+ /// Compute the mean absolute error (MAE).
+ ///
+ /// \param abs_err
+ /// Absolute errors between observations and predictions.
+ /// shape: (series, time)
+ /// \param t_msk
+ /// Temporal subsets of the whole record.
+ /// shape: (series, subsets, 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 absolute errors.
+ /// shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_MAE(
+ const xt::xtensor<double, 2>& abs_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> MAE =
+ 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 abs_err_masked = xt::where(xt::view(t_msk, xt::all(), m),
+ abs_err, NAN);
+
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ auto err = xt::view(abs_err_masked, xt::all(), b_exp[e]);
+ xt::view(MAE, xt::all(), m, e) = xt::nanmean(err, -1);
+ }
+ }
+
+ return MAE;
+ }
+
+ /// Compute the mean absolute relative error (MARE).
+ ///
+ /// \param MAE
+ /// Mean absolute errors.
+ /// shape: (series, subsets, samples)
+ /// \param mean_obs
+ /// Mean observed streamflow.
+ /// shape: (subsets, samples, series, 1)
+ /// \param n_srs
+ /// Number of prediction series.
+ /// \param n_msk
+ /// Number of temporal subsets.
+ /// \param n_exp
+ /// Number of bootstrap samples.
+ /// \return
+ /// Mean absolute relative errors.
+ /// shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_MARE(
+ const xt::xtensor<double, 3>& MAE,
+ const xt::xtensor<double, 4>& mean_obs,
+ std::size_t n_srs,
+ std::size_t n_msk,
+ std::size_t n_exp
+ )
+ {
+ // compute RMSE
+ xt::xtensor<double, 3> MARE =
+ xt::zeros<double>({n_srs, n_msk, n_exp});
+
+ for (std::size_t m = 0; m < n_msk; m++)
+ {
+ // compute variable one sample at a time
+ for (std::size_t e = 0; e < n_exp; e++)
+ {
+ xt::view(MARE, xt::all(), m, e) =
+ xt::view(MAE, xt::all(), m, e)
+ / xt::view(mean_obs, m, e, xt::all(), 0);
+ }
+ }
+
+ return MARE;
+ }
+
+ /// Compute the mean square error (MSE).
///
/// \param quad_err
/// Quadratic errors between observations and predictions.
@@ -270,9 +394,9 @@ namespace evalhyd
/// \param n_exp
/// Number of bootstrap samples.
/// \return
- /// Root-mean-square errors.
+ /// Mean square errors.
/// shape: (series, subsets, samples)
- inline xt::xtensor<double, 3> calc_RMSE(
+ inline xt::xtensor<double, 3> calc_MSE(
const xt::xtensor<double, 2>& quad_err,
const xt::xtensor<bool, 3>& t_msk,
const std::vector<xt::xkeep_slice<int>>& b_exp,
@@ -282,7 +406,7 @@ namespace evalhyd
)
{
// compute RMSE
- xt::xtensor<double, 3> RMSE =
+ xt::xtensor<double, 3> MSE =
xt::zeros<double>({n_srs, n_msk, n_exp});
for (std::size_t m = 0; m < n_msk; m++)
@@ -296,10 +420,28 @@ namespace evalhyd
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));
+ xt::view(MSE, xt::all(), m, e) = xt::nanmean(err2, -1);
}
}
+ return MSE;
+ }
+
+ /// Compute the root mean square error (RMSE).
+ ///
+ /// \param MSE
+ /// Mean square errors.
+ /// shape: (series, subsets, samples)
+ /// \return
+ /// Root mean square errors.
+ /// shape: (series, subsets, samples)
+ inline xt::xtensor<double, 3> calc_RMSE(
+ const xt::xtensor<double, 3>& MSE
+ )
+ {
+ // compute RMSE
+ auto RMSE = xt::sqrt(MSE);
+
return RMSE;
}
}
diff --git a/include/evalhyd/detail/determinist/evaluator.hpp b/include/evalhyd/detail/determinist/evaluator.hpp
index f7e58bd52bdc42ed59d9673d8601942b1c234e94..a5d5a3e0abc5d13334fc430cb2f373f4937c9228 100644
--- a/include/evalhyd/detail/determinist/evaluator.hpp
+++ b/include/evalhyd/detail/determinist/evaluator.hpp
@@ -39,6 +39,8 @@ namespace evalhyd
// members for computational elements
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> err;
+ xtl::xoptional<xt::xtensor<double, 2>, bool> abs_err;
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;
@@ -48,6 +50,9 @@ namespace evalhyd
xtl::xoptional<xt::xtensor<double, 3>, bool> bias;
// members for evaluation metrics
+ xtl::xoptional<xt::xtensor<double, 3>, bool> MAE;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> MARE;
+ xtl::xoptional<xt::xtensor<double, 3>, bool> MSE;
xtl::xoptional<xt::xtensor<double, 3>, bool> RMSE;
xtl::xoptional<xt::xtensor<double, 3>, bool> NSE;
xtl::xoptional<xt::xtensor<double, 3>, bool> KGE;
@@ -78,12 +83,34 @@ namespace evalhyd
return mean_prd.value();
};
+ xt::xtensor<double, 2> get_err()
+ {
+ if (!err.has_value())
+ {
+ err = elements::calc_err(
+ q_obs, q_prd
+ );
+ }
+ return err.value();
+ };
+
+ xt::xtensor<double, 2> get_abs_err()
+ {
+ if (!abs_err.has_value())
+ {
+ abs_err = elements::calc_abs_err(
+ get_err()
+ );
+ }
+ return abs_err.value();
+ };
+
xt::xtensor<double, 2> get_quad_err()
{
if (!quad_err.has_value())
{
quad_err = elements::calc_quad_err(
- q_obs, q_prd
+ get_err()
);
}
return quad_err.value();
@@ -196,12 +223,45 @@ namespace evalhyd
};
// methods to compute metrics
+ xt::xtensor<double, 3> get_MAE()
+ {
+ if (!MAE.has_value())
+ {
+ MAE = metrics::calc_MAE(
+ get_abs_err(), t_msk, b_exp, n_srs, n_msk, n_exp
+ );
+ }
+ return MAE.value();
+ };
+
+ xt::xtensor<double, 3> get_MARE()
+ {
+ if (!MARE.has_value())
+ {
+ MARE = metrics::calc_MARE(
+ get_MAE(), get_mean_obs(), n_srs, n_msk, n_exp
+ );
+ }
+ return MARE.value();
+ };
+
+ xt::xtensor<double, 3> get_MSE()
+ {
+ if (!MSE.has_value())
+ {
+ MSE = metrics::calc_MSE(
+ get_quad_err(), t_msk, b_exp, n_srs, n_msk, n_exp
+ );
+ }
+ return MSE.value();
+ };
+
xt::xtensor<double, 3> get_RMSE()
{
if (!RMSE.has_value())
{
RMSE = metrics::calc_RMSE(
- get_quad_err(), t_msk, b_exp, n_srs, n_msk, n_exp
+ get_MSE()
);
}
return RMSE.value();
diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 810df0ad4bfc94834d7a399848da829ac53eeea9..ef858dfcfe97d058e7351155dd711f152a9f3d1b 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -208,7 +208,8 @@ namespace evalhyd
// check that the metrics to be computed are valid
utils::check_metrics(
metrics,
- {"RMSE", "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D"}
+ {"MAE", "MARE", "MSE", "RMSE",
+ "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D"}
);
// check that optional parameters are valid
@@ -434,6 +435,24 @@ namespace evalhyd
for ( const auto& metric : metrics )
{
+ if ( metric == "MAE" )
+ {
+ r.emplace_back(
+ uncertainty::summarise_d(evaluator.get_MAE(), summary)
+ );
+ }
+ if ( metric == "MARE" )
+ {
+ r.emplace_back(
+ uncertainty::summarise_d(evaluator.get_MARE(), summary)
+ );
+ }
+ if ( metric == "MSE" )
+ {
+ r.emplace_back(
+ uncertainty::summarise_d(evaluator.get_MSE(), summary)
+ );
+ }
if ( metric == "RMSE" )
{
r.emplace_back(
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 78f07d8169b05b743a099168eab49ffbad619eb1..41df12857ebd3a5f36ac569c7562b5c4a03c1a2b 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -15,6 +15,7 @@
#include <xtensor/xmath.hpp>
#include <xtensor/xsort.hpp>
#include <xtensor/xcsv.hpp>
+#include <xtensor/xio.hpp>
#include "evalhyd/evald.hpp"
@@ -26,7 +27,8 @@ using namespace xt::placeholders; // required for `_` to work
std::vector<std::string> all_metrics_d = {
- "RMSE", "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D"
+ "MAE", "MARE", "MSE", "RMSE",
+ "NSE", "KGE", "KGE_D", "KGEPRIME", "KGEPRIME_D"
};
std::tuple<xt::xtensor<double, 2>, xt::xtensor<double, 2>> load_data_d()
@@ -60,13 +62,37 @@ TEST(DeterministTests, TestMetrics)
);
// check results on all metrics
+ xt::xtensor<double, 3> mae =
+ {{{ 265.192926}},
+ {{ 265.681672}},
+ {{ 265.70418 }},
+ {{ 265.70418 }},
+ {{ 265.836013}}};
+ EXPECT_TRUE(xt::allclose(metrics[0], mae));
+
+ xt::xtensor<double, 3> mare =
+ {{{ 0.211088}},
+ {{ 0.211477}},
+ {{ 0.211495}},
+ {{ 0.211495}},
+ {{ 0.2116 }}};
+ EXPECT_TRUE(xt::allclose(metrics[1], mare));
+
+ xt::xtensor<double, 3> mse =
+ {{{ 603782.26045 }},
+ {{ 603540.170418}},
+ {{ 604973.176849}},
+ {{ 605519.106109}},
+ {{ 606241.115756}}};
+ EXPECT_TRUE(xt::allclose(metrics[2], mse));
+
xt::xtensor<double, 3> rmse =
{{{777.034272}},
{{776.878479}},
{{777.800217}},
{{778.151082}},
{{778.61487 }}};
- EXPECT_TRUE(xt::allclose(metrics[0], rmse));
+ EXPECT_TRUE(xt::allclose(metrics[3], rmse));
xt::xtensor<double, 3> nse =
{{{0.718912}},
@@ -74,7 +100,7 @@ TEST(DeterministTests, TestMetrics)
{{0.718358}},
{{0.718104}},
{{0.717767}}};
- EXPECT_TRUE(xt::allclose(metrics[1], nse));
+ EXPECT_TRUE(xt::allclose(metrics[4], nse));
xt::xtensor<double, 3> kge =
{{{0.748088}},
@@ -82,7 +108,7 @@ TEST(DeterministTests, TestMetrics)
{{0.744111}},
{{0.743011}},
{{0.741768}}};
- EXPECT_TRUE(xt::allclose(metrics[2], kge));
+ EXPECT_TRUE(xt::allclose(metrics[5], kge));
xt::xtensor<double, 4> kge_d =
{{{{ 0.907125, 1.224429, 1.066824}}},
@@ -90,7 +116,7 @@ TEST(DeterministTests, TestMetrics)
{{{ 0.90805 , 1.228409, 1.069668}}},
{{{ 0.908256, 1.229451, 1.070558}}},
{{{ 0.908474, 1.230672, 1.071395}}}};
- EXPECT_TRUE(xt::allclose(metrics[3], kge_d));
+ EXPECT_TRUE(xt::allclose(metrics[6], kge_d));
xt::xtensor<double, 3> kgeprime =
{{{0.813141}},
@@ -98,7 +124,7 @@ TEST(DeterministTests, TestMetrics)
{{0.812032}},
{{0.811787}},
{{0.811387}}};
- EXPECT_TRUE(xt::allclose(metrics[4], kgeprime));
+ EXPECT_TRUE(xt::allclose(metrics[7], kgeprime));
xt::xtensor<double, 4> kgeprime_d =
{{{{ 0.907125, 1.147733, 1.066824}}},
@@ -106,7 +132,7 @@ TEST(DeterministTests, TestMetrics)
{{{ 0.90805 , 1.148403, 1.069668}}},
{{{ 0.908256, 1.148421, 1.070558}}},
{{{ 0.908474, 1.148663, 1.071395}}}};
- EXPECT_TRUE(xt::allclose(metrics[5], kgeprime_d));
+ EXPECT_TRUE(xt::allclose(metrics[8], kgeprime_d));
}
TEST(DeterministTests, TestTransform)