diff --git a/include/evalhyd/probabilistic.hpp b/include/evalhyd/probabilistic.hpp
index 2cc541d4b584551a00a23b3c0fb9e3903e01fc33..a24163492d40431831b68f01f95e8c733206cfd6 100644
--- a/include/evalhyd/probabilistic.hpp
+++ b/include/evalhyd/probabilistic.hpp
@@ -1,365 +1,119 @@
 #ifndef EVALHYD_PROBABILISTIC_HPP
 #define EVALHYD_PROBABILISTIC_HPP
 
-#include <xtensor/xexpression.hpp>
-#include <xtensor/xmath.hpp>
-#include <xtensor/xoperation.hpp>
-#include <xtensor/xindex_view.hpp>
-#include <xtensor/xio.hpp>
+#include <unordered_map>
+#include <unordered_set>
+#include <xtensor/xtensor.hpp>
 
 #include "utils.hpp"
+#include "probabilistic/elements.hpp"
+#include "probabilistic/brier.hpp"
 
 namespace eh = evalhyd;
 
 namespace evalhyd
 {
-    namespace detail
-    {
-        // determine whether exceedance event has occurred
-        //     q_obs shape: (1, time)
-        //     q_thr shape: (thresholds,)
-        //     returned shape: (thresholds, time)
-        xt::xtensor<double, 2> event_observed_outcome(
+    namespace probabilist {
+        /// Function allowing the evaluation of streamflow forecasts using a
+        /// range of relevant metrics.
+        ///
+        /// \param [in] metrics:
+        ///     Vector of strings for the metric(s) to be computed.
+        /// \param [in] q_obs:
+        ///     2D array of streamflow observations.
+        ///     shape: (1, time)
+        /// \param [in] q_frc:
+        ///     2D array of streamflow forecasts.
+        ///     shape: (members, time)
+        /// \param [in] q_thr:
+        ///     1D array of streamflow exceedance threshold(s).
+        ///     shape: (thresholds,)
+        /// \return
+        ///     Vector of 1D array of metrics for each threshold.
+        std::vector<xt::xtensor<double, 1>> evaluate(
+                const std::vector<std::string>& metrics,
                 const xt::xtensor<double, 2>& q_obs,
-                const xt::xtensor<double, 1>& q_thr
-        )
-        {
-            // determine observed realisation of threshold(s) exceedance
-            return xt::squeeze(eh::utils::is_above_threshold(q_obs, q_thr));
-        }
-
-        // determine forecast probability of event to occur
-        //     q_obs shape: (members, time)
-        //     q_thr shape: (thresholds,)
-        //     returned shape: (thresholds, time)
-        xt::xtensor<double, 2> event_probability_forecast(
                 const xt::xtensor<double, 2>& q_frc,
                 const xt::xtensor<double, 1>& q_thr
         )
         {
-            // determine if members have exceeded threshold(s)
-            xt::xtensor<double, 3> e_frc =
-                    eh::utils::is_above_threshold(q_frc, q_thr);
-
-            // calculate how many members have exceeded threshold(s)
-            xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
-
-            // determine correspondence between number of exceeding members
-            // and forecast probabilities of exceedance
-            // /!\ probability calculation dividing by n (the number of
-            //     members), not n+1 (the number of ranks) like other metrics
-            std::size_t n_mbr = q_frc.shape(0);
-            xt::xtensor<double, 2> p_frc = n_frc / n_mbr;
-
-            return p_frc;
+            // declare maps for memoisation purposes
+            std::unordered_map<std::string, std::vector<std::string>> elt;
+            std::unordered_map<std::string, xt::xtensor<double, 2>> mem_elt;
+            std::unordered_map<std::string, std::vector<std::string>> dep;
+            std::unordered_map<std::string, xt::xtensor<double, 1>> mem_dep;
+
+            // register potentially recurring computation elt across metrics
+            elt["bs"] = {"o_k", "y_k"};
+            elt["bss"] = {"o_k"};
+            elt["bs_crd"] = {"o_k", "y_k"};
+            elt["bs_lbd"] = {"o_k", "y_k"};
+
+            // register nested metrics (i.e. metric dependent on another metric)
+            dep["bss"] = {"bs"};
+
+            // determine required elt/dep to be pre-computed
+            std::unordered_set<std::string> req_elt = {};
+            std::unordered_set<std::string> req_dep = {};
+
+            eh::utils::find_requirements(metrics, elt, dep, req_elt, req_dep);
+
+            // pre-compute required elt
+            for (const auto& element : req_elt)
+            {
+                if ( element == "o_k" )
+                    mem_elt["o_k"] =
+                            eh::elements::event_observed_outcome(q_obs, q_thr);
+                else if ( element == "y_k" )
+                    mem_elt["y_k"] =
+                            eh::elements::event_probability_forecast(q_frc, q_thr);
+            }
+
+            // pre-compute required dep
+            for (const auto& dependency : req_dep)
+            {
+                if ( dependency == "bs" )
+                    mem_dep["bs"] = eh::bs(mem_elt);
+            }
+
+            // retrieve or compute requested metrics
+            std::vector<xt::xtensor<double, 1>> r;
+
+            for (const auto& metric : metrics)
+            {
+                // if exists, retrieve pre-computed metric
+                if (mem_dep.find(metric) != mem_dep.end())
+                {
+                    r.emplace_back(mem_dep[metric]);
+                }
+                    // else, compute relevant metric
+                else if ( metric == "bs" )
+                {
+                    r.emplace_back(eh::bs(mem_elt));
+                }
+                else if ( metric == "bss" )
+                {
+                    r.emplace_back(eh::bss(mem_elt, mem_dep));
+                }
+                else if ( metric == "bs_crd" )
+                {
+                    auto c = eh::bs_crd(mem_elt, q_frc.shape(0));
+                    r.emplace_back(xt::row(c, 0));
+                    r.emplace_back(xt::row(c, 1));
+                    r.emplace_back(xt::row(c, 2));
+                }
+                else if ( metric == "bs_lbd" )
+                {
+                    auto c = eh::bs_lbd(mem_elt);
+                    r.emplace_back(xt::row(c, 0));
+                    r.emplace_back(xt::row(c, 1));
+                    r.emplace_back(xt::row(c, 2));
+                }
+            }
+
+            return r;
         }
     }
-
-    /// Compute the Brier Score (BS).
-    ///
-    /// \param [in] q_obs:
-    ///     2D array of streamflow observations.
-    ///     shape: (1, time)
-    /// \param [in] q_frc:
-    ///     2D array of streamflow forecasts.
-    ///     shape: (members, time)
-    /// \param [in] q_thr:
-    ///     1D array of streamflow exceedance threshold(s).
-    ///     shape: (thresholds,)
-    /// \return
-    ///     1D array of Brier Score for each threshold.
-    ///     shape: (thresholds,)
-    xt::xtensor<double, 1> bs(
-            const xt::xtensor<double, 2>& q_obs,
-            const xt::xtensor<double, 2>& q_frc,
-            const xt::xtensor<double, 1>& q_thr
-    )
-    {
-        // get observed event and forecast probability of event
-        xt::xtensor<double, 2> p_obs =
-                detail::event_observed_outcome(q_obs, q_thr);
-
-        xt::xtensor<double, 2> p_frc =
-                detail::event_probability_forecast(q_frc, q_thr);
-
-        // return computed BS
-        return xt::mean(xt::square(p_obs - p_frc), -1);
-    }
-
-    /// Compute the calibration-refinement decomposition of the Brier Score
-    /// into reliability, resolution, and uncertainty.
-    ///
-    /// BS = reliability - resolution + uncertainty
-    ///
-    /// \param [in] q_obs:
-    ///     2D array of streamflow observations.
-    ///     shape: (1, time)
-    /// \param [in] q_frc:
-    ///     2D array of streamflow forecasts.
-    ///     shape: (members, time)
-    /// \param [in] q_thr:
-    ///     1D array of streamflow exceedance threshold(s).
-    ///     shape: (thresholds,)
-    /// \return
-    ///     2D array of Brier Score components (reliability, resolution,
-    ///     uncertainty) for each threshold.
-    ///     shape: (components, thresholds)
-    xt::xtensor<double, 2> bs_crd(
-            const xt::xtensor<double, 2>& q_obs,
-            const xt::xtensor<double, 2>& q_frc,
-            const xt::xtensor<double, 1>& q_thr
-    )
-    {
-        // NOTE ----------------------------------------------------------------
-        // All equations are following notations from "Wilks, D. S. (2011).
-        // Statistical methods in the atmospheric sciences. Amsterdam; Boston:
-        // Elsevier Academic Press. ISBN: 9780123850225". In particular,
-        // pp. 302-303, 332-333.
-        // ---------------------------------------------------------------------
-
-        // define some dimensions
-        std::size_t n = q_frc.shape(1);
-        std::size_t n_mbr = q_frc.shape(0);
-        std::size_t n_thr = q_thr.size();
-        std::size_t n_cmp = 3;
-
-        // 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)
-        xt::xtensor<double, 2> N_i, bar_o_i;
-        // shape: (bins,)
-        xt::xtensor<double, 1> y_i;
-
-        // declare and initialise output variable
-        // shape: (components, thresholds)
-        xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
-
-        // get observed event and forecast probability of event
-        o_k = detail::event_observed_outcome(q_obs, q_thr);
-        y_k = detail::event_probability_forecast(q_frc, q_thr);
-
-        // compute range of forecast values $y_i$
-        y_i = xt::arange<double>(double (n_mbr + 1)) / n_mbr;
-
-        // 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,
-                xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
-        );
-
-        // compute number of forecasts in each forecast bin $N_i$
-        N_i = xt::sum(msk_bins, -1);
-
-        // compute subsample relative frequency
-        // $\bar{o_i} = \frac{1}{N_i} \sum_{k \in N_i} o_k$
-        bar_o_i = xt::where(
-                N_i > 0,
-                xt::sum(
-                        xt::where(
-                                msk_bins,
-                                xt::view(o_k, xt::newaxis(), xt::all(), xt::all()),
-                                0.
-                        ),
-                        -1
-                ) / N_i,
-                0.
-        );
-
-        // 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);
-
-        // calculate reliability =
-        // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
-        xt::row(bs_components, 0) =
-                xt::sum(
-                        xt::square(
-                                xt::view(y_i, xt::all(), xt::newaxis())
-                                - bar_o_i
-                        ) * N_i,
-                        0
-                ) / n;
-
-        // calculate resolution =
-        // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
-        xt::row(bs_components, 1) =
-                xt::sum(
-                        xt::square(
-                                bar_o_i - bar_o
-                        ) * N_i,
-                        0
-                ) / n;
-
-        // calculate uncertainty = $\bar{o} (1 - \bar{o})$
-        xt::row(bs_components, 2) = bar_o * (1 - bar_o);
-
-        return bs_components;
-    }
-
-    /// Compute the likelihood-based decomposition of the Brier Score
-    /// into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
-    ///
-    /// BS = type 2 bias - discrimination + sharpness
-    ///
-    /// \param [in] q_obs:
-    ///     2D array of streamflow observations.
-    ///     shape: (1, time)
-    /// \param [in] q_frc:
-    ///     2D array of streamflow forecasts.
-    ///     shape: (members, time)
-    /// \param [in] q_thr:
-    ///     1D array of streamflow exceedance threshold(s).
-    ///     shape: (thresholds,)
-    /// \return
-    ///     2D array of Brier Score components (type 2 bias, discrimination,
-    ///     sharpness) for each threshold.
-    ///     shape: (components, thresholds)
-    xt::xtensor<double, 2> bs_lbd(
-            const xt::xtensor<double, 2>& q_obs,
-            const xt::xtensor<double, 2>& q_frc,
-            const xt::xtensor<double, 1>& q_thr
-    )
-    {
-        // NOTE ----------------------------------------------------------------
-        // All equations are following notations from "Wilks, D. S. (2011).
-        // Statistical methods in the atmospheric sciences. Amsterdam; Boston:
-        // Elsevier Academic Press. ISBN: 9780123850225". In particular,
-        // pp. 302-303, 332-333.
-        // ---------------------------------------------------------------------
-
-        // define some dimensions
-        std::size_t n = q_frc.shape(1);
-        std::size_t n_thr = q_thr.size();
-        std::size_t n_cmp = 3;
-
-        // 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)
-        xt::xtensor<double, 2> M_j, bar_y_j;
-        // shape: (bins,)
-        xt::xtensor<double, 1> o_j;
-
-        // declare and initialise output variable
-        // shape: (components, thresholds)
-        xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
-
-        // get observed event and forecast probability of event
-        o_k = detail::event_observed_outcome(q_obs, q_thr);
-        y_k = detail::event_probability_forecast(q_frc, q_thr);
-
-        // set the range of observed values $o_j$
-        o_j = {0., 1.};
-
-        // 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,
-                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);
-
-        // compute subsample relative frequency
-        // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$
-        bar_y_j = xt::where(
-                M_j > 0,
-                xt::sum(
-                        xt::where(
-                                msk_bins,
-                                xt::view(y_k, xt::newaxis(), xt::all(), xt::all()),
-                                0.
-                        ),
-                        -1
-                ) / M_j,
-                0.
-        );
-
-        // compute mean "climatology" forecast probability
-        // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
-        bar_y = xt::mean(y_k, -1);
-
-        // calculate type 2 bias =
-        // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
-        xt::row(bs_components, 0) =
-                xt::sum(
-                        xt::square(
-                                xt::view(o_j, xt::all(), xt::newaxis())
-                                - bar_y_j
-                        ) * M_j,
-                        0
-                ) / n;
-
-        // calculate discrimination =
-        // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
-        xt::row(bs_components, 1) =
-                xt::sum(
-                        xt::square(
-                                bar_y_j - bar_y
-                        ) * M_j,
-                        0
-                ) / n;
-
-        // calculate sharpness =
-        // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
-        xt::row(bs_components, 2) =
-                xt::sum(
-                        xt::square(
-                                y_k -
-                                xt::view(bar_y, xt::all(), xt::newaxis())
-                        ),
-                        -1
-                ) / n;
-
-        return bs_components;
-    }
-
-    /// Compute the Brier Skill Score (BSS).
-    ///
-    /// \param [in] q_obs:
-    ///     2D array of streamflow observations.
-    ///     shape: (1, time)
-    /// \param [in] q_frc:
-    ///     2D array of streamflow forecasts.
-    ///     shape: (members, time)
-    /// \param [in] q_thr:
-    ///     1D array of streamflow exceedance threshold(s).
-    ///     shape: (thresholds,)
-    /// \return
-    ///     1D array of Brier Skill Score for each threshold.
-    ///     shape: (thresholds,)
-    xt::xtensor<double, 1> bss(
-            const xt::xtensor<double, 2>& q_obs,
-            const xt::xtensor<double, 2>& q_frc,
-            const xt::xtensor<double, 1>& q_thr
-    )
-    {
-        // determine observed realisation of threshold(s) exceedance
-        xt::xtensor<double, 2> p_obs =
-                detail::event_observed_outcome(q_obs, q_thr);
-
-        // calculate mean observed realisation of threshold(s) exceedance
-        xt::xtensor<double, 2> p_obs_mean = xt::mean(p_obs, -1, xt::keep_dims);
-
-        // calculate reference Brier Score
-        xt::xtensor<double, 1> bs_ref = xt::mean(
-                xt::square(p_obs - p_obs_mean), -1
-        );
-
-        // return computed BSS
-        return 1 - (bs(q_obs, q_frc, q_thr) / bs_ref);
-    }
 }
 
 #endif //EVALHYD_PROBABILISTIC_HPP
