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.

In order to harmonise with `evalp`, and in view to use `evald` in `evalp`,
the subsets dimension (corresponding to the temporal masks) needed to be
explicitly included in the output of `evald` rather than allowing a
limited number of masks corresponding to the number of time series
evaluated to be used.

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.

there was none for the determinist side

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 use of `xexpression` and template parameter was showing its limits
as the data type could not be different between the different n-dim
inputs (e.g. masks had to be as double rather than bool). Since there is
no real need to be able to use `evald` on `xarray`, the genericity of
its API is reduced so that it now expects inputs to be tensor only. The
number of dimensions for the tensors becomes a template parameter.

Unlike 3de89046, this also deals with
missing predictions. This effectively pairwise remove those time steps
where either observations or predictions are missing data. This does so
without creating copies of the input data, via the use of xfunctions
stored as class members.

The unit test for missing data is updated to include NAN also in the
predictions (with a varying number of NAN across the 5 time series).

In order to avoid copies of the input data, xfunctions are generated
to replace timesteps in predictions where they are NaN in observations.
This means that with this implementation, missing predictions are not
supported.

with a customisable exponent through new *exponent* parameter

It seems that good practice for C++ applications is to only include
*public* headers in "./include" folder and keep source files and
*private* headers and implementation source files in "./src".

Since the metrics are typically summary statistics, their size is not
very big, so that using xtensor instead of xarray as a data structure is
not as critical as for input data. In turn, using xarray allows for
metrics of different sizes to be returned without unnecessary size 1
dimensions (e.g. when only one threshold is given, or when no temporal
masking is performed). So all output metrics are now returned in their
"natural" shape (e.g. 1D for mono-component metrics, 2D for
multi-component metrics), plus any additional dimension linked to
multi-thresholds, multi-masking, etc.

So that the files for the tests look like typical files `evalhyd` would
expect, rather than its transposed version.

This is because in the CLI, metrics is a sequence of unknown length,
so it was only possible to make it last positional parameters. In order
to keep the interfaces harmonised across Python/R/C++ APIs and the CLI,
metrics needed to be put last positional parameter.

i.e. `evalhyd::determinist::evaluate` becomes `evalhyd::evald`, and
`evalhyd::probabilist::evaluate` becomes `evalhyd::evalp`