avoid internal copies of pre-computed elements

include/evalhyd/probabilistic/brier.hpp

 #ifndef EVALHYD_BRIER_HPP
 #define EVALHYD_BRIER_HPP
 
+#include <unordered_map>
 #include <xtensor/xtensor.hpp>
 #include <xtensor/xmath.hpp>
 #include <xtensor/xoperation.hpp>
@@ -67,8 +68,6 @@ namespace evalhyd
         // declare internal variables
         // shape: (bins, thresholds, time)
         xt::xtensor<double, 3> msk_bins;
-        // shape: (thresholds, time)
-        xt::xtensor<double, 2> o_k, y_k;
         // shape: (thresholds,)
         xt::xtensor<double, 1> bar_o;
         // shape: (bins, thresholds)
@@ -76,13 +75,9 @@ namespace evalhyd
         // shape: (bins,)
         xt::xtensor<double, 1> y_i;
 
-        // get observed event and forecast probability of event
-        o_k = elt["o_k"];
-        y_k = elt["y_k"];
-
         // define some dimensions
-        std::size_t n = o_k.shape(1);
-        std::size_t n_thr = o_k.shape(0);
+        std::size_t n = elt["o_k"].shape(1);
+        std::size_t n_thr = elt["o_k"].shape(0);
         std::size_t n_cmp = 3;
 
         // declare and initialise output variable
@@ -95,7 +90,7 @@ namespace evalhyd
         // compute mask to subsample time steps belonging to same forecast bin
         // (where bins are defined as the range of forecast values)
         msk_bins = xt::equal(
-                y_k,
+                elt["y_k"],
                 xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
         );
 
@@ -109,7 +104,7 @@ namespace evalhyd
                 xt::sum(
                         xt::where(
                                 msk_bins,
-                                xt::view(o_k, xt::newaxis(), xt::all(), xt::all()),
+                                xt::view(elt["o_k"], xt::newaxis(), xt::all(), xt::all()),
                                 0.
                         ),
                         -1
@@ -119,7 +114,7 @@ namespace evalhyd
 
         // compute mean "climatology" relative frequency of the event
         // $\bar{o} = \frac{1}{n} \sum_{k=1}^{n} o_k$
-        bar_o = xt::mean(o_k, -1);
+        bar_o = xt::mean(elt["o_k"], -1);
 
         // calculate reliability =
         // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
@@ -172,8 +167,6 @@ namespace evalhyd
         // declare internal variables
         // shape: (bins, thresholds, time)
         xt::xtensor<double, 3> msk_bins;
-        // shape: (thresholds, time)
-        xt::xtensor<double, 2> o_k, y_k;
         // shape: (thresholds,)
         xt::xtensor<double, 1> bar_y;
         // shape: (bins, thresholds)
@@ -181,13 +174,9 @@ namespace evalhyd
         // shape: (bins,)
         xt::xtensor<double, 1> o_j;
 
-        // get observed event and forecast probability of event
-        o_k = elt["o_k"];
-        y_k = elt["y_k"];
-
         // define some dimensions
-        std::size_t n = o_k.shape(1);
-        std::size_t n_thr = o_k.shape(0);
+        std::size_t n = elt["o_k"].shape(1);
+        std::size_t n_thr = elt["o_k"].shape(0);
         std::size_t n_cmp = 3;
 
         // declare and initialise output variable
@@ -200,12 +189,12 @@ namespace evalhyd
         // compute mask to subsample time steps belonging to same observation bin
         // (where bins are defined as the range of forecast values)
         msk_bins = xt::equal(
-                o_k,
+                elt["o_k"],
                 xt::view(o_j, xt::all(), xt::newaxis(), xt::newaxis())
         );
 
         // compute number of observations in each observation bin $M_j$
-        M_j= xt::sum(msk_bins, -1);
+        M_j = xt::sum(msk_bins, -1);
 
         // compute subsample relative frequency
         // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$
@@ -214,7 +203,7 @@ namespace evalhyd
                 xt::sum(
                         xt::where(
                                 msk_bins,
-                                xt::view(y_k, xt::newaxis(), xt::all(), xt::all()),
+                                xt::view(elt["y_k"], xt::newaxis(), xt::all(), xt::all()),
                                 0.
                         ),
                         -1
@@ -224,7 +213,7 @@ namespace evalhyd
 
         // compute mean "climatology" forecast probability
         // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
-        bar_y = xt::mean(y_k, -1);
+        bar_y = xt::mean(elt["y_k"], -1);
 
         // calculate type 2 bias =
         // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
@@ -252,7 +241,7 @@ namespace evalhyd
         xt::row(bs_components, 2) =
                 xt::sum(
                         xt::square(
-                                y_k -
+                                elt["y_k"] -
                                 xt::view(bar_y, xt::all(), xt::newaxis())
                         ),
                         -1
@@ -281,17 +270,14 @@ namespace evalhyd
             std::unordered_map<std::string, xt::xtensor<double, 1>>& dep
     )
     {
-        // determine observed realisation of threshold(s) exceedance
-        xt::xtensor<double, 2> o_k = elt["o_k"];
-
         // calculate mean observed event outcome
         // $\bar{o_k} = \frac{1}{n} \sum_{k=1}^{n} o_k$
-        xt::xtensor<double, 2> bar_o_k = xt::mean(o_k, -1, xt::keep_dims);
+        xt::xtensor<double, 2> bar_o_k = xt::mean(elt["o_k"], -1, xt::keep_dims);
 
         // calculate reference Brier score(s)
         // $BS_{ref} = \frac{1}{n} \sum_{k=1}^{n} (o_k - \bar{o_k})^2$
         xt::xtensor<double, 1> bs_ref = xt::mean(
-                xt::square(o_k - bar_o_k), -1
+                xt::square(elt["o_k"] - bar_o_k), -1
         );
 
         // return computed Brier skill score(s)