add pre-computation flow transforms (sqrt, inv, log)

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(

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));
 }