diff --git a/include/evalhyd/probabilistic/brier.hpp b/include/evalhyd/probabilistic/brier.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ade5fbfa3ed293ffa694b52894c8df6009ca8471
--- /dev/null
+++ b/include/evalhyd/probabilistic/brier.hpp
@@ -0,0 +1,307 @@
+#ifndef EVALHYD_BRIER_HPP
+#define EVALHYD_BRIER_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xoperation.hpp>
+
+#include "elements.hpp"
+
+namespace eh = evalhyd;
+
+namespace evalhyd
+{
+    // NOTE --------------------------------------------------------------------
+    // All equations in metrics below are following notations from
+    // "Wilks, D. S. (2011). Statistical methods in the atmospheric sciences.
+    // Amsterdam; Boston: Elsevier Academic Press. ISBN: 9780123850225".
+    // In particular, pp. 302-303, 332-333.
+    // -------------------------------------------------------------------------
+
+    /// Compute the Brier score (BS).
+    ///
+    /// \param [in] elt:
+    ///     Map of pre-computed elements, including required elements:
+    ///         o_k:
+    ///             Observed event outcome.
+    ///             shape: (thresholds, time)
+    ///         y_k:
+    ///             Event probability forecast.
+    ///             shape: (thresholds, time)
+    /// \return
+    ///     1D array of Brier Score for each threshold.
+    ///     shape: (thresholds,)
+    xt::xtensor<double, 1> bs(
+            std::unordered_map<std::string, xt::xtensor<double, 2>>& elt
+    )
+    {
+        // return computed Brier score(s)
+        // $BS = \frac{1}{n} \sum_{k=1}^{n} (o_k - y_k)^2$
+        return xt::mean(xt::square(elt["o_k"] - elt["y_k"]), -1);
+    }
+
+    /// Compute the calibration-refinement decomposition of the Brier score
+    /// into reliability, resolution, and uncertainty.
+    ///
+    /// BS = reliability - resolution + uncertainty
+    ///
+    /// \param [in] elt:
+    ///     Map of pre-computed elements, including required elements:
+    ///         o_k:
+    ///             Observed event outcome.
+    ///             shape: (thresholds, time)
+    ///         y_k:
+    ///             Event probability forecast.
+    ///             shape: (thresholds, time)
+    /// \param [in] n_mbr:
+    ///     Number of forecast members.
+    /// \return
+    ///     2D array of Brier score components (reliability, resolution,
+    ///     uncertainty) for each threshold.
+    ///     shape: (components, thresholds)
+    xt::xtensor<double, 2> bs_crd(
+            std::unordered_map<std::string, xt::xtensor<double, 2>>& elt,
+            std::size_t n_mbr
+    )
+    {
+        // 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)
+        xt::xtensor<double, 2> N_i, bar_o_i;
+        // 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_cmp = 3;
+
+        // declare and initialise output variable
+        // shape: (components, thresholds)
+        xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
+
+        // compute range of forecast values $y_i$
+        y_i = xt::arange<double>(double (n_mbr + 1)) / n_mbr;
+
+        // 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,
+                xt::view(y_i, xt::all(), xt::newaxis(), xt::newaxis())
+        );
+
+        // compute number of forecasts in each forecast bin $N_i$
+        N_i = xt::sum(msk_bins, -1);
+
+        // compute subsample relative frequency
+        // $\bar{o_i} = \frac{1}{N_i} \sum_{k \in N_i} o_k$
+        bar_o_i = xt::where(
+                N_i > 0,
+                xt::sum(
+                        xt::where(
+                                msk_bins,
+                                xt::view(o_k, xt::newaxis(), xt::all(), xt::all()),
+                                0.
+                        ),
+                        -1
+                ) / N_i,
+                0.
+        );
+
+        // 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);
+
+        // calculate reliability =
+        // $\frac{1}{n} \sum_{i=1}^{I} N_i (y_i - \bar{o_i})^2$
+        xt::row(bs_components, 0) =
+                xt::sum(
+                        xt::square(
+                                xt::view(y_i, xt::all(), xt::newaxis())
+                                - bar_o_i
+                        ) * N_i,
+                        0
+                ) / n;
+
+        // calculate resolution =
+        // $\frac{1}{n} \sum_{i=1}^{I} N_i (\bar{o_i} - \bar{o})^2$
+        xt::row(bs_components, 1) =
+                xt::sum(
+                        xt::square(
+                                bar_o_i - bar_o
+                        ) * N_i,
+                        0
+                ) / n;
+
+        // calculate uncertainty = $\bar{o} (1 - \bar{o})$
+        xt::row(bs_components, 2) = bar_o * (1 - bar_o);
+
+        return bs_components;
+    }
+
+    /// Compute the likelihood-based decomposition of the Brier score
+    /// into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
+    ///
+    /// BS = type 2 bias - discrimination + sharpness
+    ///
+    /// \param [in] elt:
+    ///     Map of pre-computed elements, including required elements:
+    ///         o_k:
+    ///             Observed event outcome.
+    ///             shape: (thresholds, time)
+    ///         y_k:
+    ///             Event probability forecast.
+    ///             shape: (thresholds, time)
+    /// \return
+    ///     2D array of Brier score components (type 2 bias, discrimination,
+    ///     sharpness) for each threshold.
+    ///     shape: (components, thresholds)
+    xt::xtensor<double, 2> bs_lbd(
+            std::unordered_map<std::string, xt::xtensor<double, 2>>& elt
+    )
+    {
+        // 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)
+        xt::xtensor<double, 2> M_j, bar_y_j;
+        // 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_cmp = 3;
+
+        // declare and initialise output variable
+        // shape: (components, thresholds)
+        xt::xtensor<double, 2> bs_components = xt::zeros<double>({n_cmp, n_thr});
+
+        // get observed event and forecast probability of event
+        o_k = elt["o_k"];
+        y_k = elt["y_k"];
+
+        // set the range of observed values $o_j$
+        o_j = {0., 1.};
+
+        // 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,
+                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);
+
+        // compute subsample relative frequency
+        // $\bar{y_j} = \frac{1}{M_j} \sum_{k \in M_j} y_k$
+        bar_y_j = xt::where(
+                M_j > 0,
+                xt::sum(
+                        xt::where(
+                                msk_bins,
+                                xt::view(y_k, xt::newaxis(), xt::all(), xt::all()),
+                                0.
+                        ),
+                        -1
+                ) / M_j,
+                0.
+        );
+
+        // compute mean "climatology" forecast probability
+        // $\bar{y} = \frac{1}{n} \sum_{k=1}^{n} y_k$
+        bar_y = xt::mean(y_k, -1);
+
+        // calculate type 2 bias =
+        // $\frac{1}{n} \sum_{j=1}^{J} M_j (o_j - \bar{y_j})^2$
+        xt::row(bs_components, 0) =
+                xt::sum(
+                        xt::square(
+                                xt::view(o_j, xt::all(), xt::newaxis())
+                                - bar_y_j
+                        ) * M_j,
+                        0
+                ) / n;
+
+        // calculate discrimination =
+        // $\frac{1}{n} \sum_{j=1}^{J} M_j (\bar{y_j} - \bar{y})^2$
+        xt::row(bs_components, 1) =
+                xt::sum(
+                        xt::square(
+                                bar_y_j - bar_y
+                        ) * M_j,
+                        0
+                ) / n;
+
+        // calculate sharpness =
+        // $\frac{1}{n} \sum_{k=1}^{n} (\bar{y_k} - \bar{y})^2$
+        xt::row(bs_components, 2) =
+                xt::sum(
+                        xt::square(
+                                y_k -
+                                xt::view(bar_y, xt::all(), xt::newaxis())
+                        ),
+                        -1
+                ) / n;
+
+        return bs_components;
+    }
+
+    /// Compute the Brier skill score (BSS).
+    ///
+    /// \param [in] elt:
+    ///     Map of pre-computed elements, including required elements:
+    ///         o_k:
+    ///             Observed event outcome.
+    ///             shape: (thresholds, time)
+    /// \param [in] dep:
+    ///     Map of pre-computed metrics, including required metrics:
+    ///         bs:
+    ///             Brier score(s).
+    ///             shape: (thresholds,)
+    /// \return
+    ///     1D array of Brier skill score for each threshold.
+    ///     shape: (thresholds,)
+    xt::xtensor<double, 1> bss(
+            std::unordered_map<std::string, xt::xtensor<double, 2>>& elt,
+            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);
+
+        // 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
+        );
+
+        // return computed Brier skill score(s)
+        // $BSS = 1 - \frac{BS}{BS_{ref}}
+        return 1 - (dep["bs"] / bs_ref);
+    }
+}
+
+#endif //EVALHYD_BRIER_HPP
diff --git a/include/evalhyd/probabilistic/elements.hpp b/include/evalhyd/probabilistic/elements.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9253f9018e51244f3319f402ec0294f99b5bb21c
--- /dev/null
+++ b/include/evalhyd/probabilistic/elements.hpp
@@ -0,0 +1,99 @@
+#ifndef EVALHYD_ELEMENTS_HPP
+#define EVALHYD_ELEMENTS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include "xtensor/xview.hpp"
+#include <xtensor/xmath.hpp>
+
+
+namespace eh = evalhyd;
+
+namespace evalhyd
+{
+    namespace elements
+    {
+        namespace detail
+        {
+            // determine whether flows `q` are greater than given threshold(s) `thr`
+            //     q shape: (1+, time)
+            //     thr shape: (thresholds,)
+            //     returned shape: (thresholds, 1+, time)
+            xt::xtensor<double, 3> is_above_threshold(
+                    const xt::xtensor<double, 2>& q,
+                    const xt::xtensor<double, 1>& thr
+            )
+            {
+                // return a boolean-like array
+                return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
+            }
+
+            // determine whether flows `q` are strictly lower than given threshold(s) `thr`
+            //     q shape: (1+, time)
+            //     thr shape: (thresholds,)
+            //     returned shape: (thresholds, 1+, time)
+            xt::xtensor<double, 3> is_below_threshold(
+                    const xt::xtensor<double, 2>& q,
+                    const xt::xtensor<double, 1>& thr
+            )
+            {
+                // return a boolean-like array
+                return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
+            }
+        }
+
+        /// Determine whether exceedance event has occurred.
+        ///
+        /// \param [in] q_obs:
+        ///     2D array of streamflow observations.
+        ///     shape: (1, time)
+        /// \param [in] q_thr:
+        ///     1D array of streamflow exceedance threshold(s).
+        ///     shape: (thresholds,)
+        /// \return
+        ///     2D array of event observed outcome
+        ///     shape: (thresholds, time)
+        xt::xtensor<double, 2> event_observed_outcome(
+                const xt::xtensor<double, 2>& q_obs,
+                const xt::xtensor<double, 1>& q_thr
+        )
+        {
+            // determine observed realisation of threshold(s) exceedance
+            return xt::squeeze(detail::is_above_threshold(q_obs, q_thr));
+        }
+
+        /// Determine forecast probability of event to occur.
+        ///
+        /// \param [in] q_frc:
+        ///     2D array of streamflow forecasts.
+        ///     shape: (members, time)
+        /// \param [in] q_thr:
+        ///     1D array of streamflow exceedance threshold(s).
+        ///     shape: (thresholds,)
+        /// \return
+        ///     2D array of event probability forecast
+        ///     shape: (thresholds, time)
+        xt::xtensor<double, 2> event_probability_forecast(
+                const xt::xtensor<double, 2>& q_frc,
+                const xt::xtensor<double, 1>& q_thr
+        )
+        {
+            // determine if members have exceeded threshold(s)
+            xt::xtensor<double, 3> e_frc =
+                    detail::is_above_threshold(q_frc, q_thr);
+
+            // calculate how many members have exceeded threshold(s)
+            xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
+
+            // determine correspondence between number of exceeding members
+            // and forecast probabilities of exceedance
+            // /!\ probability calculation dividing by n (the number of
+            //     members), not n+1 (the number of ranks) like in other metrics
+            std::size_t n_mbr = q_frc.shape(0);
+            xt::xtensor<double, 2> p_frc = n_frc / n_mbr;
+
+            return p_frc;
+        }
+    }
+}
+
+#endif //EVALHYD_ELEMENTS_HPP
diff --git a/include/evalhyd/utils.hpp b/include/evalhyd/utils.hpp
index d8bbb815fc050877bfcb072772cd4050e775fb17..a2c8b9cd78d036bb491082b0f79b108707271036 100644
--- a/include/evalhyd/utils.hpp
+++ b/include/evalhyd/utils.hpp
@@ -1,40 +1,53 @@
 #ifndef EVALHYD_UTILS_HPP
 #define EVALHYD_UTILS_HPP
 
