diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp index 5e629998c12de5bfbe6bf6490214cfb6c770d446..0bebd065b0f38b469463c7675b605ec174fed19c 100644 --- a/include/evalhyd/evald.hpp +++ b/include/evalhyd/evald.hpp @@ -3,8 +3,10 @@ #include <unordered_map> #include <vector> +#include <stdexcept> #include <xtensor/xexpression.hpp> #include <xtensor/xarray.hpp> +#include <xtensor/xscalar.hpp> #include "../../src/utils.hpp" #include "../../src/determinist/evaluator.hpp" @@ -19,26 +21,52 @@ namespace evalhyd /// /// :Template Parameters: /// - /// A: Any inherited class of :cpp:`xt::xcontainer<...>` + /// A: Any inherited class of ``xt::xcontainer<...>`` /// The type of data structure for the input data (e.g. - /// :cpp:`xt::xtensor<...>`, :cpp:`xt::xarray<...>`). + /// ``xt::xtensor<...>``, ``xt::xarray<...>``). /// /// :Parameters: /// - /// q_obs: :cpp:`xt::xexpression<A>` + /// q_obs: ``xt::xexpression<A>`` /// Streamflow observations. /// shape: (1+, time) /// - /// q_prd: :cpp:`xt::xexpression<A>` + /// q_prd: ``xt::xexpression<A>`` /// Streamflow predictions. /// shape: (1+, time) /// - /// metrics: :cpp:`std::vector<std::string>` + /// metrics: ``std::vector<std::string>`` /// The sequence of evaluation metrics to be computed. /// + /// transform: ``std::string``, optional + /// The transformation to apply to both streamflow observations and + /// predictions prior to the calculation of the *metrics*. The options + /// are listed in the table below. + /// + /// ========================= ================================= + /// transformations details + /// ========================= ================================= + /// ``"inv"`` The reciprocal function + /// **f(Q) = 1/Q** is applied. + /// ``"sqrt"`` The square root function + /// **f(Q) = √Q** is applied. + /// ``"log"`` The natural logarithm function + /// **f(Q) = ln(Q)** is applied. + /// ========================= ================================= + /// + /// 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) + /// <https://doi.org/10.1016/j.jhydrol.2011.11.055>`_. + /// /// :Returns: /// - /// :cpp:`std::vector<xt::xarray<double>>` + /// ``std::vector<xt::xarray<double>>`` /// The sequence of evaluation metrics computed /// in the same order as given in *metrics*. /// @@ -56,16 +84,66 @@ namespace evalhyd /// /// evalhyd::evald(obs, prd, {"NSE"}); /// + /// .. code-block:: c++ + /// + /// evalhyd::evald(obs, prd, {"NSE"}, "sqrt"); + /// + /// .. code-block:: c++ + /// + /// evalhyd::evald(obs, prd, {"NSE"}, "log", 0.05); + /// /// \endrst template <class A> std::vector<xt::xarray<double>> evald( const xt::xexpression<A>& q_obs, const xt::xexpression<A>& q_prd, - const std::vector<std::string>& metrics + const std::vector<std::string>& metrics, + const std::string& transform="none", + xt::xscalar<double> epsilon=-9 ) { - const A& obs = q_obs.derived_cast(); - const A& prd = q_prd.derived_cast(); + A obs; + A prd; + + // apply streamflow transformation if requested + if ( transform == "none" ) + { + obs = std::move(q_obs.derived_cast()); + prd = std::move(q_prd.derived_cast()); + } + else if ( transform == "sqrt" ) + { + obs = xt::sqrt(q_obs.derived_cast()); + prd = xt::sqrt(q_prd.derived_cast()); + } + else if ( transform == "inv" ) + { + if (epsilon == -9) + // determine an epsilon value to avoid zero divide + // following recommendation in Pushpalatha et al. (2012) + // doi.org/10.1016/j.jhydrol.2011.11.055 + epsilon = xt::mean(q_obs.derived_cast())() * 0.01; + + obs = 1. / (q_obs.derived_cast() + epsilon); + prd = 1. / (q_prd.derived_cast() + epsilon); + } + else if ( transform == "log" ) + { + if (epsilon == -9) + // determine an epsilon value to avoid zero divide + // following recommendation in Pushpalatha et al. (2012) + // doi.org/10.1016/j.jhydrol.2011.11.055 + epsilon = xt::mean(q_obs.derived_cast())() * 0.01; + + obs = xt::log(q_obs.derived_cast() + epsilon); + prd = xt::log(q_prd.derived_cast() + epsilon); + } + else + { + throw std::runtime_error( + "invalid streamflow transformation: " + transform + ); + } // check that the metrics to be computed are valid utils::check_metrics( diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp index a0078c59a616a1d5472d82fa1ee3b8ccbbeccbb7..088f3b09aead2e1ee6f96d92ef5e75ac29451f8d 100644 --- a/tests/test_determinist.cpp +++ b/tests/test_determinist.cpp @@ -8,7 +8,62 @@ #include "evalhyd/evald.hpp" -TEST(DeterministTests, TestNSE1D2D) { +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)); +} + +TEST(DeterministTests, TestNSE_1d2d) { // read in data std::ifstream ifs; ifs.open("./data/q_obs.csv"); @@ -49,7 +104,7 @@ TEST(DeterministTests, TestNSE1D2D) { EXPECT_TRUE(xt::allclose(metrics_1d[0], nse[0])); } -TEST(DeterministTests, TestMetrics) { +TEST(DeterministTests, TestNSE_transform) { // read in data std::ifstream ifs; ifs.open("./data/q_obs.csv"); @@ -64,42 +119,46 @@ TEST(DeterministTests, TestMetrics) { xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0); - // compute scores (both with 2D and 1D tensors) + // compute and check results on square-rooted streamflow series std::vector<xt::xarray<double>> metrics = evalhyd::evald<xt::xtensor<double, 2>>( - observed, predicted_2d, {"RMSE", "NSE", "KGE", "KGEPRIME"} + observed, predicted_2d, {"NSE"}, "sqrt" ); - // 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_sqrt = + {{0.882817}, + {0.883023}, + {0.883019}, + {0.883029}, + {0.882972}}; + EXPECT_TRUE(xt::allclose(metrics[0], nse_sqrt)); + + // compute and check results on inverted streamflow series + metrics = + evalhyd::evald<xt::xtensor<double, 2>>( + observed, predicted_2d, {"NSE"}, "inv" + ); + + xt::xtensor<double, 2> nse_inv = + {{0.737323}, + {0.737404}, + {0.737429}, + {0.737546}, + {0.737595}}; + EXPECT_TRUE(xt::allclose(metrics[0], nse_inv)); + + // compute and check results on square-rooted streamflow series + metrics = + evalhyd::evald<xt::xtensor<double, 2>>( + observed, predicted_2d, {"NSE"}, "log" + ); + + xt::xtensor<double, 2> nse_log = + {{0.893344}, + {0.893523}, + {0.893585}, + {0.893758}, + {0.893793}}; + EXPECT_TRUE(xt::allclose(metrics[0], nse_log)); - 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)); }