return squeezed xarray

Since the metrics are typically summary statistics, their size is not
very big, so that using xtensor instead of xarray as a data structure is
not as critical as for input data. In turn, using xarray allows for
metrics of different sizes to be returned without unnecessary size 1
dimensions (e.g. when only one threshold is given, or when no temporal
masking is performed). So all output metrics are now returned in their
"natural" shape (e.g. 1D for mono-component metrics, 2D for
multi-component metrics), plus any additional dimension linked to
multi-thresholds, multi-masking, etc.

include/evalhyd/determinist.hpp

@@ -4,6 +4,7 @@
 #include <unordered_map>
 #include <vector>
 #include <xtensor/xexpression.hpp>
+#include <xtensor/xarray.hpp>
 
 #include "utils.hpp"
 #include "determinist/evaluator.hpp"
@@ -26,7 +27,7 @@ namespace evalhyd
     /// \return
     ///     Vector of 1D array of metrics for each time series.
     template <class A>
-    std::vector<A> evald(
+    std::vector<xt::xarray<double>> evald(
             const xt::xexpression<A>& q_obs,
             const xt::xexpression<A>& q_prd,
             const std::vector<std::string>& metrics
@@ -73,7 +74,7 @@ namespace evalhyd
         }
 
         // retrieve or compute requested metrics
-        std::vector<A> r;
+        std::vector<xt::xarray<double>> r;
 
         for ( const auto& metric : metrics )
         {
@@ -82,7 +83,7 @@ namespace evalhyd
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_NSE();
-                r.emplace_back(evaluator.NSE);
+                r.emplace_back(xt::squeeze(evaluator.NSE));
             }
         }
 

include/evalhyd/probabilist.hpp

@@ -5,6 +5,7 @@
 #include <unordered_map>
 #include <vector>
 #include <xtensor/xtensor.hpp>
+#include <xtensor/xarray.hpp>
 #include <xtensor/xview.hpp>
 
 #include "utils.hpp"
@@ -39,7 +40,7 @@ namespace evalhyd
     ///     shape: (1+, time)
     /// \return
     ///     Vector of 2D array of metrics for each threshold.
-    std::vector<xt::xtensor<double, 3>> evalp(
+    std::vector<xt::xarray<double>> evalp(
             const xt::xtensor<double, 2>& q_obs,
             const xt::xtensor<double, 2>& q_prd,
             const std::vector<std::string>& metrics,
@@ -107,7 +108,7 @@ namespace evalhyd
         }
 
         // retrieve or compute requested metrics
-        std::vector<xt::xtensor<double, 3>> r;
+        std::vector<xt::xarray<double>> r;
 
         for (const auto& metric : metrics)
         {
@@ -116,42 +117,42 @@ namespace evalhyd
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_BS();
-                r.emplace_back(evaluator.BS);
+                r.emplace_back(xt::squeeze(evaluator.BS));
             }
             else if ( metric == "BSS" )
             {
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_BSS();
-                r.emplace_back(evaluator.BSS);
+                r.emplace_back(xt::squeeze(evaluator.BSS));
             }
             else if ( metric == "BS_CRD" )
             {
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_BS_CRD();
-                r.emplace_back(evaluator.BS_CRD);
+                r.emplace_back(xt::squeeze(evaluator.BS_CRD));
             }
             else if ( metric == "BS_LBD" )
             {
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_BS_LBD();
-                r.emplace_back(evaluator.BS_LBD);
+                r.emplace_back(xt::squeeze(evaluator.BS_LBD));
             }
             else if ( metric == "QS" )
             {
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_QS();
-                r.emplace_back(evaluator.QS);
+                r.emplace_back(xt::squeeze(evaluator.QS));
             }
             else if ( metric == "CRPS" )
             {
                 if (std::find(req_dep.begin(), req_dep.end(), metric)
                         == req_dep.end())
                     evaluator.calc_CRPS();
-                r.emplace_back(evaluator.CRPS);
+                r.emplace_back(xt::squeeze(evaluator.CRPS));
             }
         }
 

include/evalhyd/probabilist/evaluator.h

@@ -55,12 +55,12 @@ namespace evalhyd
             xt::xtensor<double, 2> crps;
 
             // members for evaluation metrics
