resolves #10

resolves #7

the energy score (ES) is a multi-variate (multi-site)
generalisation of the CRPS

new metric is called CRPS_FROM_ECDF

new metric is called CRPS_FROM_BS, and existing metric is renamed
CRPS_FROM_QS

where data completeness corresponds to the actual number of time steps
available for each subset/sample, once missing data, mask, and bootstrap
have been applied

Low flow events were defined as the complement of high flow events,
meaning that Brier scores were symmetric. But it is more consistent to
include the threshold value in both the definitions of low flow and
high flow events since we either study one or the other, and rarely
both at the same time, especially not using the same threshold values.
Plus, the choice of including the threshold in high flow events rather
than low flow ones was really arbitrary.

Note, only the unit tests for contingency table-based metrics are
impacted because threshold-based metrics are using evants="high".

in `xt::stddev``, when passing a scalar, it is passed to *ddof*
instead of *axes*, so to compute the standard deviation on a given axis,
the axis must be passed in curly braces

because of the random process in the rank computations, and the
different input tensor sizes

replace clumsy xt::sum ... == xt::scalar, by xt::all

while Brier-based scores are symmetric, Contigency Table-based metrics
are not, so that a definition how what an "event" means is required.
A new optional parameter *events* is added to `evalp` (taking either
"high" or "low" as value).

so that specifying the template arguments become optional

Ultimately, the objective is for the user to be able to get the raw
sampled metric values, or the mean and standard deviation of the sampled
metric values, or a series of quantiles of the sampled metric values.
There are still problems with the standard deviation on rtensor, and
the computation of the quantiles does not work on n-dim expressions yet.
So the second and third options are not possible yet, so only the raw
values can be returned. Nonetheless, the machinery and the choice of
where to introduce the summary functionality could be implemented,
which is the purpose of this commit. A new parameter of the bootstrap
experiment called "summary" is added: it can be given a value of 0 (to
get the raw values). In the future, it would also take a value of 1 for
mean+std, and 2 for quantiles.

The bootstrapping method is based on a non-overlapping block sampling
with replacement, where the blocks are years of data. The number of
samples and the sample length (i.e the number of year blocks) are both
customisable.

The method is accessible both for deterministic and probabilistic
evaluation where a new axis is added. For now, the metrics for all the
samples are returned, but in the future, some summary statistics would
be implemented (e.g. quantiles or mean/standard deviation).

/!\ For determinist evaluation, the n-dimensional functionality became
    untenable such that the number of dimensions was fixed and
    restricted to 2D tensors.

New unit tests are included to test both the bootstrapping generator
and the numerical results obtained with the bootstrapping turned on.

An earlier implementation of the masking conditions assumed that the
conditions on streamflow would only be on the observations, but this is
not always the case. For example, reliability scores cannot be done on
the observed streamflow and need to be performed on the predicted
streamflow. So this is now possible as the condition syntax is changed
and now *q_obs*/*q_prd_median*/*q_prd_mean* in place of *q*.

In order to be an accessible parameter type for Python (and hopefully
R in the future) bindings, `std::string` needed to be replaced with
`std::array<char, 32>`, a fixed-length string type.

This functionality is inherited from `evalhyd-cli`. It allows the user
to provide conditions as strings to specify how to generate temporal
subsets. Conditions can be based on observed streamflow values (e.g.
q>800, q<=120) or on time indices (e.g. to select particular events).

This functionality is made available both for determinist and
probabilist evaluation, unlike in  `evalhyd-cli` where it was only
available for probabilist evaluation.

This is documented in the docstrings, and new unit tests are written.

The general approach is to "eliminate" the time steps where observations
or predictions are missing as early as possible in the algorithm. The
best approach seemed to update the user-provided temporal masks to
also mask those time steps with missing data.

An alternative approach would have been to create a view on the
observations and predictions, e.g. using something like
`xt::view(obs, ..., xt::drop(...))`, but this produces a non-contiguous
view which cannot be sorted with `xt::sort` later to determine the
quantiles.

This is documented in `evalp` docstring and new unit tests are added.

resulting in only checking the first metric (i.e. BS) repeatedly

since `xt::allclose` does not have a *equal_nan* like `xt::isclose`
(see https://github.com/xtensor-stack/xtensor/issues/1995), the
check is a bit more convoluted than before...