+#include <unordered_map>
+#include <unordered_set>
 #include <xtensor/xtensor.hpp>
-#include <xtensor/xbroadcast.hpp>
-#include <xtensor/xadapt.hpp>
-#include <xtensor/xview.hpp>
-#include <xtensor/xmath.hpp>
 
 namespace evalhyd
 {
     namespace utils
     {
-        // determine whether flows `q` are greater than given threshold(s) `thr`
-        //     q shape: (1+, time)
-        //     thr shape: (thresholds,)
-        //     returned shape: (thresholds, 1+, time)
-        xt::xtensor<double, 3> is_above_threshold(
-                const xt::xtensor<double, 2>& q,
-                const xt::xtensor<double, 1>& thr
+        /// Procedure to determine based on a list of metrics which elements
+        /// and which metrics (and their associated elements) require to be
+        /// pre-computed for memoisation purposes.
+        ///
+        /// \param [in] metrics:
+        ///     Vector of strings for the metric(s) to be computed.
+        /// \param [in] elements:
+        ///     Map between metrics and their required computation elements.
+        /// \param [in] dependencies:
+        ///     Map between metrics and their dependencies.
+        /// \param [out] required_elements:
+        ///     Set of elements identified as required to be pre-computed.
+        /// \param [out] required_dependencies:
+        ///     Set of metrics identified as required to be pre-computed.
+        void find_requirements (
+                const std::vector<std::string>& metrics,
+                std::unordered_map<std::string, std::vector<std::string>>& elements,
+                std::unordered_map<std::string, std::vector<std::string>>& dependencies,
+                std::unordered_set<std::string>& required_elements,
+                std::unordered_set<std::string>& required_dependencies
         )
         {
-            // return a boolean-like array
-            return q >= xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
-        }
-
-        // determine whether flows `q` are strictly lower than given threshold(s) `thr`
-        //     q shape: (1+, time)
-        //     thr shape: (thresholds,)
-        //     returned shape: (thresholds, 1+, time)
-        xt::xtensor<double, 3> is_below_threshold(
-                const xt::xtensor<double, 2>& q,
-                const xt::xtensor<double, 1>& thr
-        )
-        {
-            // return a boolean-like array
-            return q < xt::view(thr, xt::all(), xt::newaxis(), xt::newaxis());
+            for (const auto& metric : metrics)
+            {
+                // add elements to pre-computation set
+                for (const auto& element : elements[metric])
+                    required_elements.insert(element);
+                // add metric dependencies to pre-computation set
+                if (dependencies.find(metric) != dependencies.end())
+                {
+                    for (const auto& dependency : dependencies[metric])
+                    {
+                        required_dependencies.insert(dependency);
+                        // add dependency elements to pre-computation set
+                        for (const auto& element : elements[dependency])
+                            required_elements.insert(element);
+                    }
+                }
+            }
         }
     }
 }