diff --git a/include/evalhyd/detail/determinist/elements.hpp b/include/evalhyd/detail/determinist/elements.hpp
index ab4f8955c59fb7a1bd2d1bbc4c10007f2fb526c1..74007923355a95f16279450a2fed03d9c6c68969 100644
--- a/include/evalhyd/detail/determinist/elements.hpp
+++ b/include/evalhyd/detail/determinist/elements.hpp
@@ -17,26 +17,26 @@ namespace evalhyd
     {
         namespace elements
         {
-            // Compute the mean of the observations.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (1, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Mean observed streamflow.
-            //     shape: (subsets, samples, series, 1)
+            /// Compute the mean of the observations.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (1, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Mean observed streamflow.
+            ///     shape: (subsets, samples, series, 1)
             template <class XD2>
             inline xt::xtensor<double, 4> calc_mean_obs(
                     const XD2& q_obs,
@@ -69,26 +69,26 @@ namespace evalhyd
                 return mean_obs;
             }
 
-            // Compute the mean of the predictions.
-            //
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Mean predicted streamflow.
-            //     shape: (subsets, samples, series, 1)
+            /// Compute the mean of the predictions.
+            ///
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Mean predicted streamflow.
+            ///     shape: (subsets, samples, series, 1)
             template <class XD2>
             inline xt::xtensor<double, 4> calc_mean_prd(
                     const XD2& q_prd,
@@ -121,17 +121,17 @@ namespace evalhyd
                 return mean_prd;
             }
 
-            // Compute the quadratic error between observations and predictions.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (1, time)
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \return
-            //     Quadratic errors between observations and predictions.
-            //     shape: (series, time)
+            /// Compute the quadratic error between observations and predictions.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (1, time)
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \return
+            ///     Quadratic errors between observations and predictions.
+            ///     shape: (series, time)
             template <class XD2>
             inline xt::xtensor<double, 2> calc_quad_err(
                     const XD2& q_obs,
@@ -141,28 +141,28 @@ namespace evalhyd
                 return xt::square(q_obs - q_prd);
             }
 
-            // Compute the quadratic error between observations and mean observation.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (series, time)
-            // \param mean_obs
-            //     Mean observed streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_tim
-            //     Number of time steps.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Quadratic errors between observations and mean observation.
-            //     shape: (subsets, samples, series, time)
+            /// Compute the quadratic error between observations and mean observation.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (series, time)
+            /// \param mean_obs
+            ///     Mean observed streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_tim
+            ///     Number of time steps.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Quadratic errors between observations and mean observation.
+            ///     shape: (subsets, samples, series, time)
             template <class XD2>
             inline xt::xtensor<double, 4> calc_quad_obs(
                     const XD2& q_obs,
@@ -194,28 +194,28 @@ namespace evalhyd
                 return quad_obs;
             }
 
-            // Compute the quadratic error between predictions and mean prediction.
-            //
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \param mean_prd
-            //     Mean predicted streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_tim
-            //     Number of time steps.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Quadratic errors between predictions and mean prediction.
-            //     shape: (subsets, samples, series, time)
+            /// Compute the quadratic error between predictions and mean prediction.
+            ///
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \param mean_prd
+            ///     Mean predicted streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_tim
+            ///     Number of time steps.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Quadratic errors between predictions and mean prediction.
+            ///     shape: (subsets, samples, series, time)
             template <class XD2>
             inline xt::xtensor<double, 4> calc_quad_prd(
                     const XD2& q_prd,
@@ -247,41 +247,41 @@ namespace evalhyd
                 return quad_prd;
             }
 
-            // Compute the Pearson correlation coefficient.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (series, time)
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \param mean_obs
-            //     Mean observed streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param mean_prd
-            //     Mean predicted streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param quad_obs
-            //     Quadratic errors between observations and mean observation.
-            //     shape: (subsets, samples, series, time)
-            // \param quad_prd
-            //     Quadratic errors between predictions and mean prediction.
-            //     shape: (subsets, samples, series, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Pearson correlation coefficients.
-            //     shape: (subsets, samples, series)
+            /// Compute the Pearson correlation coefficient.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (series, time)
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \param mean_obs
+            ///     Mean observed streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param mean_prd
+            ///     Mean predicted streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param quad_obs
+            ///     Quadratic errors between observations and mean observation.
+            ///     shape: (subsets, samples, series, time)
+            /// \param quad_prd
+            ///     Quadratic errors between predictions and mean prediction.
+            ///     shape: (subsets, samples, series, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Pearson correlation coefficients.
+            ///     shape: (subsets, samples, series)
             template <class XD2>
             inline xt::xtensor<double, 3> calc_r_pearson(
                     const XD2& q_obs,
@@ -333,35 +333,35 @@ namespace evalhyd
                 return r_pearson;
             }
 
-            // Compute alpha.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (series, time)
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \param mean_obs
-            //     Mean observed streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param mean_prd
-            //     Mean predicted streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Alphas, ratios of standard deviations.
-            //     shape: (subsets, samples, series)
+            /// Compute alpha.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (series, time)
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \param mean_obs
+            ///     Mean observed streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param mean_prd
+            ///     Mean predicted streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Alphas, ratios of standard deviations.
+            ///     shape: (subsets, samples, series)
             template <class XD2>
             inline xt::xtensor<double, 3> calc_alpha(
                     const XD2& q_obs,
@@ -400,29 +400,29 @@ namespace evalhyd
                 return alpha;
             }
 
-            // Compute the bias.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (series, time)
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (series, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Biases.
-            //     shape: (subsets, samples, series)
+            /// Compute the bias.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (series, time)
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (series, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Biases.
+            ///     shape: (subsets, samples, series)
             template <class XD2>
             inline xt::xtensor<double, 3> calc_bias(
                     const XD2& q_obs,
diff --git a/include/evalhyd/detail/determinist/metrics.hpp b/include/evalhyd/detail/determinist/metrics.hpp
index 5a2d502b255eb4c708344663c6b318ae5c5c83c7..8ee318f0cfadff8f4c4933c79c85f1fbf3c5a6d1 100644
--- a/include/evalhyd/detail/determinist/metrics.hpp
+++ b/include/evalhyd/detail/determinist/metrics.hpp
@@ -16,26 +16,26 @@ namespace evalhyd
     {
         namespace metrics
         {
-            // Compute the root-mean-square error (RMSE).
-            //
-            // \param quad_err
-            //     Quadratic errors between observations and predictions.
-            //     shape: (series, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Root-mean-square errors.
-            //     shape: (series, subsets, samples)
+            /// Compute the root-mean-square error (RMSE).
+            ///
+            /// \param quad_err
+            ///     Quadratic errors between observations and predictions.
+            ///     shape: (series, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Root-mean-square errors.
+            ///     shape: (series, subsets, samples)
             inline xt::xtensor<double, 3> calc_RMSE(
                     const xt::xtensor<double, 2>& quad_err,
                     const xt::xtensor<bool, 3>& t_msk,
@@ -66,29 +66,29 @@ namespace evalhyd
                 return RMSE;
             }
 
-            // Compute the Nash-Sutcliffe Efficiency (NSE).
-            //
-            // \param quad_err
-            //     Quadratic errors between observations and predictions.
-            //     shape: (series, time)
-            // \param quad_obs
-            //     Quadratic errors between observations and mean observation.
-            //     shape: (subsets, samples, series, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (subsets, series, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Nash-Sutcliffe efficiencies.
-            //     shape: (series, subsets, samples)
+            /// Compute the Nash-Sutcliffe Efficiency (NSE).
+            ///
+            /// \param quad_err
+            ///     Quadratic errors between observations and predictions.
+            ///     shape: (series, time)
+            /// \param quad_obs
+            ///     Quadratic errors between observations and mean observation.
+            ///     shape: (subsets, samples, series, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (subsets, series, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Nash-Sutcliffe efficiencies.
+            ///     shape: (series, subsets, samples)
             inline xt::xtensor<double, 3> calc_NSE(
                     const xt::xtensor<double, 2>& quad_err,
                     const xt::xtensor<double, 4>& quad_obs,
@@ -127,26 +127,26 @@ namespace evalhyd
                 return NSE;
             }
 
-            // Compute the Kling-Gupta Efficiency (KGE).
-            //
-            // \param r_pearson
-            //     Pearson correlation coefficients.
-            //     shape: (subsets, samples, series)
-            // \param alpha
-            //     Alphas, ratios of standard deviations.
-            //     shape: (subsets, samples, series)
-            // \param bias
-            //     Biases.
-            //     shape: (subsets, samples, series)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Kling-Gupta efficiencies.
-            //     shape: (series, subsets, samples)
+            /// Compute the Kling-Gupta Efficiency (KGE).
+            ///
+            /// \param r_pearson
+            ///     Pearson correlation coefficients.
+            ///     shape: (subsets, samples, series)
+            /// \param alpha
+            ///     Alphas, ratios of standard deviations.
+            ///     shape: (subsets, samples, series)
+            /// \param bias
+            ///     Biases.
+            ///     shape: (subsets, samples, series)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Kling-Gupta efficiencies.
+            ///     shape: (series, subsets, samples)
             inline xt::xtensor<double, 3> calc_KGE(
                     const xt::xtensor<double, 3>& r_pearson,
                     const xt::xtensor<double, 3>& alpha,
@@ -175,32 +175,32 @@ namespace evalhyd
                 return KGE;
             }
 
-            // Compute the modified Kling-Gupta Efficiency (KGEPRIME).
-            //
-            // \param mean_obs
-            //     Mean observed streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param mean_prd
-            //     Mean predicted streamflow.
-            //     shape: (subsets, samples, series, 1)
-            // \param r_pearson
-            //     Pearson correlation coefficients.
-            //     shape: (subsets, samples, series)
-            // \param alpha
-            //     Alphas, ratios of standard deviations.
-            //     shape: (subsets, samples, series)
-            // \param bias
-            //     Biases.
-            //     shape: (subsets, samples, series)
-            // \param n_srs
-            //     Number of prediction series.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Modified Kling-Gupta efficiencies.
-            //     shape: (series, subsets, samples)
+            /// Compute the modified Kling-Gupta Efficiency (KGEPRIME).
+            ///
+            /// \param mean_obs
+            ///     Mean observed streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param mean_prd
+            ///     Mean predicted streamflow.
+            ///     shape: (subsets, samples, series, 1)
+            /// \param r_pearson
+            ///     Pearson correlation coefficients.
+            ///     shape: (subsets, samples, series)
+            /// \param alpha
+            ///     Alphas, ratios of standard deviations.
+            ///     shape: (subsets, samples, series)
+            /// \param bias
+            ///     Biases.
+            ///     shape: (subsets, samples, series)
+            /// \param n_srs
+            ///     Number of prediction series.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Modified Kling-Gupta efficiencies.
+            ///     shape: (series, subsets, samples)
             inline xt::xtensor<double, 3> calc_KGEPRIME(
                     const xt::xtensor<double, 4>& mean_obs,
                     const xt::xtensor<double, 4>& mean_prd,
diff --git a/include/evalhyd/detail/probabilist/brier.hpp b/include/evalhyd/detail/probabilist/brier.hpp
index 1c6824311004e4281f0a914fa4995bb4fd4eca1e..945fe7532ae7d25c8e864d7372f64a79a2881cf4 100644
--- a/include/evalhyd/detail/probabilist/brier.hpp
+++ b/include/evalhyd/detail/probabilist/brier.hpp
@@ -23,17 +23,17 @@ namespace evalhyd
     {
         namespace elements
         {
-            // Determine observed realisation of threshold(s) exceedance.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (sites, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \return
-            //     Event observed outcome.
-            //     shape: (sites, thresholds, time)
+            /// Determine observed realisation of threshold(s) exceedance.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (sites, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \return
+            ///     Event observed outcome.
+            ///     shape: (sites, thresholds, time)
             template<class XD2>
             inline xt::xtensor<double, 3> calc_o_k(
                     const XD2& q_obs,
@@ -56,30 +56,30 @@ namespace evalhyd
 
             }
 
-            // Determine mean observed realisation of threshold(s) exceedance.
-            //
-            // \param o_k
-            //     Event observed outcome.
-            //     shape: (sites, thresholds, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Mean event observed outcome.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
+            /// Determine mean observed realisation of threshold(s) exceedance.
+            ///
+            /// \param o_k
+            ///     Event observed outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Mean event observed outcome.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
             inline xt::xtensor<double, 5> calc_bar_o(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<bool, 4>& t_msk,
@@ -122,17 +122,17 @@ namespace evalhyd
                 return bar_o;
             }
 
-            // Determine number of forecast members exceeding threshold(s)
-            //
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (sites, lead times, members, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \return
-            //     Number of forecast members exceeding threshold(s).
-            //     shape: (sites, lead times, thresholds, time)
+            /// Determine number of forecast members exceeding threshold(s)
+            ///
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (sites, lead times, members, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \return
+            ///     Number of forecast members exceeding threshold(s).
+            ///     shape: (sites, lead times, thresholds, time)
             template<class XD4, class XD2>
             inline xt::xtensor<double, 4> calc_sum_f_k(
                     const XD4& q_prd,
@@ -164,16 +164,16 @@ namespace evalhyd
                 }
             }
 
-            // Determine forecast probability of threshold(s) exceedance to occur.
-            //
-            // \param sum_f_k
-            //     Number of forecast members exceeding threshold(s).
-            //     shape: (sites, lead times, thresholds, time)
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \return
-            //     Event probability forecast.
-            //     shape: (sites, lead times, thresholds, time)
+            /// Determine forecast probability of threshold(s) exceedance to occur.
+            ///
+            /// \param sum_f_k
+            ///     Number of forecast members exceeding threshold(s).
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \return
+            ///     Event probability forecast.
+            ///     shape: (sites, lead times, thresholds, time)
             inline xt::xtensor<double, 4> calc_y_k(
                     const xt::xtensor<double, 4>& sum_f_k,
                     std::size_t n_mbr
@@ -188,17 +188,17 @@ namespace evalhyd
 
         namespace intermediate
         {
-            // Compute the Brier score for each time step.
-            //
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param y_k
-            //     Event probability forecast.
-            //     shape: (sites, lead times, thresholds, time)
-            // \return
-            //     Brier score for each threshold for each time step.
-            //     shape: (sites, lead times, thresholds, time)
+            /// Compute the Brier score for each time step.
+            ///
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param y_k
+            ///     Event probability forecast.
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \return
+            ///     Brier score for each threshold for each time step.
+            ///     shape: (sites, lead times, thresholds, time)
             inline xt::xtensor<double, 4> calc_bs(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<double, 4>& y_k
@@ -216,33 +216,33 @@ namespace evalhyd
 
         namespace metrics
         {
-            // Compute the Brier score (BS).
-            //
-            // \param bs
-            //     Brier score for each threshold for each time step.
-            //     shape: (sites, lead times, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Brier score for each subset and for each threshold.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
+            /// Compute the Brier score (BS).
+            ///
+            /// \param bs
+            ///     Brier score for each threshold for each time step.
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Brier score for each subset and for each threshold.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 5> calc_BS(
                     const xt::xtensor<double, 4>& bs,
@@ -298,45 +298,45 @@ namespace evalhyd
                 return BS;
             }
 
-            // Compute the calibration-refinement decomposition of the Brier score
-            // into reliability, resolution, and uncertainty.
-            //
-            // BS = reliability - resolution + uncertainty
-            //
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param y_k
-            //     Event probability forecast.
-            //     shape: (sites, lead times, thresholds, time)
-            // \param bar_o
-            //     Mean event observed outcome.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Brier score components (reliability, resolution, uncertainty)
-            //     for each subset and for each threshold.
-            //     shape: (sites, lead times, subsets, samples, thresholds, components)
+            /// Compute the calibration-refinement decomposition of the Brier score
+            /// into reliability, resolution, and uncertainty.
+            ///
+            /// BS = reliability - resolution + uncertainty
+            ///
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param y_k
+            ///     Event probability forecast.
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param bar_o
+            ///     Mean event observed outcome.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Brier score components (reliability, resolution, uncertainty)
+            ///     for each subset and for each threshold.
+            ///     shape: (sites, lead times, subsets, samples, thresholds, components)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_BS_CRD(
                     const XD2& q_thr,
@@ -486,40 +486,40 @@ namespace evalhyd
                 return BS_CRD;
             }
 
-            // Compute the likelihood-base rate decomposition of the Brier score
-            // into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
-            //
-            // BS = type 2 bias - discrimination + sharpness
-            //
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param y_k
-            //     Event probability forecast.
-            //     shape: (sites, lead times, thresholds, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Brier score components (type 2 bias, discrimination, sharpness)
-            //     for each subset and for each threshold.
-            //     shape: (sites, lead times, subsets, samples, thresholds, components)
+            /// Compute the likelihood-base rate decomposition of the Brier score
+            /// into type 2 bias, discrimination, and sharpness (a.k.a. refinement).
+            ///
+            /// BS = type 2 bias - discrimination + sharpness
+            ///
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param y_k
+            ///     Event probability forecast.
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Brier score components (type 2 bias, discrimination, sharpness)
+            ///     for each subset and for each threshold.
+            ///     shape: (sites, lead times, subsets, samples, thresholds, components)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_BS_LBD(
                     const XD2& q_thr,
@@ -679,39 +679,39 @@ namespace evalhyd
                 return BS_LBD;
             }
 
-            // Compute the Brier skill score (BSS).
-            //
-            // \param bs
-            //     Brier score for each threshold for each time step.
-            //     shape: (sites, lead times, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param bar_o
-            //     Mean event observed outcome.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Brier skill score for each subset and for each threshold.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
+            /// Compute the Brier skill score (BSS).
+            ///
+            /// \param bs
+            ///     Brier score for each threshold for each time step.
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param bar_o
+            ///     Mean event observed outcome.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Brier skill score for each subset and for each threshold.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 5> calc_BSS(
                     const xt::xtensor<double, 4>& bs,
diff --git a/include/evalhyd/detail/probabilist/contingency.hpp b/include/evalhyd/detail/probabilist/contingency.hpp
index 6d07a059a4335474c6faaa9734f267b83d1df88c..dbae492e084e5008a4c714742338c9a5fa9ca84a 100644
--- a/include/evalhyd/detail/probabilist/contingency.hpp
+++ b/include/evalhyd/detail/probabilist/contingency.hpp
@@ -35,16 +35,16 @@ namespace evalhyd
             //             +-----+-----+
             //
 
-            // Determine alerts based on forecast.
-            //
-            // \param sum_f_k
-            //     Number of forecast members exceeding threshold(s).
-            //     shape: (sites, lead times, thresholds, time)
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \return
-            //     Alerts based on forecast.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Determine alerts based on forecast.
+            ///
+            /// \param sum_f_k
+            ///     Number of forecast members exceeding threshold(s).
+            ///     shape: (sites, lead times, thresholds, time)
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \return
+            ///     Alerts based on forecast.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_a_k(
                     const xt::xtensor<double, 3>& sum_f_k,
                     std::size_t n_mbr
@@ -61,17 +61,17 @@ namespace evalhyd
                                  xt::all(), xt::newaxis(), xt::newaxis());
             }
 
-            // Determine hits ('a' in contingency table).
-            //
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param a_k
-            //     Alerts based on forecast.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Hits.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Determine hits ('a' in contingency table).
+            ///
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param a_k
+            ///     Alerts based on forecast.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Hits.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_ct_a(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<double, 5>& a_k
@@ -83,17 +83,17 @@ namespace evalhyd
                        && xt::equal(a_k, 1.);
             }
 
-            // Determine false alarms ('b' in contingency table).
-            //
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param a_k
-            //     Alerts based on forecast.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     False alarms.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Determine false alarms ('b' in contingency table).
+            ///
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param a_k
+            ///     Alerts based on forecast.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     False alarms.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_ct_b(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<double, 5>& a_k
@@ -105,17 +105,17 @@ namespace evalhyd
                        && xt::equal(a_k, 1.);
             }
 
-            // Determine misses ('c' in contingency table).
-            //
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param a_k
-            //     Alerts based on forecast.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Misses.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Determine misses ('c' in contingency table).
+            ///
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param a_k
+            ///     Alerts based on forecast.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Misses.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_ct_c(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<double, 5>& a_k
@@ -127,17 +127,17 @@ namespace evalhyd
                        && xt::equal(a_k, 0.);
             }
 
-            // Determine correct rejections ('d' in contingency table).
-            //
-            // \param o_k
-            //     Observed event outcome.
-            //     shape: (sites, thresholds, time)
-            // \param a_k
-            //     Alerts based on forecast.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Correct rejections.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Determine correct rejections ('d' in contingency table).
+            ///
+            /// \param o_k
+            ///     Observed event outcome.
+            ///     shape: (sites, thresholds, time)
+            /// \param a_k
+            ///     Alerts based on forecast.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Correct rejections.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_ct_d(
                     const xt::xtensor<double, 3>& o_k,
                     const xt::xtensor<double, 5>& a_k
@@ -152,17 +152,17 @@ namespace evalhyd
 
         namespace intermediate
         {
-            // Compute the probability of detection for each time step.
-            //
-            // \param ct_a
-            //     Hits.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param ct_c
-            //     Misses.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Probability of detection for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Compute the probability of detection for each time step.
+            ///
+            /// \param ct_a
+            ///     Hits.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param ct_c
+            ///     Misses.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Probability of detection for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_pod(
                     const xt::xtensor<double, 5>& ct_a,
                     const xt::xtensor<double, 5>& ct_c
@@ -171,17 +171,17 @@ namespace evalhyd
                 return ct_a / (ct_a + ct_c);
             }
 
-            // Compute the probability of false detection for each time step.
-            //
-            // \param ct_b
-            //     False alarms.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param ct_d
-            //     Correct rejections.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Probability of false detection for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Compute the probability of false detection for each time step.
+            ///
+            /// \param ct_b
+            ///     False alarms.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param ct_d
+            ///     Correct rejections.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Probability of false detection for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_pofd(
                     const xt::xtensor<double, 5>& ct_b,
                     const xt::xtensor<double, 5>& ct_d
@@ -190,17 +190,17 @@ namespace evalhyd
                 return ct_b / (ct_b + ct_d);
             }
 
-            // Compute the false alarm ratio for each time step.
-            //
-            // \param ct_a
-            //     Hits.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param ct_b
-            //     False alarms.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     False alarm ratio for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Compute the false alarm ratio for each time step.
+            ///
+            /// \param ct_a
+            ///     Hits.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param ct_b
+            ///     False alarms.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     False alarm ratio for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_far(
                     const xt::xtensor<double, 5>& ct_a,
                     const xt::xtensor<double, 5>& ct_b
@@ -209,20 +209,20 @@ namespace evalhyd
                 return ct_b / (ct_a + ct_b);
             }
 
-            // Compute the critical success index for each time step.
-            //
-            // \param ct_a
-            //     Hits.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param ct_b
-            //     False alarms.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param ct_c
-            //     Misses.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \return
-            //     Critical success index for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
+            /// Compute the critical success index for each time step.
+            ///
+            /// \param ct_a
+            ///     Hits.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param ct_b
+            ///     False alarms.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param ct_c
+            ///     Misses.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \return
+            ///     Critical success index for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
             inline xt::xtensor<double, 5> calc_csi(
                     const xt::xtensor<double, 5>& ct_a,
                     const xt::xtensor<double, 5>& ct_b,
@@ -300,36 +300,36 @@ namespace evalhyd
                 }
             }
 
-            // Compute the probability of detection (POD),
-            // also known as 'hit rate'.
-            //
-            // \param pod
-            //     Probability of detection for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Probabilities of detection.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// Compute the probability of detection (POD),
+            /// also known as 'hit rate'.
+            ///
+            /// \param pod
+            ///     Probability of detection for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Probabilities of detection.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_POD(
                     const xt::xtensor<double, 5>& pod,
@@ -350,36 +350,36 @@ namespace evalhyd
                 );
             }
 
-            // Compute the probability of detection (POFD),
-            // also known as 'false alarm rate'.
-            //
-            // \param pofd
-            //     Probability of false detection for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Probabilities of false detection.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// Compute the probability of detection (POFD),
+            /// also known as 'false alarm rate'.
+            ///
+            /// \param pofd
+            ///     Probability of false detection for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Probabilities of false detection.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_POFD(
                     const xt::xtensor<double, 5>& pofd,
@@ -400,35 +400,35 @@ namespace evalhyd
                 );
             }
 
-            // Compute the false alarm ratio (FAR).
-            //
-            // \param far
-            //     False alarm ratio for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     False alarm ratios.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// Compute the false alarm ratio (FAR).
+            ///
+            /// \param far
+            ///     False alarm ratio for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     False alarm ratios.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_FAR(
                     const xt::xtensor<double, 5>& far,
@@ -449,35 +449,35 @@ namespace evalhyd
                 );
             }
 
