add power function as flow transform

with a customisable exponent through new *exponent* parameter

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(

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