-            xt::xtensor<double, 3> BS;
+            xt::xtensor<double, 2> BS;
             xt::xtensor<double, 3> BS_CRD;
             xt::xtensor<double, 3> BS_LBD;
-            xt::xtensor<double, 3> BSS;
-			xt::xtensor<double, 3> QS;
-            xt::xtensor<double, 3> CRPS;
+            xt::xtensor<double, 2> BSS;
+			xt::xtensor<double, 2> QS;
+            xt::xtensor<double, 1> CRPS;
 
             // methods to compute derived data
             void calc_q_qnt();

include/evalhyd/probabilist/evaluator_brier.cpp

@@ -50,7 +50,7 @@ namespace evalhyd
         {
             // initialise output variable
             // shape: (subsets, thresholds, 1)
-            BS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}});
+            BS = xt::zeros<double>({n_msk, n_thr});
 
             // compute variable one mask at a time to minimise memory imprint
             for (int m = 0; m < n_msk; m++)
@@ -61,8 +61,8 @@ namespace evalhyd
 
                 // calculate the mean over the time steps
                 // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
-                xt::view(BS, m, xt::all(), xt::all()) =
-                        xt::nanmean(bs_masked, -1, xt::keep_dims);
+                xt::view(BS, m, xt::all()) =
+                        xt::nanmean(bs_masked, -1);
             }
         }
 
@@ -306,7 +306,7 @@ namespace evalhyd
         {
             // declare and initialise output variable
             // shape: (subsets, thresholds, components)
-            BSS = xt::zeros<double>({n_msk, n_thr, std::size_t {1}});
+            BSS = xt::zeros<double>({n_msk, n_thr});
 
             // compute variable one mask at a time to minimise memory imprint
             for (int m = 0; m < n_msk; m++)
@@ -332,8 +332,8 @@ namespace evalhyd
 
                 // compute Brier skill score(s)
                 // $BSS = \frac{1}{n} \sum_{k=1}^{n} 1 - \frac{bs}{bs_{ref}}
-                xt::view(BSS, m, xt::all(), xt::all()) =
-                        xt::nanmean(1 - (bs_masked / bs_ref), -1, xt::keep_dims);
+                xt::view(BSS, m, xt::all()) =
+                        xt::nanmean(1 - (bs_masked / bs_ref), -1);
             }
         }
     }

include/evalhyd/probabilist/evaluator_quantiles.cpp

@@ -55,7 +55,7 @@ namespace evalhyd
         {
             // initialise output variable
             // shape: (subsets, quantiles, 1)
-            QS = xt::zeros<double>({n_msk, n_mbr, std::size_t {1}});
+            QS = xt::zeros<double>({n_msk, n_mbr});
 
             // compute variable one mask at a time to minimise memory imprint
             for (int m = 0; m < n_msk; m++)
@@ -66,8 +66,8 @@ namespace evalhyd
 
                 // calculate the mean over the time steps
                 // $QS = \frac{1}{n} \sum_{k=1}^{n} qs$
-                xt::view(QS, m, xt::all(), xt::all()) =
-                        xt::nanmean(qs_masked, -1, xt::keep_dims);
+                xt::view(QS, m, xt::all()) =
+                        xt::nanmean(qs_masked, -1);
             }
         }
         
@@ -111,7 +111,7 @@ namespace evalhyd
         {
             // initialise output variable
             // shape: (subsets, thresholds, 1)
-            CRPS = xt::zeros<double>({n_msk, std::size_t {1}, std::size_t {1}});
+            CRPS = xt::zeros<double>({n_msk});
 
             // compute variable one mask at a time to minimise memory imprint
             for (int m = 0; m < n_msk; m++)
@@ -122,8 +122,8 @@ namespace evalhyd
 
                 // calculate the mean over the time steps
                 // $CRPS = \frac{1}{n} \sum_{k=1}^{n} crps$
-                xt::view(CRPS, m, xt::all(), xt::all()) =
-                        xt::nanmean(crps_masked, -1, xt::keep_dims);
+                xt::view(CRPS, m) =
+                        xt::nanmean(crps_masked, -1);
             }
         }
     }

tests/test_determinist.cpp

@@ -3,6 +3,7 @@
 #include <vector>
 #include <gtest/gtest.h>
 #include <xtensor/xtensor.hpp>
