add "sites" axis to thresholds dimensions

different thresholds may be required of different sites, e.g.
if based on streamflow statistics, which are intrinsically site-specific

include/evalhyd/evalp.hpp

@@ -34,7 +34,7 @@ namespace evalhyd
     ///
     ///    q_thr: ``xt::xtensor<double, 2>``, optional
     ///       Streamflow exceedance threshold(s).
-    ///       shape: (thresholds,)
+    ///       shape: (sites, thresholds)
     ///
     ///    t_msk: ``xt::xtensor<bool, 3>``, optional
     ///       Mask(s) used to generate temporal subsets of the whole streamflow
@@ -62,7 +62,7 @@ namespace evalhyd
     ///       xt::xtensor<double, 4> prd = {{{{ 5.3, 4.2, 5.7, 2.3, 3.1 },
     ///                                       { 4.3, 4.2, 4.7, 4.3, 3.3 },
     ///                                       { 5.3, 5.2, 5.7, 2.3, 3.9 }}}};
-    ///       xt::xtensor<double, 1> thr = { 4.7, 4.3, 5.5, 2.7, 4.1 };
+    ///       xt::xtensor<double, 2> thr = {{ 4.7, 4.3, 5.5, 2.7, 4.1 }};
     ///    
     ///       evalhyd::evalp(obs, prd, {"BS"}, thr);
     ///
@@ -81,7 +81,7 @@ namespace evalhyd
             const xt::xtensor<double, 2>& q_obs,
             const xt::xtensor<double, 4>& q_prd,
             const std::vector<std::string>& metrics,
-            const xt::xtensor<double, 1>& q_thr = {},
+            const xt::xtensor<double, 2>& q_thr = {},
             const xt::xtensor<bool, 3>& t_msk = {}
     )
     {
@@ -98,7 +98,7 @@ namespace evalhyd
         std::size_t n_sit = q_prd.shape(0);
         std::size_t n_ltm = q_prd.shape(1);
         std::size_t n_mbr = q_prd.shape(2);
-        std::size_t n_thr = q_thr.size();
+        std::size_t n_thr = q_thr.shape(1);
         std::size_t n_msk = t_msk.size() < 1 ? 1 : t_msk.shape(0);
 
         // register metrics number of dimensions
@@ -147,10 +147,11 @@ namespace evalhyd
             // instantiate probabilist evaluator
             const auto q_obs_v = xt::view(q_obs, s, xt::all());
             const auto q_prd_v = xt::view(q_prd, s, l, xt::all(), xt::all());
+            const auto q_thr_v = xt::view(q_thr, s, xt::all());
             const auto t_msk_v = xt::view(t_msk, s, xt::all(), xt::all());
 
             eh::probabilist::Evaluator evaluator(
-                    q_obs_v, q_prd_v, q_thr, t_msk_v
+                    q_obs_v, q_prd_v, q_thr_v, t_msk_v
             );
 
             // pre-compute required elt

tests/test_probabilist.cpp

@@ -22,7 +22,7 @@ TEST(ProbabilistTests, TestBrier)
     ifs.close();
 
     // compute scores
-    xt::xtensor<double, 1> thresholds = {690, 534, 445};
+    xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
 
     std::vector<xt::xarray<double>> metrics =
             evalhyd::evalp(
@@ -73,16 +73,13 @@ TEST(ProbabilistTests, TestQuantiles)
     ifs.close();
 
     // compute scores
-    xt::xtensor<double, 1> thresholds = {690, 534, 445};
-
     std::vector<xt::xarray<double>> metrics =
             evalhyd::evalp(
                     // shape: (sites [1], time [t])
                     xt::view(observed, xt::newaxis(), xt::all()),
                     // shape: (sites [1], lead times [1], members [m], time [t])
                     xt::view(predicted, xt::newaxis(), xt::newaxis(), xt::all(), xt::all()),
-                    {"QS", "CRPS"},
-                    thresholds
+                    {"QS", "CRPS"}
             );
 
     // check results
@@ -126,7 +123,7 @@ TEST(ProbabilistTests, TestMasks)
     xt::view(masks, 0, 0, xt::range(0, 20)) = 0;
 
     // compute scores using masks to subset whole record
-    xt::xtensor<double, 1> thresholds = {690, 534, 445};
+    xt::xtensor<double, 2> thresholds = {{690, 534, 445}};
     std::vector<std::string> metrics =
             {"BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS"};