diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp index 0bebd065b0f38b469463c7675b605ec174fed19c..cb21e5ff2062ed3e84a4d6db555f7d3b96788636 100644 --- a/include/evalhyd/evald.hpp +++ b/include/evalhyd/evald.hpp @@ -46,22 +46,31 @@ namespace evalhyd /// ========================= ================================= /// transformations details /// ========================= ================================= - /// ``"inv"`` The reciprocal function - /// **f(Q) = 1/Q** is applied. /// ``"sqrt"`` The square root function /// **f(Q) = √Q** is applied. + /// ``"pow"`` The power function + /// **f(Q) = Q^n** is applied (where + /// **n** can be set through the + /// *exponent* parameter). + /// ``"inv"`` The reciprocal function + /// **f(Q) = 1/Q** is applied. /// ``"log"`` The natural logarithm function /// **f(Q) = ln(Q)** is applied. /// ========================= ================================= /// + /// exponent: ``xt::xscalar<double>``, optional + /// The value of the exponent n to use when the *transform* is the + /// power function. If not provided (or set to default value 1), + /// the streamflow observations and predictions remain untransformed. + /// /// epsilon: ``xt::xscalar<double>``, optional /// The value of the small constant ε to add to both the streamflow /// observations and predictions prior to the calculation of the /// *metrics* when the *transform* is the reciprocal function or - /// the natural logarithm since neither are defined for 0. If not - /// provided, set to default value equal to one hundredth of the - /// mean of the streamflow observations, as recommended by - /// `Pushpalatha et al. (2012) + /// the natural logarithm (since neither are defined for 0). If not + /// provided (or set to default value -9), one hundredth of the mean + /// of the streamflow observations is used as value for epsilon, as + /// recommended by `Pushpalatha et al. (2012) /// <https://doi.org/10.1016/j.jhydrol.2011.11.055>`_. /// /// :Returns: @@ -90,7 +99,11 @@ namespace evalhyd /// /// .. code-block:: c++ /// - /// evalhyd::evald(obs, prd, {"NSE"}, "log", 0.05); + /// evalhyd::evald(obs, prd, {"NSE"}, "pow", 0.2); + /// + /// .. code-block:: c++ + /// + /// evalhyd::evald(obs, prd, {"NSE"}, "log", 1, 0.05); /// /// \endrst template <class A> @@ -99,6 +112,7 @@ namespace evalhyd const xt::xexpression<A>& q_prd, const std::vector<std::string>& metrics, const std::string& transform="none", + const xt::xscalar<double> exponent=1, xt::xscalar<double> epsilon=-9 ) { @@ -106,10 +120,11 @@ namespace evalhyd A prd; // apply streamflow transformation if requested - if ( transform == "none" ) - { - obs = std::move(q_obs.derived_cast()); - prd = std::move(q_prd.derived_cast()); + if ( transform == "none" ) { + obs = std::move( + q_obs.derived_cast()); + prd = std::move( + q_prd.derived_cast()); } else if ( transform == "sqrt" ) { @@ -138,6 +153,14 @@ namespace evalhyd obs = xt::log(q_obs.derived_cast() + epsilon); prd = xt::log(q_prd.derived_cast() + epsilon); } + else if ( transform == "pow" ) + { + if ( exponent != 1 ) + { + obs = xt::pow(q_obs.derived_cast(), exponent); + prd = xt::pow(q_prd.derived_cast(), exponent); + } + } else { throw std::runtime_error( diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp index 088f3b09aead2e1ee6f96d92ef5e75ac29451f8d..dc7689ed2a40840c132123449b8aa823a837bad6 100644 --- a/tests/test_determinist.cpp +++ b/tests/test_determinist.cpp @@ -9,58 +9,56 @@ #include "evalhyd/evald.hpp" TEST(DeterministTests, TestMetrics) { -// read in data -std::ifstream ifs; -ifs.open("./data/q_obs.csv"); -xt::xtensor<double, 2> observed =xt::load_csv<int>(ifs); -ifs.close(); - -ifs.open("./data/q_prd.csv"); -xt::xtensor<double, 2> predicted_2d = xt::view( - xt::load_csv<double>(ifs), xt::range(0, 5), xt::all() -); -ifs.close(); - -xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0); - -// compute scores (both with 2D and 1D tensors) -std::vector<xt::xarray<double>> metrics = - evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted_2d, {"RMSE", "NSE", "KGE", "KGEPRIME"} - ); - -// check results on all metrics -xt::xtensor<double, 2> rmse = - {{777.034272}, - {776.878479}, - {777.800217}, - {778.151082}, - {778.61487}}; -EXPECT_TRUE(xt::allclose(metrics[0], rmse)); - -xt::xtensor<double, 2> nse = - {{0.718912}, - {0.719025}, - {0.718358}, - {0.718104}, - {0.717767}}; -EXPECT_TRUE(xt::allclose(metrics[1], nse)); - -xt::xtensor<double, 2> kge = - {{0.748088}, - {0.746106}, - {0.744111}, - {0.743011}, - {0.741768}}; -EXPECT_TRUE(xt::allclose(metrics[2], kge)); - -xt::xtensor<double, 2> kgeprime = - {{0.813141}, - {0.812775}, - {0.812032}, - {0.811787}, - {0.811387}}; -EXPECT_TRUE(xt::allclose(metrics[3], kgeprime)); + // read in data + std::ifstream ifs; + ifs.open("./data/q_obs.csv"); + xt::xtensor<double, 2> observed =xt::load_csv<int>(ifs); + ifs.close(); + + ifs.open("./data/q_prd.csv"); + xt::xtensor<double, 2> predicted = xt::view( + xt::load_csv<double>(ifs), xt::range(0, 5), xt::all() + ); + ifs.close(); + + // compute scores (both with 2D and 1D tensors) + std::vector<xt::xarray<double>> metrics = + evalhyd::evald<xt::xtensor<double, 2>>( + observed, predicted, {"RMSE", "NSE", "KGE", "KGEPRIME"} + ); + + // check results on all metrics + xt::xtensor<double, 2> rmse = + {{777.034272}, + {776.878479}, + {777.800217}, + {778.151082}, + {778.61487}}; + EXPECT_TRUE(xt::allclose(metrics[0], rmse)); + + xt::xtensor<double, 2> nse = + {{0.718912}, + {0.719025}, + {0.718358}, + {0.718104}, + {0.717767}}; + EXPECT_TRUE(xt::allclose(metrics[1], nse)); + + xt::xtensor<double, 2> kge = + {{0.748088}, + {0.746106}, + {0.744111}, + {0.743011}, + {0.741768}}; + EXPECT_TRUE(xt::allclose(metrics[2], kge)); + + xt::xtensor<double, 2> kgeprime = + {{0.813141}, + {0.812775}, + {0.812032}, + {0.811787}, + {0.811387}}; + EXPECT_TRUE(xt::allclose(metrics[3], kgeprime)); } TEST(DeterministTests, TestNSE_1d2d) { @@ -112,17 +110,15 @@ TEST(DeterministTests, TestNSE_transform) { ifs.close(); ifs.open("./data/q_prd.csv"); - xt::xtensor<double, 2> predicted_2d = xt::view( + xt::xtensor<double, 2> predicted = xt::view( xt::load_csv<double>(ifs), xt::range(0, 5), xt::all() ); ifs.close(); - xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0); - // compute and check results on square-rooted streamflow series std::vector<xt::xarray<double>> metrics = evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted_2d, {"NSE"}, "sqrt" + observed, predicted, {"NSE"}, "sqrt" ); xt::xtensor<double, 2> nse_sqrt = @@ -136,7 +132,7 @@ TEST(DeterministTests, TestNSE_transform) { // compute and check results on inverted streamflow series metrics = evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted_2d, {"NSE"}, "inv" + observed, predicted, {"NSE"}, "inv" ); xt::xtensor<double, 2> nse_inv = @@ -150,7 +146,7 @@ TEST(DeterministTests, TestNSE_transform) { // compute and check results on square-rooted streamflow series metrics = evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted_2d, {"NSE"}, "log" + observed, predicted, {"NSE"}, "log" ); xt::xtensor<double, 2> nse_log = @@ -161,4 +157,18 @@ TEST(DeterministTests, TestNSE_transform) { {0.893793}}; EXPECT_TRUE(xt::allclose(metrics[0], nse_log)); + // compute and check results on power-transformed streamflow series + metrics = + evalhyd::evald<xt::xtensor<double, 2>>( + observed, predicted, {"NSE"}, "pow", 0.2 + ); + + xt::xtensor<double, 2> nse_pow = + {{0.899207}, + {0.899395}, + {0.899451}, + {0.899578}, + {0.899588}}; + EXPECT_TRUE(xt::allclose(metrics[0], nse_pow)); + }