-            // Compute the critical success index (CSI).
-            //
-            // \param csi
-            //     Critical success index for each time step.
-            //     shape: (sites, lead times, levels, thresholds, time)
-            // \param q_thr
-            //     Streamflow exceedance threshold(s).
-            //     shape: (sites, thresholds)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_thr
-            //     Number of thresholds.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Critical success indices.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// Compute the critical success index (CSI).
+            ///
+            /// \param csi
+            ///     Critical success index for each time step.
+            ///     shape: (sites, lead times, levels, thresholds, time)
+            /// \param q_thr
+            ///     Streamflow exceedance threshold(s).
+            ///     shape: (sites, thresholds)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_thr
+            ///     Number of thresholds.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Critical success indices.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 6> calc_CSI(
                     const xt::xtensor<double, 5>& csi,
@@ -498,17 +498,17 @@ namespace evalhyd
                 );
             }
 
-            // Compute the relative operating characteristic skill score (ROCSS).
-            //
-            // \param POD
-            //     Probabilities of detection.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
-            // \param POFD
-            //     Probabilities of false detection.
-            //     shape: (sites, lead times, subsets, samples, levels, thresholds)
-            // \return
-            //     ROC skill scores.
-            //     shape: (sites, lead times, subsets, samples, thresholds)
+            /// Compute the relative operating characteristic skill score (ROCSS).
+            ///
+            /// \param POD
+            ///     Probabilities of detection.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// \param POFD
+            ///     Probabilities of false detection.
+            ///     shape: (sites, lead times, subsets, samples, levels, thresholds)
+            /// \return
+            ///     ROC skill scores.
+            ///     shape: (sites, lead times, subsets, samples, thresholds)
             template <class XD2>
             inline xt::xtensor<double, 5> calc_ROCSS(
                     const xt::xtensor<double, 6>& POD,
diff --git a/include/evalhyd/detail/probabilist/quantiles.hpp b/include/evalhyd/detail/probabilist/quantiles.hpp
index ca3d668d941c4cb274c3d171c04284e971438f00..471f8b3af6c68bdc0ca051e0c76e53675649bffa 100644
--- a/include/evalhyd/detail/probabilist/quantiles.hpp
+++ b/include/evalhyd/detail/probabilist/quantiles.hpp
@@ -23,14 +23,14 @@ namespace evalhyd
     {
         namespace elements
         {
-            // Compute the forecast quantiles from the ensemble members.
-            //
-            // \param q_prd
-            //     Streamflow predictions.
-            //     shape: (sites, lead times, members, time)
-            // \return
-            //     Streamflow forecast quantiles.
-            //     shape: (sites, lead times, quantiles, time)
+            /// Compute the forecast quantiles from the ensemble members.
+            ///
+            /// \param q_prd
+            ///     Streamflow predictions.
+            ///     shape: (sites, lead times, members, time)
+            /// \return
+            ///     Streamflow forecast quantiles.
+            ///     shape: (sites, lead times, quantiles, time)
             template <class XD4>
             inline xt::xtensor<double, 4> calc_q_qnt(
                     const XD4& q_prd
@@ -42,17 +42,17 @@ namespace evalhyd
 
         namespace intermediate
         {
-            // Compute the quantile scores for each time step.
-            //
-            // \param q_obs
-            //     Streamflow observations.
-            //     shape: (sites, time)
-            // \param q_qnt
-            //     Streamflow quantiles.
-            //     shape: (sites, lead times, quantiles, time)
-            // \return
-            //     Quantile scores for each time step.
-            //     shape: (sites, lead times, quantiles, time)
+            /// Compute the quantile scores for each time step.
+            ///
+            /// \param q_obs
+            ///     Streamflow observations.
+            ///     shape: (sites, time)
+            /// \param q_qnt
+            ///     Streamflow quantiles.
+            ///     shape: (sites, lead times, quantiles, time)
+            /// \return
+            ///     Quantile scores for each time step.
+            ///     shape: (sites, lead times, quantiles, time)
             template <class XD2>
             inline xt::xtensor<double, 4> calc_qs(
                     const XD2 &q_obs,
@@ -82,20 +82,20 @@ namespace evalhyd
                 return qs;
             }
 
-            // Compute the continuous rank probability score(s) based
-            // on quantile scores for each time step, and integrating using the
-            // trapezoidal rule.
-            //
-            // /!\ The number of quantiles must be sufficiently large so that the
-            //     cumulative distribution is smooth enough for the numerical
-            //     integration to be accurate.
-            //
-            // \param qs
-            //     Quantile scores for each time step.
-            //     shape: (sites, lead times, quantiles, time)
-            // \return
-            //     CRPS for each time step.
-            //     shape: (sites, lead times, time)
+            /// Compute the continuous rank probability score(s) based
+            /// on quantile scores for each time step, and integrating using the
+            /// trapezoidal rule.
+            ///
+            /// /!\ The number of quantiles must be sufficiently large so that the
+            ///     cumulative distribution is smooth enough for the numerical
+            ///     integration to be accurate.
+            ///
+            /// \param qs
+            ///     Quantile scores for each time step.
+            ///     shape: (sites, lead times, quantiles, time)
+            /// \return
+            ///     CRPS for each time step.
+            ///     shape: (sites, lead times, time)
             inline xt::xtensor<double, 3> calc_crps(
                     const xt::xtensor<double, 4>& qs,
                     std::size_t n_mbr
@@ -109,30 +109,30 @@ namespace evalhyd
 
         namespace metrics
         {
-            // Compute the quantile score (QS).
-            //
-            // \param qs
-            //     Quantile scores for each time step.
-            //     shape: (sites, lead times, quantiles, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_mbr
-            //     Number of ensemble members.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     Quantile scores.
-            //     shape: (sites, lead times, subsets, samples, quantiles)
+            /// Compute the quantile score (QS).
+            ///
+            /// \param qs
+            ///     Quantile scores for each time step.
+            ///     shape: (sites, lead times, quantiles, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_mbr
+            ///     Number of ensemble members.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     Quantile scores.
+            ///     shape: (sites, lead times, subsets, samples, quantiles)
             inline xt::xtensor<double, 5> calc_QS(
                     const xt::xtensor<double, 4>& qs,
                     const xt::xtensor<bool, 4>& t_msk,
@@ -178,29 +178,29 @@ namespace evalhyd
                 return QS;
             }
 
-            // Compute the continuous rank probability score (CRPS) based
-            // on quantile scores.
-            //
-            // \param crps
-            //     CRPS for each time step.
-            //     shape: (sites, lead times, time)
-            // \param t_msk
-            //     Temporal subsets of the whole record.
-            //     shape: (sites, lead times, subsets, time)
-            // \param b_exp
-            //     Boostrap samples.
-            //     shape: (samples, time slice)
-            // \param n_sit
-            //     Number of sites.
-            // \param n_ldt
-            //     Number of lead times.
-            // \param n_msk
-            //     Number of temporal subsets.
-            // \param n_exp
-            //     Number of bootstrap samples.
-            // \return
-            //     CRPS.
-            //     shape: (sites, lead times, subsets, samples)
+            /// Compute the continuous rank probability score (CRPS) based
+            /// on quantile scores.
+            ///
+            /// \param crps
+            ///     CRPS for each time step.
+            ///     shape: (sites, lead times, time)
+            /// \param t_msk
+            ///     Temporal subsets of the whole record.
+            ///     shape: (sites, lead times, subsets, time)
+            /// \param b_exp
+            ///     Boostrap samples.
+            ///     shape: (samples, time slice)
+            /// \param n_sit
+            ///     Number of sites.
+            /// \param n_ldt
+            ///     Number of lead times.
+            /// \param n_msk
+            ///     Number of temporal subsets.
+            /// \param n_exp
+            ///     Number of bootstrap samples.
+            /// \return
+            ///     CRPS.
+            ///     shape: (sites, lead times, subsets, samples)
             inline xt::xtensor<double, 4> calc_CRPS(
                     const xt::xtensor<double, 3>& crps,
                     const xt::xtensor<bool, 4>& t_msk,