+#include <xtensor/xarray.hpp>
 #include <xtensor/xview.hpp>
 #include <xtensor/xmanipulation.hpp>
 #include <xtensor/xcsv.hpp>
@@ -27,19 +28,19 @@ TEST(DeterministTests, TestNSE) {
     xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0);
 
     // compute scores (both with 2D and 1D tensors)
-    std::vector<xt::xtensor<double, 2>> metrics_2d =
+    std::vector<xt::xarray<double>> metrics_2d =
             evalhyd::evald<xt::xtensor<double, 2>>(
                     observed_2d, predicted_2d, {"NSE"}
             );
 
-    std::vector<xt::xtensor<double, 1>> metrics_1d =
+    std::vector<xt::xarray<double>> metrics_1d =
             evalhyd::evald<xt::xtensor<double, 1>>(
                     observed_1d, predicted_1d, {"NSE"}
             );
 
     // check results (both with 2D and 1D tensors)
-    xt::xtensor<double, 2> nse_2d =
-            {{0.71891219}, {0.7190249} , {0.71835777}, {0.71810361}, {0.71776748}};
+    xt::xtensor<double, 1> nse_2d =
+            {0.71891219, 0.7190249, 0.71835777, 0.71810361, 0.71776748};
     EXPECT_TRUE(xt::allclose(metrics_2d[0], nse_2d));
 
     xt::xtensor<double, 1> nse_1d = {0.71891219};

tests/test_probabilist.cpp

@@ -3,6 +3,7 @@
 #include <vector>
 #include <gtest/gtest.h>
 #include <xtensor/xtensor.hpp>
+#include <xtensor/xarray.hpp>
 #include <xtensor/xmanipulation.hpp>
 #include <xtensor/xcsv.hpp>
 
@@ -22,7 +23,7 @@ TEST(ProbabilistTests, TestBrier) {
     // compute scores
     xt::xtensor<double, 1> thresholds = {690, 534, 445};
 
-    std::vector<xt::xtensor<double, 3>> metrics =
+    std::vector<xt::xarray<double>> metrics =
             evalhyd::evalp(
                     observed, predicted, {"BS", "BSS", "BS_CRD", "BS_LBD"},
                     thresholds
@@ -30,30 +31,26 @@ TEST(ProbabilistTests, TestBrier) {
 
     // check results
     // Brier scores
-    xt::xtensor<double, 3> bs =
-            {{{0.10615136},
-              {0.07395622},
-              {0.08669186}}};
+    xt::xtensor<double, 1> bs =
+            {0.10615136, 0.07395622, 0.08669186};
     EXPECT_TRUE(xt::allclose(metrics[0], bs));
 
     // Brier skill scores
-    xt::xtensor<double, 3> bss =
-            {{{0.5705594},
-              {0.6661165},
-              {0.5635126}}};
+    xt::xtensor<double, 1> bss =
+            {0.5705594, 0.6661165, 0.5635126};
     EXPECT_TRUE(xt::allclose(metrics[1], bss));
 
     // Brier calibration-refinement decompositions
-    xt::xtensor<double, 3> bs_crd =
-            {{{0.011411758, 0.1524456, 0.2471852},
-              {0.005532413, 0.1530793, 0.2215031},
-              {0.010139431, 0.1220601, 0.1986125}}};
+    xt::xtensor<double, 2> bs_crd =
+            {{0.011411758, 0.1524456, 0.2471852},
+             {0.005532413, 0.1530793, 0.2215031},
+             {0.010139431, 0.1220601, 0.1986125}};
     EXPECT_TRUE(xt::allclose(metrics[2], bs_crd));
 
     // Brier likelihood-base rate decompositions
-    xt::xtensor<double, 3> bs_lbd =
-            {{{0.012159881, 0.1506234, 0.2446149},
-              {0.008031746, 0.1473869, 0.2133114},
-              {0.017191279, 0.1048221, 0.1743227}}};
+    xt::xtensor<double, 2> bs_lbd =
+            {{0.012159881, 0.1506234, 0.2446149},
+             {0.008031746, 0.1473869, 0.2133114},
+             {0.017191279, 0.1048221, 0.1743227}};
     EXPECT_TRUE(xt::allclose(metrics[3], bs_lbd));
 }