From 0336dbe6bb59ffa46bba46fcaf70f68a22d9bdc7 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Tue, 20 Dec 2022 16:16:37 -0500
Subject: [PATCH 01/26] Add functions to partition uncertainties according to
 method from Hawkins & Sutton (2009).

---
 xclim/ensembles/__init__.py              |   2 +
 xclim/ensembles/_filters.py              | 147 ++++++++++++++++++
 xclim/ensembles/_partitioning.py         | 186 +++++++++++++++++++++++
 xclim/testing/tests/test_partitioning.py |  18 +++
 4 files changed, 353 insertions(+)
 create mode 100644 xclim/ensembles/_filters.py
 create mode 100644 xclim/ensembles/_partitioning.py
 create mode 100644 xclim/testing/tests/test_partitioning.py

diff --git a/xclim/ensembles/__init__.py b/xclim/ensembles/__init__.py
index d77ba220a..440aeb16a 100644
--- a/xclim/ensembles/__init__.py
+++ b/xclim/ensembles/__init__.py
@@ -8,5 +8,7 @@
 from __future__ import annotations
 
 from ._base import create_ensemble, ensemble_mean_std_max_min, ensemble_percentiles
+from ._filters import concat_hist, model_in_all_scens, single_member
+from ._partitioning import hawkins_sutton
 from ._reduce import kkz_reduce_ensemble, kmeans_reduce_ensemble, plot_rsqprofile
 from ._robustness import change_significance, robustness_coefficient
diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
new file mode 100644
index 000000000..ff88f0cf2
--- /dev/null
+++ b/xclim/ensembles/_filters.py
@@ -0,0 +1,147 @@
+import numpy as np
+import xarray as xr
+
+
+def concat_hist(da, **hist):
+    """Concatenate historical scenario with future scenarios along time.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Input data where the historical scenario is stored alongside other, future, scenarios.
+    hist: {str: str}
+      Mapping of the scenario dimension name to the historical scenario coordinate, e.g. `scen="historical"`.
+
+    Returns
+    -------
+    xr.DataArray
+      Data with the historical scenario is stacked in time before each one of the other scenarios.
+
+    Notes
+    -----
+    Data goes from:
+
+        +----------+----------------------------------+
+        | Scenario | Time                             +
+        +==========+==================================+
+        | hist     | hhhhhhhhhhhhhhhh                 |
+        +----------+----------------------------------+
+        | scen1    |                 111111111111     |
+        +----------+----------------------------------+
+        | scen2    |                 222222222222     |
+        +----------+----------------------------------+
+
+    to:
+        +----------+----------------------------------+
+        | Scenario | Time                             +
+        +==========+==================================+
+        | scen1    | hhhhhhhhhhhhhhhh111111111111     |
+        +----------+----------------------------------+
+        | scen2    | hhhhhhhhhhhhhhhh222222222222     |
+        +----------+----------------------------------+
+    """
+    if len(hist) > 1:
+        raise ValueError("Too many values in hist scenario.")
+
+    ((dim, name),) = hist.items()
+
+    h = da.sel(**hist).dropna("time", how="all")
+    ens = da.drop_sel(**hist)
+
+    index = ens[dim]
+    bare = ens.drop(dim).dropna("time", how="all")
+
+    return xr.concat([h, bare], dim="time").assign_coords({dim: index})
+
+
+def model_in_all_scens(da, dimensions=None):
+    """Return data with only simulations that have at least one member in each scenario.
+
+    Parameters
+    ----------
+    da: xr.DataArray
+      Input data with dimensions for time, member, model and scenario.
+    dimensions: dict
+      Mapping from original dimension names to standard dimension names: scenario, model, member.
+
+    Returns
+    -------
+    xr.DataArray
+      Data for models that have values for all scenarios.
+
+    Notes
+    -----
+    In the following example, `GCM_C` would be filtered out from the data because it has no member for `scen2`.
+
+    +-------+-------+-------+
+    | Model | Members       |
+    +-------+---------------+
+    |       | scen1 | scen2 |
+    +=======+=======+=======+
+    | GCM_A | 1,2,3 | 1,2,3 |
+    +-------+-------+-------+
+    | GCM_B | 1     | 2,3   |
+    +-------+-------+-------+
+    | GCM_C | 1,2,3 |       |
+    +-------+-------+-------+
+    """
+    if dimensions is None:
+        dimensions = {}
+
+    da = da.rename(reverse_dict(dimensions))
+
+    ok = da.notnull().any("time").any("member").all("scenario")
+
+    return da.sel(model=ok).rename(dimensions)
+
+
+def single_member(da, dimensions=None):
+    """Return data for a single member per model.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Input data with dimensions for time, member, model and scenario.
+    dimensions: dict
+      Mapping from original dimension names to standard dimension names: scenario, model, member.
+
+    Returns
+    -------
+    xr.DataArray
+      Data with only one member per model.
+
+    Notes
+    -----
+    In the following example, the original members would be filtered to return only the first member found for each
+    scenario.
+
+    +-------+-------+-------+----+-------+-------+
+    | Model | Members       |    | Selected      |
+    +-------+---------------+----+---------------+
+    |       | scen1 | scen2 |    | scen1 | scen2 |
+    +=======+=======+=======+====+=======+=======+
+    | GCM_A | 1,2,3 | 1,2,3 |    | 1     | 1     |
+    +-------+-------+-------+----+-------+-------+
+    | GCM_B | 1,2   | 2,3   |    | 1     | 2     |
+    +-------+-------+-------+----+-------+-------+
+    """
+    if dimensions is None:
+        dimensions = {}
+
+    da = da.rename(reverse_dict(dimensions))
+
+    # Stack by simulation specifications - drop simulations with missing values
+    full = da.stack(i=("scenario", "model", "member")).dropna("i", how="any")
+
+    # Pick first run with data
+    s = full.i.to_series()
+    s[:] = np.arange(len(s))
+    i = s.unstack().T.min().to_list()
+
+    out = full.isel(i=i).unstack().squeeze()
+    return out.rename(dimensions)
+
+
+def reverse_dict(d):
+    """Reverse dictionary."""
+    return {v: k for (k, v) in d.items()}
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
new file mode 100644
index 000000000..cc88060ea
--- /dev/null
+++ b/xclim/ensembles/_partitioning.py
@@ -0,0 +1,186 @@
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from ._filters import reverse_dict
+
+"""
+Uncertainty Partitioning
+========================
+
+This module implements methods and tools meant to partition climate projection uncertainties into different components:
+natural variability, GHG scenario and climate models.
+
+Partitioning algorithms:
+
+ - `hawkins_sutton`
+"""
+
+"""
+# References for other more recent algorithms
+
+Yip, S., Ferro, C. A. T., Stephenson, D. B., and Hawkins, E. (2011). A Simple, Coherent Framework for Partitioning
+Uncertainty in Climate Predictions. Journal of Climate 24, 17, 4634-4643, doi:10.1175/2011JCLI4085.1
+
+Northrop, P. J., & Chandler, R. E. (2014). Quantifying sources of uncertainty in projections of future climate.
+Journal of Climate, 27(23), 8793–8808, doi:10.1175/JCLI-D-14-00265.1
+
+Goldenson, N., Mauger, G., Leung, L. R., Bitz, C. M., & Rhines, A. (2018). Effects of ensemble configuration on
+estimates of regional climate uncertainties. Geophysical Research Letters, 45, 926– 934.
+https://doi.org/10.1002/2017GL076297
+
+Lehner, F., Deser, C., Maher, N., Marotzke, J., Fischer, E. M., Brunner, L., Knutti, R., and Hawkins,
+E. (2020). Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6, Earth Syst. Dynam.,
+11, 491–508, https://doi.org/10.5194/esd-11-491-2020.
+
+Evin, G., Hingray, B., Blanchet, J., Eckert, N., Morin, S., & Verfaillie, D. (2019). Partitioning Uncertainty
+Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation, Journal of Climate, 32(8),
+2423-2440, https://doi.org/10.1175/JCLI-D-18-0606.1
+"""
+
+# Default dimension names
+_dims = dict(model="model", scenario="scenario", member="member")
+
+
+def hawkins_sutton(
+    da: xr.DataArray,
+    sm: xr.DataArray = None,
+    weights: xr.DataArray = None,
+    baseline: tuple = ("1971", "2000"),
+    kind: str = "+",
+    dimensions: dict = None,
+):
+    """Return the mean and partitioned variance of an ensemble.
+
+    Algorithm based on Hawkins and Sutton (2009). Input data should meet the following requirements:
+      - annual frequency;
+      - covers the baseline and future period;
+      - defined over time, scenario and model dimensions;
+      - the same models should be available for each scenario.
+
+    Parameters
+    ----------
+    da: xr.DataArray
+      Time series over dimensions of 'time', 'scenario' and 'model'. If other dimension names are used, provide
+      the mapping with the `dimensions` argument.
+    sm: xr.DataArray
+      Smoothed time series over time. By default, the method uses a fourth order polynomial. Use this argument to
+      use other smoothing strategies, such as another polynomial order, or a LOESS curve.
+    weights: xr.DataArray
+      Weights to be applied to individual models.
+    baseline: [str, str]
+      Start and end year of the reference period.
+    kind: {'+', '*'}
+      Whether the mean over the reference period should be substracted (+) or divided by (*).
+    dimensions: dict
+      Mapping from original dimension names to standard dimension names: scenario, model, member.
+
+    Returns
+    -------
+    xr.DataArray, xr.DataArray
+      The mean and the components of variance of the ensemble. These components are coordinates along the
+      `uncertainty` dimension: `variability`, `model`, `scenario`, and `total`.
+
+    References
+    ----------
+    Hawkins, E., & Sutton, R. (2009) The Potential to Narrow Uncertainty in Regional Climate Predictions. Bulletin of
+      the American Meteorological Society, 90(8), 1095–1108. doi:10.1175/2009BAMS2607.1
+    Hawkins, E., and R. Sutton (2011) The potential to narrow uncertainty in projections of regional precipitation
+      change. Climate Dyn., 37, 407–418, doi:10.1007/s00382-010-0810-6.
+    """
+    if xr.infer_freq(da.time)[0] not in ["A", "Y"]:
+        raise ValueError("This algorithm expects annual time series.")
+
+    if dimensions is None:
+        dimensions = {}
+
+    da = da.rename(reverse_dict(dimensions))
+
+    if "member" in da.dims and len(da.member) > 1:
+        raise ValueError("There should only be one member per model.")
+
+    # Confirm the same models have data for all scenarios
+    check = da.notnull().any("time").all("scenario")
+    if not check.all():
+        raise ValueError(f"Some models are missing data for some scenarios: \n {check}")
+
+    if weights is None:
+        weights = xr.ones_like(da.model, float)
+
+    # Perform analysis 1950 onward
+    da = da.sel(time=slice("1950", None))
+
+    if sm is None:
+        # Fit 4th order polynomial to smooth natural fluctuations
+        # Note that the order of the polynomial has a substantial influence on the results.
+        fit = da.polyfit(dim="time", deg=4, skipna=True)
+        sm = xr.polyval(coord=da.time, coeffs=fit.polyfit_coefficients).where(
+            da.notnull()
+        )
+
+    # Decadal mean residuals
+    res = (da - sm).rolling(time=10).mean(center=True)
+
+    # Individual model variance after 2000: V
+    # Note that the historical data is the same for all scenarios.
+    nv_u = (
+        res.sel(time=slice("2000", None))
+        .var(dim=("scenario", "time"))
+        .weighted(weights)
+        .mean("model")
+    )
+
+    # Compute baseline average
+    ref = sm.sel(time=slice(*baseline)).mean(dim="time")
+
+    # Remove baseline average from smoothed time series
+    if kind == "+":
+        sm -= ref
+    elif kind == "*":
+        sm /= ref
+    else:
+        raise ValueError(kind)
+
+    # Model uncertainty: M(t)
+    model_u = sm.weighted(weights).var(dim="model").mean(dim="scenario")
+
+    # Scenario uncertainty: S(t)
+    scenario_u = sm.weighted(weights).mean(dim="model").var(dim="scenario")
+
+    # Total uncertainty: T(t)
+    total = nv_u + scenario_u + model_u
+
+    # Create output array with the uncertainty components
+    u = pd.Index(["variability", "model", "scenario", "total"], name="uncertainty")
+    uncertainty = xr.concat([nv_u, model_u, scenario_u, total], dim=u)
+
+    # Mean projection: G(t)
+    g = sm.weighted(weights).mean(dim="model").mean(dim="scenario")
+
+    return g, uncertainty
+
+
+def hawkins_sutton_09_weighting(da, obs, baseline=("1971", "2000")):
+    """Return weights according to the ability of models to simulate observed climate change.
+
+    Weights are computed by comparing the 2000 value to the baseline mean: w_m = 1 / (x_{obs} + | x_{m,
+    2000} - x_obs | )
+
+    Parameters
+    ----------
+    da: xr.DataArray
+      Input data over the historical period. Should have a time and model dimension.
+    obs: float
+      Observed change.
+    baseline: (str, str)
+      Baseline start and end year.
+
+    Returns
+    -------
+    xr.DataArray
+      Weights over the model dimension.
+    """
+    mm = da.sel(time=slice(*baseline)).mean("time")
+    xm = da.sel(time=baseline[1]) - mm
+    xm = xm.drop("time").squeeze()
+    return 1 / (obs + np.abs(xm - obs))
diff --git a/xclim/testing/tests/test_partitioning.py b/xclim/testing/tests/test_partitioning.py
new file mode 100644
index 000000000..f219ccb40
--- /dev/null
+++ b/xclim/testing/tests/test_partitioning.py
@@ -0,0 +1,18 @@
+from xclim.ensembles import (
+    concat_hist,
+    hawkins_sutton,
+    model_in_all_scens,
+    single_member,
+)
+
+
+def test_hawkins_sutton(open_dataset):
+    """Just a smoke test - looking for data for a hard validation."""
+    dims = {"member": "run", "scenario": "scen"}
+    da = open_dataset(
+        "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
+    ).pr
+    da1 = model_in_all_scens(da, dimensions=dims)
+    dac = concat_hist(da1, scen="historical")
+    das = single_member(dac, dimensions=dims)
+    hawkins_sutton(das, dimensions=dims)

From 486eed4ad8d11d87d362f38d96e1da52a4feb6eb Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Tue, 20 Dec 2022 16:23:10 -0500
Subject: [PATCH 02/26] update HISTORY

---
 HISTORY.rst | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/HISTORY.rst b/HISTORY.rst
index bbf6ea63f..7f09d1200 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,15 @@
 History
 =======
 
+0.41.0 (unreleased)
+-------------------
+Contributors to this version: David Huard (:user:`huard`),
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* `ensembles.hawkins_sutton` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`)
+
+
 0.40.0 (unreleased)
 -------------------
 Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), David Huard (:user:`huard`), Juliette Lavoie (:user:`juliettelavoie`).

From 8f90f6c77663bf90af521a8a475a82ca08625db1 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Tue, 31 Jan 2023 15:02:39 +0000
Subject: [PATCH 03/26] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 xclim/ensembles/_filters.py              | 2 ++
 xclim/ensembles/_partitioning.py         | 2 ++
 xclim/testing/tests/test_partitioning.py | 2 ++
 3 files changed, 6 insertions(+)

diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
index ff88f0cf2..3cc30b46b 100644
--- a/xclim/ensembles/_filters.py
+++ b/xclim/ensembles/_filters.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 import xarray as xr
 
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index cc88060ea..745cdcb29 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import numpy as np
 import pandas as pd
 import xarray as xr
diff --git a/xclim/testing/tests/test_partitioning.py b/xclim/testing/tests/test_partitioning.py
index f219ccb40..ed1b3e80e 100644
--- a/xclim/testing/tests/test_partitioning.py
+++ b/xclim/testing/tests/test_partitioning.py
@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 from xclim.ensembles import (
     concat_hist,
     hawkins_sutton,

From c99a632047b2f340a2a634ac85285f5d5ae9855c Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Fri, 3 Feb 2023 13:13:20 -0500
Subject: [PATCH 04/26] move references to bibtex. Add module and function to
 api docs

---
 docs/api.rst                     |   9 ++-
 docs/references.bib              | 108 +++++++++++++++++++++++++++++++
 xclim/ensembles/_partitioning.py |  38 ++++++-----
 3 files changed, 138 insertions(+), 17 deletions(-)

diff --git a/docs/api.rst b/docs/api.rst
index 8487c5590..44f42157a 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -39,13 +39,20 @@ Ensembles module
     :noindex:
 
 .. automodule:: xclim.ensembles._robustness
-    :noindeX:
+    :noindex:
 
 .. autofunction:: xclim.ensembles.change_significance
     :noindex:
 .. autofunction:: xclim.ensembles.robustness_coefficient
     :noindex:
 
+.. automodule:: xclim.ensembles._partitioning
+    :noindex:
+
+.. autofunction:: xclim.ensembles.hawkins_sutton
+    :noindex:
+
+
 Indicator Tools
 ===============
 
diff --git a/docs/references.bib b/docs/references.bib
index 70780a565..5d3c6e059 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -1868,3 +1868,111 @@ @article{vrac_multivariate_2018
 	author = {Vrac, M.},
 	year = {2018},
 }
+
+@article{yip_2011,
+abstract = {A simple and coherent framework for partitioning uncertainty in multimodel climate ensembles is presented. The analysis of variance (ANOVA) is used to decompose a measure of total variation additively into scenario uncertainty, model uncertainty, and internal variability. This approach requires fewer assumptions than existing methods and can be easily used to quantify uncertainty related to model–scenario interaction—the contribution to model uncertainty arising from the variation across scenarios of model deviations from the ensemble mean. Uncertainty in global mean surface air temperature is quantified as a function of lead time for a subset of the Coupled Model Intercomparison Project phase 3 ensemble and results largely agree with those published by other authors: scenario uncertainty dominates beyond 2050 and internal variability remains approximately constant over the twenty-first century. Both elements of model uncertainty, due to scenario-independent and scenario-dependent deviations from the ensemble mean, are found to increase with time. Estimates of model deviations that arise as by-products of the framework reveal significant differences between models that could lead to a deeper understanding of the sources of uncertainty in multimodel ensembles. For example, three models show a diverging pattern over the twenty-first century, while another model exhibits an unusually large variation among its scenario-dependent deviations.},
+author = {Yip, Stan and Ferro, Christopher A. T. and Stephenson, David B. and Hawkins, Ed},
+doi = {10.1175/2011JCLI4085.1},
+issn = {0894-8755},
+journal = {Journal of Climate},
+mendeley-groups = {xclim},
+month = {sep},
+number = {17},
+pages = {4634--4643},
+title = {{A Simple, Coherent Framework for Partitioning Uncertainty in Climate Predictions}},
+url = {https://journals.ametsoc.org/doi/10.1175/2011JCLI4085.1},
+volume = {24},
+year = {2011}
+}
+
+@article{northrop_2014,
+abstract = {A simple statistical model is used to partition uncertainty from different sources, in projections of future climate from multimodel ensembles. Three major sources of uncertainty are considered: the choice of climate model, the choice of emissions scenario, and the internal variability of the modeled climate system. The relative contributions of these sources are quantified for mid- and late-twenty-first-century climate projections, using data from 23 coupled atmosphere–ocean general circulation models obtained from phase 3 of the Coupled Model Intercomparison Project (CMIP3). Similar investigations have been carried out recently by other authors but within a statistical framework for which the unbalanced nature of the data and the small number (three) of scenarios involved are potentially problematic. Here, a Bayesian analysis is used to overcome these difficulties. Global and regional analyses of surface air temperature and precipitation are performed. It is found that the relative contributions to uncertainty depend on the climate variable considered, as well as the region and time horizon. As expected, the uncertainty due to the choice of emissions scenario becomes more important toward the end of the twenty-first century. However, for midcentury temperature, model internal variability makes a large contribution in high-latitude regions. For midcentury precipitation, model internal variability is even more important and this persists in some regions into the late century. Implications for the design of climate model experiments are discussed.},
+author = {Northrop, Paul J. and Chandler, Richard E.},
+doi = {10.1175/JCLI-D-14-00265.1},
+issn = {0894-8755},
+journal = {Journal of Climate},
+mendeley-groups = {xclim},
+month = {dec},
+number = {23},
+pages = {8793--8808},
+title = {{Quantifying Sources of Uncertainty in Projections of Future Climate*}},
+url = {http://journals.ametsoc.org/doi/10.1175/JCLI-D-14-00265.1},
+volume = {27},
+year = {2014}
+}
+
+@article{goldenson_2018,
+author = {Goldenson, N. and Mauger, G. and Leung, L. R. and Bitz, C. M. and Rhines, A.},
+doi = {10.1002/2017GL076297},
+issn = {0094-8276},
+journal = {Geophysical Research Letters},
+mendeley-groups = {xclim},
+month = {jan},
+number = {2},
+pages = {926--934},
+title = {{Effects of Ensemble Configuration on Estimates of Regional Climate Uncertainties}},
+url = {https://onlinelibrary.wiley.com/doi/10.1002/2017GL076297},
+volume = {45},
+year = {2018}
+}
+
+@article{lehner_2020,
+author = {Lehner, Flavio and Deser, Clara and Maher, Nicola and Marotzke, Jochem and Fischer, Erich M. and Brunner, Lukas and Knutti, Reto and Hawkins, Ed},
+doi = {10.5194/esd-11-491-2020},
+issn = {2190-4987},
+journal = {Earth System Dynamics},
+mendeley-groups = {xclim},
+month = {may},
+number = {2},
+pages = {491--508},
+title = {{Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6}},
+url = {https://esd.copernicus.org/articles/11/491/2020/},
+volume = {11},
+year = {2020}
+}
+
+@article{evin_2019,
+abstract = {The quantification of uncertainty sources in ensembles of climate projections obtained from combinations of different scenarios and climate and impact models is a key issue in climate impact studies. The small size of the ensembles of simulation chains and their incomplete sampling of scenario and climate model combinations makes the analysis difficult. In the popular single-time ANOVA approach for instance, a precise estimate of internal variability requires multiple members for each simulation chain (e.g., each emission scenario–climate model combination), but multiple members are typically available for a few chains only. In most ensembles also, a precise partition of model uncertainty components is not possible because the matrix of available scenario/models combinations is incomplete (i.e., projections are missing for many scenario–model combinations). The method we present here, based on data augmentation and Bayesian techniques, overcomes such limitations and makes the statistical analysis possible for single-member and incomplete ensembles. It provides unbiased estimates of climate change responses of all simulation chains and of all uncertainty variables. It additionally propagates uncertainty due to missing information in the estimates. This approach is illustrated for projections of regional precipitation and temperature for four mountain massifs in France. It is applicable for any kind of ensemble of climate projections, including those produced from ad hoc impact models.},
+author = {Evin, Guillaume and Hingray, Benoit and Blanchet, Juliette and Eckert, Nicolas and Morin, Samuel and Verfaillie, Deborah},
+doi = {10.1175/JCLI-D-18-0606.1},
+issn = {0894-8755},
+journal = {Journal of Climate},
+mendeley-groups = {xclim},
+month = {apr},
+number = {8},
+pages = {2423--2440},
+title = {{Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation}},
+url = {http://journals.ametsoc.org/doi/10.1175/JCLI-D-18-0606.1},
+volume = {32},
+year = {2019}
+}
+
+@article{hawkins_2009,
+author = {Hawkins, Ed and Sutton, Rowan},
+doi = {10.1175/2009BAMS2607.1},
+issn = {0003-0007},
+journal = {Bulletin of the American Meteorological Society},
+mendeley-groups = {xclim},
+month = {aug},
+number = {8},
+pages = {1095--1108},
+title = {{The Potential to Narrow Uncertainty in Regional Climate Predictions}},
+url = {https://journals.ametsoc.org/doi/10.1175/2009BAMS2607.1},
+volume = {90},
+year = {2009}
+}
+
+@article{hawkins_2011,
+author = {Hawkins, Ed and Sutton, Rowan},
+doi = {10.1007/s00382-010-0810-6},
+issn = {0930-7575},
+journal = {Climate Dynamics},
+mendeley-groups = {xclim},
+month = {jul},
+number = {1-2},
+pages = {407--418},
+title = {{The potential to narrow uncertainty in projections of regional precipitation change}},
+url = {http://link.springer.com/10.1007/s00382-010-0810-6},
+volume = {37},
+year = {2011}
+}
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index 745cdcb29..1d658e4bc 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -1,3 +1,13 @@
+# noqa: D205,D400
+"""
+Uncertainty Partitioning
+========================
+
+This module implements methods and tools meant to partition climate projection uncertainties into different components:
+natural variability, GHG scenario and climate models.
+"""
+
+
 from __future__ import annotations
 
 import numpy as np
@@ -7,19 +17,11 @@
 from ._filters import reverse_dict
 
 """
-Uncertainty Partitioning
-========================
-
-This module implements methods and tools meant to partition climate projection uncertainties into different components:
-natural variability, GHG scenario and climate models.
-
-Partitioning algorithms:
+Implemented partitioning algorithms:
 
  - `hawkins_sutton`
-"""
 
-"""
-# References for other more recent algorithms
+# References for other more recent algorithms that could be added here.
 
 Yip, S., Ferro, C. A. T., Stephenson, D. B., and Hawkins, E. (2011). A Simple, Coherent Framework for Partitioning
 Uncertainty in Climate Predictions. Journal of Climate 24, 17, 4634-4643, doi:10.1175/2011JCLI4085.1
@@ -38,6 +40,13 @@
 Evin, G., Hingray, B., Blanchet, J., Eckert, N., Morin, S., & Verfaillie, D. (2019). Partitioning Uncertainty
 Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation, Journal of Climate, 32(8),
 2423-2440, https://doi.org/10.1175/JCLI-D-18-0606.1
+
+Related bixtex entries:
+ - yip_2011
+ - northrop_2014
+ - goldenson_2018
+ - lehner_2020
+ - evin_2019
 """
 
 # Default dimension names
@@ -52,9 +61,9 @@ def hawkins_sutton(
     kind: str = "+",
     dimensions: dict = None,
 ):
-    """Return the mean and partitioned variance of an ensemble.
+    """Return the mean and partitioned variance of an ensemble based on method from Hawkins & Sutton.
 
-    Algorithm based on Hawkins and Sutton (2009). Input data should meet the following requirements:
+    Algorithm based on :cite:t:`hawkins_2009`. Input data should meet the following requirements:
       - annual frequency;
       - covers the baseline and future period;
       - defined over time, scenario and model dimensions;
@@ -85,10 +94,7 @@ def hawkins_sutton(
 
     References
     ----------
-    Hawkins, E., & Sutton, R. (2009) The Potential to Narrow Uncertainty in Regional Climate Predictions. Bulletin of
-      the American Meteorological Society, 90(8), 1095–1108. doi:10.1175/2009BAMS2607.1
-    Hawkins, E., and R. Sutton (2011) The potential to narrow uncertainty in projections of regional precipitation
-      change. Climate Dyn., 37, 407–418, doi:10.1007/s00382-010-0810-6.
+    :cite:cts:`hawkins_2009,hawkins_2011`
     """
     if xr.infer_freq(da.time)[0] not in ["A", "Y"]:
         raise ValueError("This algorithm expects annual time series.")

From 56ad771d2e729079ecc1934f8f8670668311b75c Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Tue, 28 Feb 2023 09:30:01 -0500
Subject: [PATCH 05/26] Apply suggestions from code review

Co-authored-by: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com>
Co-authored-by: Pascal Bourgault <bourgault.pascal@ouranos.ca>
---
 HISTORY.rst                      | 2 +-
 xclim/ensembles/_partitioning.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/HISTORY.rst b/HISTORY.rst
index dcdfa5e07..1b6a9b461 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -4,7 +4,7 @@ History
 
 0.41.0 (unreleased)
 -------------------
-
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Ludwig Lierhammer (:user:`ludwiglierhammer`), Éric Dupuis (:user:`coxipi`).
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index 1d658e4bc..443556eb5 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -127,7 +127,7 @@ def hawkins_sutton(
         )
 
     # Decadal mean residuals
-    res = (da - sm).rolling(time=10).mean(center=True)
+    res = (da - sm).rolling(time=10, center=True).mean()
 
     # Individual model variance after 2000: V
     # Note that the historical data is the same for all scenarios.

From e6ac6f292cb6fddaadc27f50d3279402cd29f681 Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Wed, 8 Mar 2023 15:36:51 -0500
Subject: [PATCH 06/26] update CHANGES.rst

---
 CHANGES.rst | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index b4bdf2f45..25c152198 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -4,7 +4,11 @@ Changelog
 
 0.42.0 (unreleased)
 -------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), David Huard (:user:`huard`).
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
@@ -24,9 +28,6 @@ Internal changes
 0.41.0 (2023-02-28)
 -------------------
 Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Ludwig Lierhammer (:user:`ludwiglierhammer`), Éric Dupuis (:user:`coxipi`).
-New features and enhancements
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-* ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -36,7 +37,7 @@ New indicators
 ^^^^^^^^^^^^^^
 * New indices and indicators for converting from snow water equivalent to snow depth (``snw_to_snd``) and snow depth to snow water equivalent (``snd_to_snw``) using snow density [kg/m^3]. (:pull:`1271`).
 * New indices and indicators for determining upwelling radiation (`shortwave_upwelling_radiation_from_net_downwelling` and `longwave_upwelling_radiation_from_net_downwelling`; CF variables `rsus` and `rlus`) from net and downwelling radiation (shortwave: `rss` and `rsds`; longwave: `rls` and `rlds`). (:pull:`1271`).
-* New indice and indicator ``{snd | snw}_season_{length | start | end}`` which generalize ``snow_cover_duration`` and  ``continuous_snow_cover_{start | end}`` to allow using these functions with variable `snw` (:pull:`1275`).
+* New indice and indicator ``{snd | snw}_season_{length | start | end}`` which generalize ``snow_cover_duration`` and ``continuous_snow_cover_{start | end}`` to allow using these functions with variable `snw` (:pull:`1275`).
 * New indice and indicator (``dryness_index``) for estimating soil humidity classifications for winegrowing regions (based on Riou et al. (1994)). (:issue:`355`, :pull:`1235`).
 * ``ensembles.change_significance`` now supports Mann-whitney U-test and flexible ``realization``. (:pull:`1285`).
 
@@ -54,7 +55,7 @@ Breaking changes
 
 Bug fixes
 ^^^^^^^^^
-* ``build_indicator_module_from_yaml`` now accepts a ``reload`` argument. When re-building a module that already exists,  ``reload=True`` removes all previous indicator before creating the new ones. (:issue:`1192`,:pull:`1284`).
+* ``build_indicator_module_from_yaml`` now accepts a ``reload`` argument. When re-building a module that already exists, ``reload=True`` removes all previous indicator before creating the new ones. (:issue:`1192`,:pull:`1284`).
 * The test for french translations of official indicators was fixed and translations for CFFWIS indices, FFDI, KDBI, DF and Jetstream metric woollings have been added or fixed. (:pull:`1271`).
 * ``use_ufunc`` in ``windowed_run_count`` is now supplied with argument ``freq`` to warn users that the 1d method does not support resampling after run length operations (:issue:`1279`, :pull:`1291`).
 * ``{snd|snw}_max_doy`` now avoids an error due to `xr.argmax` when there are all-NaN slices. (:pull:`1277`).
@@ -488,7 +489,7 @@ Internal changes
 ^^^^^^^^^^^^^^^^
 * Added a CI hook in ``.pre-commit-config.yaml`` to perform automated `pre-commit` corrections with GitHub CI. (:pull:`965`).
 * Adjusted CI hooks to fail earlier if `lint` checks fail. (:pull:`972`).
-* `TrainAdjust` and `Adjust` object have a new `skip_input_checks` keyword arg to their `train` and  `adjust` methods. When `True`, all unit-, calendar- and coordinate-related input checks are skipped. This is an ugly solution to disappearing attributes when using `xr.map_blocks` with dask. (:pull:`964`).
+* `TrainAdjust` and `Adjust` object have a new `skip_input_checks` keyword arg to their `train` and `adjust` methods. When `True`, all unit-, calendar- and coordinate-related input checks are skipped. This is an ugly solution to disappearing attributes when using `xr.map_blocks` with dask. (:pull:`964`).
 * Some slow tests were marked `slow` to help speed up the standard test ensemble. (:pull:`969`).
     - Tox testing ensemble now also reports slowest tests using the ``--durations`` flag.
 * `pint` no longer emits warnings about redefined units when the `logging` module is loaded. (:issue:`990`, :pull:`991`).
@@ -645,7 +646,7 @@ New features and enhancements
 * ``xclim.core.utils._cal_perc`` is now only a proxy for ``xc.core.utils.nan_calc_percentiles`` with some axis moves.
 * `xclim` now implements many data quality assurance flags (``xclim.core.dataflags``) for temperature and precipitation based on `ICCLIM documentation guidelines <https://www.ecad.eu/documents/atbd.pdf>`_. These checks include the following:
     - Temperature (variables: ``tas``, ``tasmin``, ``tasmax``): ``tasmax_below_tasmin``, ``tas_exceeds_tasmax``, ``tas_below_tasmin``, ``temperature_extremely_low`` (`thresh="-90 degC"`), ``temperature_extremely_high`` (`thresh="60 degC"`).
-    - Precipitation-specific (variables: ``pr``, ``prsn``, ):  ``negative_accumulation_values``, ``very_large_precipitation_events`` (`thresh="300 mm d-1"`).
+    - Precipitation-specific (variables: ``pr``, ``prsn``, ): ``negative_accumulation_values``, ``very_large_precipitation_events`` (`thresh="300 mm d-1"`).
     - Wind-specific (variables: ``sfcWind``, ``wsgsmax``/``sfcWindMax``): ``wind_values_outside_of_bounds``
     - Generic: ``outside_n_standard_deviations_of_climatology``, ``values_repeating_for_n_or_more_days``, ``values_op_thresh_repeating_for_n_or_more_days``, ``percentage_values_outside_of_bounds``.
 
@@ -762,8 +763,8 @@ New features and enhancements
 * ``humidex`` can be computed using relative humidity instead of dewpoint temperature.
 * New ``sdba.construct_moving_yearly_window`` and ``sdba.unpack_moving_yearly_window`` for moving window adjustments.
 * New ``sdba.adjustment.NpdfTransform`` which is an adaptation of Alex Cannon's version of Pitié's *N-dimensional probability density function transform*. Uses new ``sdba.utils.rand_rot_matrix``. *Experimental, subject to changes.*
-* New ``sdba.processing.standardize``, ``.unstandardize`` and  ``.reordering``. All of them, tools needed to replicate Cannon's MBCn algorithm.
-* New ``sdba.processing.escore``, backed by  ``sdba.nbutils._escore`` to evaluate the performance of the N pdf transform.
+* New ``sdba.processing.standardize``, ``.unstandardize`` and ``.reordering``. All of them, tools needed to replicate Cannon's MBCn algorithm.
+* New ``sdba.processing.escore``, backed by ``sdba.nbutils._escore`` to evaluate the performance of the N pdf transform.
 * New function ``xclim.indices.clausius_clapeyron_scaled_precipitation`` can be used to scale precipitation according to changes in mean temperature.
 * Percentile based indices gained a ``bootstrap`` argument that applies a bootstrapping algorithm to reduce biases on exceedance frequencies computed over *in base* and *out of base* periods. *Experimental, subject to changes.*
 * Added a `.zenodo.json` file for collecting and maintaining author order and tracking ORCIDs.
@@ -1215,7 +1216,7 @@ v0.15.x (2020-03-12)
 --------------------
 * Improvement in FWI: Vectorization of DC, DMC and FFMC with numba and small code refactoring for better maintainability.
 * Added example notebook for creating a catalog of selected indices
-* Added `growing_season_end`, `last_spring_frost`, `dry_days`,  `hot_spell_frequency`, `hot_spell_max_length`, and `maximum_consecutive_frost_free_days` indices.
+* Added `growing_season_end`, `last_spring_frost`, `dry_days`, `hot_spell_frequency`, `hot_spell_max_length`, and `maximum_consecutive_frost_free_days` indices.
 * Dropped use of `fiona.crs` class in lieu of the newer pyproj CRS handler for `subset_shape` operations.
 * Complete internal reorganization of xclim.
 * Internationalization of xclim : add `locales` submodule for localized metadata.

From d4c27c2725665d88319f504a06aa1df3dbacf8b9 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Mon, 3 Apr 2023 17:01:24 -0400
Subject: [PATCH 07/26] changes based on review.

---
 xclim/ensembles/_partitioning.py         | 40 +++++++++---------------
 xclim/testing/tests/test_partitioning.py | 20 +++++++-----
 2 files changed, 27 insertions(+), 33 deletions(-)

diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index 443556eb5..6cc26d12b 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -14,8 +14,6 @@
 import pandas as pd
 import xarray as xr
 
-from ._filters import reverse_dict
-
 """
 Implemented partitioning algorithms:
 
@@ -49,9 +47,6 @@
  - evin_2019
 """
 
-# Default dimension names
-_dims = dict(model="model", scenario="scenario", member="member")
-
 
 def hawkins_sutton(
     da: xr.DataArray,
@@ -59,21 +54,18 @@ def hawkins_sutton(
     weights: xr.DataArray = None,
     baseline: tuple = ("1971", "2000"),
     kind: str = "+",
-    dimensions: dict = None,
 ):
-    """Return the mean and partitioned variance of an ensemble based on method from Hawkins & Sutton.
+    """Return the mean and partitioned variance of an ensemble based on method from Hawkins & Sutton (2009).
 
-    Algorithm based on :cite:t:`hawkins_2009`. Input data should meet the following requirements:
-      - annual frequency;
+    To reproduce results from :cite:t:`hawkins_2009`, input data should meet the following requirements:
+      - annual frequency, starting in 1950;
       - covers the baseline and future period;
-      - defined over time, scenario and model dimensions;
-      - the same models should be available for each scenario.
+      - same models available for each scenario.
 
     Parameters
     ----------
     da: xr.DataArray
-      Time series over dimensions of 'time', 'scenario' and 'model'. If other dimension names are used, provide
-      the mapping with the `dimensions` argument.
+      Time series over dimensions 'time', 'scenario' and 'model'.
     sm: xr.DataArray
       Smoothed time series over time. By default, the method uses a fourth order polynomial. Use this argument to
       use other smoothing strategies, such as another polynomial order, or a LOESS curve.
@@ -83,8 +75,6 @@ def hawkins_sutton(
       Start and end year of the reference period.
     kind: {'+', '*'}
       Whether the mean over the reference period should be substracted (+) or divided by (*).
-    dimensions: dict
-      Mapping from original dimension names to standard dimension names: scenario, model, member.
 
     Returns
     -------
@@ -92,6 +82,12 @@ def hawkins_sutton(
       The mean and the components of variance of the ensemble. These components are coordinates along the
       `uncertainty` dimension: `variability`, `model`, `scenario`, and `total`.
 
+    Notes
+    -----
+    To prepare input data, make sure `da` has dimensions `time`, `scenario` and `model`,
+    e.g. `da.rename({"scen": "scenario"})`. Also, to match the paper's methodology, select years from 1950 onward,
+    e.g. `da.sel(time=slice("1950", None))`.
+
     References
     ----------
     :cite:cts:`hawkins_2009,hawkins_2011`
@@ -99,13 +95,10 @@ def hawkins_sutton(
     if xr.infer_freq(da.time)[0] not in ["A", "Y"]:
         raise ValueError("This algorithm expects annual time series.")
 
-    if dimensions is None:
-        dimensions = {}
-
-    da = da.rename(reverse_dict(dimensions))
-
-    if "member" in da.dims and len(da.member) > 1:
-        raise ValueError("There should only be one member per model.")
+    if not {"time", "scenario", "model"}.issubset(da.dims):
+        raise ValueError(
+            "DataArray dimensions should include 'time', 'scenario' and 'model'."
+        )
 
     # Confirm the same models have data for all scenarios
     check = da.notnull().any("time").all("scenario")
@@ -115,9 +108,6 @@ def hawkins_sutton(
     if weights is None:
         weights = xr.ones_like(da.model, float)
 
-    # Perform analysis 1950 onward
-    da = da.sel(time=slice("1950", None))
-
     if sm is None:
         # Fit 4th order polynomial to smooth natural fluctuations
         # Note that the order of the polynomial has a substantial influence on the results.
diff --git a/xclim/testing/tests/test_partitioning.py b/xclim/testing/tests/test_partitioning.py
index ed1b3e80e..49e0682ea 100644
--- a/xclim/testing/tests/test_partitioning.py
+++ b/xclim/testing/tests/test_partitioning.py
@@ -10,11 +10,15 @@
 
 def test_hawkins_sutton(open_dataset):
     """Just a smoke test - looking for data for a hard validation."""
-    dims = {"member": "run", "scenario": "scen"}
-    da = open_dataset(
-        "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
-    ).pr
-    da1 = model_in_all_scens(da, dimensions=dims)
-    dac = concat_hist(da1, scen="historical")
-    das = single_member(dac, dimensions=dims)
-    hawkins_sutton(das, dimensions=dims)
+    dims = {"run": "member", "scen": "scenario"}
+    da = (
+        open_dataset(
+            "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
+        )
+        .pr.sel(time=slice("1950", None))
+        .rename(dims)
+    )
+    da1 = model_in_all_scens(da)
+    dac = concat_hist(da1, scenario="historical")
+    das = single_member(dac)
+    hawkins_sutton(das)

From 380ac6868b6828484a55d706eb46af998b4a984a Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 24 May 2023 18:24:06 +0000
Subject: [PATCH 08/26] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 xclim/ensembles/__init__.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/xclim/ensembles/__init__.py b/xclim/ensembles/__init__.py
index d273be080..91c31b0d3 100644
--- a/xclim/ensembles/__init__.py
+++ b/xclim/ensembles/__init__.py
@@ -18,5 +18,4 @@
     make_criteria,
     plot_rsqprofile,
 )
-
 from ._robustness import change_significance, robustness_coefficient

From c342ba4d8f92d7ee0eee8b1641a1cc95a9a4deca Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 14 Jun 2023 16:42:59 -0400
Subject: [PATCH 09/26] start writing test for hawkins-sutton

---
 tests/test_partitioning.py               | 55 ++++++++++++++++++++++++
 xclim/ensembles/_partitioning.py         |  2 +-
 xclim/testing/tests/test_partitioning.py | 24 -----------
 3 files changed, 56 insertions(+), 25 deletions(-)
 create mode 100644 tests/test_partitioning.py
 delete mode 100644 xclim/testing/tests/test_partitioning.py

diff --git a/tests/test_partitioning.py b/tests/test_partitioning.py
new file mode 100644
index 000000000..59b34de70
--- /dev/null
+++ b/tests/test_partitioning.py
@@ -0,0 +1,55 @@
+from __future__ import annotations
+
+import numpy as np
+import xarray as xr
+
+from xclim.ensembles import (
+    concat_hist,
+    hawkins_sutton,
+    model_in_all_scens,
+    single_member,
+)
+
+
+def test_hawkins_sutton_smoke(open_dataset):
+    """Just a smoke test - looking for data for a hard validation."""
+    dims = {"run": "member", "scen": "scenario"}
+    da = (
+        open_dataset(
+            "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
+        )
+        .pr.sel(time=slice("1950", None))
+        .rename(dims)
+    )
+    da1 = model_in_all_scens(da)
+    dac = concat_hist(da1, scenario="historical")
+    das = single_member(dac)
+    hawkins_sutton(das)
+
+
+def test_hawkins_sutton_synthetic():
+    """Create synthetic data to test logic of Hawkins-Sutton's implementation."""
+    r = np.random.randn(4, 11, 60)
+    # Time, scenario, model
+    sm = np.arange(10, 41, 10)
+    mm = np.arange(-5, 6, 1)
+    mean = mm[np.newaxis, :] + sm[:, np.newaxis]
+
+    # Here the scenarios don't change over time, so there should be no model variability (since it's relative to the
+    # reference period.
+    x = r + mean[:, :, np.newaxis]
+    time = xr.date_range("1970-01-01", periods=60, freq="Y")
+    da = xr.DataArray(x, dims=("scenario", "model", "time"), coords={"time": time})
+    m, v = hawkins_sutton(da)
+
+    np.testing.assert_array_almost_equal(m, 0, decimal=1)
+
+    vm = v.mean(dim="time")
+
+    np.testing.assert_array_almost_equal(vm.sel(uncertainty="scenario"), 0, decimal=2)
+    np.testing.assert_array_almost_equal(
+        vm.sel(uncertainty="variability"), 0, decimal=1
+    )
+    np.testing.assert_array_almost_equal(
+        vm.sel(uncertainty="model"), vm.sel(uncertainty="variability"), decimal=1
+    )
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index 6cc26d12b..49050adc0 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -70,7 +70,7 @@ def hawkins_sutton(
       Smoothed time series over time. By default, the method uses a fourth order polynomial. Use this argument to
       use other smoothing strategies, such as another polynomial order, or a LOESS curve.
     weights: xr.DataArray
-      Weights to be applied to individual models.
+      Weights to be applied to individual models. Should have `model` dimension.
     baseline: [str, str]
       Start and end year of the reference period.
     kind: {'+', '*'}
diff --git a/xclim/testing/tests/test_partitioning.py b/xclim/testing/tests/test_partitioning.py
deleted file mode 100644
index 49e0682ea..000000000
--- a/xclim/testing/tests/test_partitioning.py
+++ /dev/null
@@ -1,24 +0,0 @@
-from __future__ import annotations
-
-from xclim.ensembles import (
-    concat_hist,
-    hawkins_sutton,
-    model_in_all_scens,
-    single_member,
-)
-
-
-def test_hawkins_sutton(open_dataset):
-    """Just a smoke test - looking for data for a hard validation."""
-    dims = {"run": "member", "scen": "scenario"}
-    da = (
-        open_dataset(
-            "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
-        )
-        .pr.sel(time=slice("1950", None))
-        .rename(dims)
-    )
-    da1 = model_in_all_scens(da)
-    dac = concat_hist(da1, scenario="historical")
-    das = single_member(dac)
-    hawkins_sutton(das)

From e716db6bacfb3cc33ff6d5e841e1411354d554ae Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Tue, 20 Jun 2023 09:57:53 -0400
Subject: [PATCH 10/26] Synthetic test

---
 CHANGES.rst                      |  2 --
 tests/test_partitioning.py       | 51 ++++++++++++++++++++++----------
 xclim/ensembles/_partitioning.py | 26 ++++++++--------
 3 files changed, 50 insertions(+), 29 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 4acf88008..3156ae8a7 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -2,7 +2,6 @@
 Changelog
 =========
 
-
 v0.45.0 (unreleased)
 --------------------
 * ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
@@ -11,7 +10,6 @@ v0.44.0 (unreleased)
 --------------------
 Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Ludwig Lierhammer (:user:`ludwiglierhammer`), David Huard (:user:`huard`).
 
-
 Announcements
 ^^^^^^^^^^^^^
 * `xclim: xarray-based climate data analytics` has been published in the Journal of Open Source Software (`DOI:10.21105/joss.05415 <https://doi.org/10.21105/joss.05415>`_). Users can now make use of the `Cite this repository` button in the sidebar for academic purposes. Many thanks to our core developers and user base for their fine contributions over the years! (:issue:`95`, :pull:`250`).
diff --git a/tests/test_partitioning.py b/tests/test_partitioning.py
index 59b34de70..134b8b23b 100644
--- a/tests/test_partitioning.py
+++ b/tests/test_partitioning.py
@@ -12,7 +12,7 @@
 
 
 def test_hawkins_sutton_smoke(open_dataset):
-    """Just a smoke test - looking for data for a hard validation."""
+    """Just a smoke test."""
     dims = {"run": "member", "scen": "scenario"}
     da = (
         open_dataset(
@@ -28,28 +28,49 @@ def test_hawkins_sutton_smoke(open_dataset):
 
 
 def test_hawkins_sutton_synthetic():
-    """Create synthetic data to test logic of Hawkins-Sutton's implementation."""
-    r = np.random.randn(4, 11, 60)
-    # Time, scenario, model
-    sm = np.arange(10, 41, 10)
-    mm = np.arange(-5, 6, 1)
-    mean = mm[np.newaxis, :] + sm[:, np.newaxis]
+    """Test logic of Hawkins-Sutton's implementation using synthetic data."""
 
+    # Time, scenario, model
     # Here the scenarios don't change over time, so there should be no model variability (since it's relative to the
     # reference period.
+    sm = np.arange(10, 41, 10)  # Scenario mean
+    mm = np.arange(-6, 7, 1)  # Model mean
+    mean = mm[np.newaxis, :] + sm[:, np.newaxis]
+
+    # Natural variability
+    r = np.random.randn(4, 13, 60)
+
     x = r + mean[:, :, np.newaxis]
     time = xr.date_range("1970-01-01", periods=60, freq="Y")
     da = xr.DataArray(x, dims=("scenario", "model", "time"), coords={"time": time})
     m, v = hawkins_sutton(da)
+    # Mean uncertainty over time
+    vm = v.mean(dim="time")
 
-    np.testing.assert_array_almost_equal(m, 0, decimal=1)
+    # Check that the mean relative to the baseline is zero
+    np.testing.assert_array_almost_equal(m.mean(dim="time"), 0, decimal=1)
 
-    vm = v.mean(dim="time")
+    # Check that the scenario uncertainty is zero
+    np.testing.assert_array_almost_equal(vm.sel(uncertainty="scenario"), 0, decimal=1)
 
-    np.testing.assert_array_almost_equal(vm.sel(uncertainty="scenario"), 0, decimal=2)
-    np.testing.assert_array_almost_equal(
-        vm.sel(uncertainty="variability"), 0, decimal=1
-    )
-    np.testing.assert_array_almost_equal(
-        vm.sel(uncertainty="model"), vm.sel(uncertainty="variability"), decimal=1
+    # Check that model uncertainty > variability
+    assert vm.sel(uncertainty="model") > vm.sel(uncertainty="variability")
+
+    # Smoke test with polynomial of order 2
+    fit = da.polyfit(dim="time", deg=2, skipna=True)
+    sm = xr.polyval(coord=da.time, coeffs=fit.polyfit_coefficients).where(da.notnull())
+    hawkins_sutton(da, sm=sm)
+
+    # Test with a multiplicative variable and time evolving scenarios
+    r = np.random.randn(4, 13, 60) + np.arange(60)
+    x = r + mean[:, :, np.newaxis]
+    da = xr.DataArray(x, dims=("scenario", "model", "time"), coords={"time": time})
+    m, v = hawkins_sutton(da, kind="*")
+    su = v.sel(uncertainty="scenario")
+    # We expect the scenario uncertainty to grow over time
+    # The scenarios all have the same absolute slope, but since their reference mean is different, the relative increase
+    # is not the same and this creates a spread over time across "relative" scenarios.
+    assert (
+        su.sel(time=slice("2020", None)).mean()
+        > su.sel(time=slice("2000", "2010")).mean()
     )
diff --git a/xclim/ensembles/_partitioning.py b/xclim/ensembles/_partitioning.py
index 49050adc0..72f4475ce 100644
--- a/xclim/ensembles/_partitioning.py
+++ b/xclim/ensembles/_partitioning.py
@@ -39,6 +39,9 @@
 Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation, Journal of Climate, 32(8),
 2423-2440, https://doi.org/10.1175/JCLI-D-18-0606.1
 
+Beigi E, Tsai FT-C, Singh VP, Kao S-C. Bayesian Hierarchical Model Uncertainty Quantification for Future Hydroclimate
+Projections in Southern Hills-Gulf Region, USA. Water. 2019; 11(2):268. https://doi.org/10.3390/w11020268
+
 Related bixtex entries:
  - yip_2011
  - northrop_2014
@@ -57,18 +60,14 @@ def hawkins_sutton(
 ):
     """Return the mean and partitioned variance of an ensemble based on method from Hawkins & Sutton (2009).
 
-    To reproduce results from :cite:t:`hawkins_2009`, input data should meet the following requirements:
-      - annual frequency, starting in 1950;
-      - covers the baseline and future period;
-      - same models available for each scenario.
-
     Parameters
     ----------
     da: xr.DataArray
-      Time series over dimensions 'time', 'scenario' and 'model'.
+      Time series with dimensions 'time', 'scenario' and 'model'.
     sm: xr.DataArray
-      Smoothed time series over time. By default, the method uses a fourth order polynomial. Use this argument to
-      use other smoothing strategies, such as another polynomial order, or a LOESS curve.
+      Smoothed time series over time, with the same dimensions as `da`. By default, this is estimated using a 4th order
+      polynomial. Results are sensitive to the choice of smoothing function, use this to set another polynomial
+      order, or a LOESS curve.
     weights: xr.DataArray
       Weights to be applied to individual models. Should have `model` dimension.
     baseline: [str, str]
@@ -79,14 +78,17 @@ def hawkins_sutton(
     Returns
     -------
     xr.DataArray, xr.DataArray
-      The mean and the components of variance of the ensemble. These components are coordinates along the
-      `uncertainty` dimension: `variability`, `model`, `scenario`, and `total`.
+      The mean relative to the baseline, and the components of variance of the ensemble. These components are
+      coordinates along the `uncertainty` dimension: `variability`, `model`, `scenario`, and `total`.
 
     Notes
     -----
     To prepare input data, make sure `da` has dimensions `time`, `scenario` and `model`,
-    e.g. `da.rename({"scen": "scenario"})`. Also, to match the paper's methodology, select years from 1950 onward,
-    e.g. `da.sel(time=slice("1950", None))`.
+    e.g. `da.rename({"scen": "scenario"})`.
+
+    To reproduce results from :cite:t:`hawkins_2009`, input data should meet the following requirements:
+      - annual time series starting in 1950 and ending in 2100;
+      - the same models are available for all scenarios.
 
     References
     ----------

From bbc0f6ce21500ccabd026a43681c3475a1732199 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 21 Jun 2023 10:33:54 -0400
Subject: [PATCH 11/26] make filter functions private

---
 tests/test_filters.py       |  0
 tests/test_partitioning.py  | 14 +++++---------
 xclim/ensembles/__init__.py |  1 -
 xclim/ensembles/_filters.py |  8 +++++---
 4 files changed, 10 insertions(+), 13 deletions(-)
 create mode 100644 tests/test_filters.py

diff --git a/tests/test_filters.py b/tests/test_filters.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/test_partitioning.py b/tests/test_partitioning.py
index 134b8b23b..e5e9e7a58 100644
--- a/tests/test_partitioning.py
+++ b/tests/test_partitioning.py
@@ -3,12 +3,8 @@
 import numpy as np
 import xarray as xr
 
-from xclim.ensembles import (
-    concat_hist,
-    hawkins_sutton,
-    model_in_all_scens,
-    single_member,
-)
+from xclim.ensembles import hawkins_sutton
+from xclim.ensembles._filters import _concat_hist, _model_in_all_scens, _single_member
 
 
 def test_hawkins_sutton_smoke(open_dataset):
@@ -21,9 +17,9 @@ def test_hawkins_sutton_smoke(open_dataset):
         .pr.sel(time=slice("1950", None))
         .rename(dims)
     )
-    da1 = model_in_all_scens(da)
-    dac = concat_hist(da1, scenario="historical")
-    das = single_member(dac)
+    da1 = _model_in_all_scens(da)
+    dac = _concat_hist(da1, scenario="historical")
+    das = _single_member(dac)
     hawkins_sutton(das)
 
 
diff --git a/xclim/ensembles/__init__.py b/xclim/ensembles/__init__.py
index 91c31b0d3..57c3e8d0f 100644
--- a/xclim/ensembles/__init__.py
+++ b/xclim/ensembles/__init__.py
@@ -10,7 +10,6 @@
 from __future__ import annotations
 
 from ._base import create_ensemble, ensemble_mean_std_max_min, ensemble_percentiles
-from ._filters import concat_hist, model_in_all_scens, single_member
 from ._partitioning import hawkins_sutton
 from ._reduce import (
     kkz_reduce_ensemble,
diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
index 3cc30b46b..d2f9b53ba 100644
--- a/xclim/ensembles/_filters.py
+++ b/xclim/ensembles/_filters.py
@@ -4,7 +4,7 @@
 import xarray as xr
 
 
-def concat_hist(da, **hist):
+def _concat_hist(da, **hist):
     """Concatenate historical scenario with future scenarios along time.
 
     Parameters
@@ -45,8 +45,10 @@ def concat_hist(da, **hist):
     if len(hist) > 1:
         raise ValueError("Too many values in hist scenario.")
 
+    # Scenario dimension, and name of the historical scenario
     ((dim, name),) = hist.items()
 
+    # Select historical scenario and drop it from the data
     h = da.sel(**hist).dropna("time", how="all")
     ens = da.drop_sel(**hist)
 
@@ -56,7 +58,7 @@ def concat_hist(da, **hist):
     return xr.concat([h, bare], dim="time").assign_coords({dim: index})
 
 
-def model_in_all_scens(da, dimensions=None):
+def _model_in_all_scens(da, dimensions=None):
     """Return data with only simulations that have at least one member in each scenario.
 
     Parameters
@@ -97,7 +99,7 @@ def model_in_all_scens(da, dimensions=None):
     return da.sel(model=ok).rename(dimensions)
 
 
-def single_member(da, dimensions=None):
+def _single_member(da, dimensions=None):
     """Return data for a single member per model.
 
     Parameters

From 83e02c5052cf595dc7ccd227431718d9f3cf9903 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 21 Jun 2023 10:44:19 -0400
Subject: [PATCH 12/26] docstring clarifications

---
 xclim/ensembles/_filters.py | 70 ++++++++++++++++++-------------------
 1 file changed, 35 insertions(+), 35 deletions(-)

diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
index d2f9b53ba..ff47bbcc3 100644
--- a/xclim/ensembles/_filters.py
+++ b/xclim/ensembles/_filters.py
@@ -12,7 +12,7 @@ def _concat_hist(da, **hist):
     da : xr.DataArray
       Input data where the historical scenario is stored alongside other, future, scenarios.
     hist: {str: str}
-      Mapping of the scenario dimension name to the historical scenario coordinate, e.g. `scen="historical"`.
+      Mapping of the scenario dimension name to the historical scenario coordinate, e.g. `scenario="historical"`.
 
     Returns
     -------
@@ -23,23 +23,23 @@ def _concat_hist(da, **hist):
     -----
     Data goes from:
 
-        +----------+----------------------------------+
-        | Scenario | Time                             +
-        +==========+==================================+
-        | hist     | hhhhhhhhhhhhhhhh                 |
-        +----------+----------------------------------+
-        | scen1    |                 111111111111     |
-        +----------+----------------------------------+
-        | scen2    |                 222222222222     |
-        +----------+----------------------------------+
+        +------------+----------------------------------+
+        | scenario   | time                             |
+        +============+==================================+
+        | historical | hhhhhhhhhhhhhhhh                 |
+        +------------+----------------------------------+
+        | ssp245     |                 111111111111     |
+        +------------+----------------------------------+
+        | ssp370     |                 222222222222     |
+        +------------+----------------------------------+
 
     to:
         +----------+----------------------------------+
-        | Scenario | Time                             +
+        | scenario | time                             |
         +==========+==================================+
-        | scen1    | hhhhhhhhhhhhhhhh111111111111     |
+        | ssp245   | hhhhhhhhhhhhhhhh111111111111     |
         +----------+----------------------------------+
-        | scen2    | hhhhhhhhhhhhhhhh222222222222     |
+        | ssp370   | hhhhhhhhhhhhhhhh222222222222     |
         +----------+----------------------------------+
     """
     if len(hist) > 1:
@@ -75,19 +75,19 @@ def _model_in_all_scens(da, dimensions=None):
 
     Notes
     -----
-    In the following example, `GCM_C` would be filtered out from the data because it has no member for `scen2`.
-
-    +-------+-------+-------+
-    | Model | Members       |
-    +-------+---------------+
-    |       | scen1 | scen2 |
-    +=======+=======+=======+
-    | GCM_A | 1,2,3 | 1,2,3 |
-    +-------+-------+-------+
-    | GCM_B | 1     | 2,3   |
-    +-------+-------+-------+
-    | GCM_C | 1,2,3 |       |
-    +-------+-------+-------+
+    In the following example, model `C` would be filtered out from the data because it has no member for `ssp370`.
+
+    +-------+--------+--------+
+    | model | members         |
+    +-------+-----------------+
+    |       | ssp245 | ssp370 |
+    +=======+========+========+
+    | A     | 1,2,3  | 1,2,3  |
+    +-------+--------+--------+
+    | B     | 1      | 2,3    |
+    +-------+--------+--------+
+    | C     | 1,2,3  |        |
+    +-------+--------+--------+
     """
     if dimensions is None:
         dimensions = {}
@@ -119,15 +119,15 @@ def _single_member(da, dimensions=None):
     In the following example, the original members would be filtered to return only the first member found for each
     scenario.
 
-    +-------+-------+-------+----+-------+-------+
-    | Model | Members       |    | Selected      |
-    +-------+---------------+----+---------------+
-    |       | scen1 | scen2 |    | scen1 | scen2 |
-    +=======+=======+=======+====+=======+=======+
-    | GCM_A | 1,2,3 | 1,2,3 |    | 1     | 1     |
-    +-------+-------+-------+----+-------+-------+
-    | GCM_B | 1,2   | 2,3   |    | 1     | 2     |
-    +-------+-------+-------+----+-------+-------+
+    +-------+--------+--------+----+--------+--------+
+    | model | member          |    | Selected        |
+    +-------+-----------------+----+-----------------+
+    |       | ssp245 | ssp370 |    | ssp245 | ssp370 |
+    +=======+========+========+====+========+========+
+    | A     | 1,2,3  | 1,2,3  |    | 1      | 1      |
+    +-------+--------+--------+----+--------+--------+
+    | B     | 1,2    | 2,3    |    | 1      | 2      |
+    +-------+--------+--------+----+--------+--------+
     """
     if dimensions is None:
         dimensions = {}

From 15b2f057bf88f8de6923f1bf5fa678dec0e4cbf6 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 21 Jun 2023 11:35:57 -0400
Subject: [PATCH 13/26] Update CHANGES.rst

Co-authored-by: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com>
---
 CHANGES.rst | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 08a5200c6..9235d2815 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -4,6 +4,10 @@ Changelog
 
 v0.45.0 (unreleased)
 --------------------
+Contributors to this version: David Huard (:user:`huard`).
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * ``ensembles.hawkins_sutton`` method to partition the uncertainty sources in a climate projection ensemble. (:issue:`771`, :pull:`1262`).
 
 v0.44.0 (unreleased)

From 0f540c78594b9b0d97d48500b3736e6b38f626a1 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 21 Jun 2023 12:05:45 -0400
Subject: [PATCH 14/26] table rst markup fix

---
 xclim/ensembles/_filters.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
index ff47bbcc3..9c7f51330 100644
--- a/xclim/ensembles/_filters.py
+++ b/xclim/ensembles/_filters.py
@@ -26,20 +26,20 @@ def _concat_hist(da, **hist):
         +------------+----------------------------------+
         | scenario   | time                             |
         +============+==================================+
-        | historical | hhhhhhhhhhhhhhhh                 |
+        | historical | ``hhhhhhhhhhhhhhhh------------`` |
         +------------+----------------------------------+
-        | ssp245     |                 111111111111     |
+        | ssp245     | ``----------------111111111111`` |
         +------------+----------------------------------+
-        | ssp370     |                 222222222222     |
+        | ssp370     | ``----------------222222222222`` |
         +------------+----------------------------------+
 
     to:
         +----------+----------------------------------+
         | scenario | time                             |
         +==========+==================================+
-        | ssp245   | hhhhhhhhhhhhhhhh111111111111     |
+        | ssp245   | ``hhhhhhhhhhhhhhhh111111111111`` |
         +----------+----------------------------------+
-        | ssp370   | hhhhhhhhhhhhhhhh222222222222     |
+        | ssp370   | ``hhhhhhhhhhhhhhhh222222222222`` |
         +----------+----------------------------------+
     """
     if len(hist) > 1:

From e9585c3e5ce28d683ed2b18c3ea3128c8ebc5c3e Mon Sep 17 00:00:00 2001
From: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com>
Date: Wed, 21 Jun 2023 12:45:59 -0400
Subject: [PATCH 15/26] add newline for formatting

---
 xclim/ensembles/_filters.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/xclim/ensembles/_filters.py b/xclim/ensembles/_filters.py
index 9c7f51330..ef0e8f969 100644
--- a/xclim/ensembles/_filters.py
+++ b/xclim/ensembles/_filters.py
@@ -34,6 +34,7 @@ def _concat_hist(da, **hist):
         +------------+----------------------------------+
 
     to:
+
         +----------+----------------------------------+
         | scenario | time                             |
         +==========+==================================+

From ecc62886379af576be850a6ab78e7779fecdd4d7 Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Tue, 27 Jun 2023 13:00:04 -0400
Subject: [PATCH 16/26] increase relative tolerance for error threshold

---
 tests/test_stats.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_stats.py b/tests/test_stats.py
index 07a44a892..da4e21999 100644
--- a/tests/test_stats.py
+++ b/tests/test_stats.py
@@ -141,7 +141,7 @@ def test_fit(self, fitda):
         # Test with MM
         pm = stats.fit(fitda, "lognorm", method="MM")
         mm, mv = lognorm(*pm.values).stats()
-        np.testing.assert_allclose(np.exp(2 + 1 / 2), mm, rtol=0.5)
+        np.testing.assert_allclose(np.exp(2 + 1 / 2), mm, rtol=0.65)
 
 
 def test_weibull_min_fit(weibull_min):

From 4c1c7f7d02d70711df0d1925e4521918b08b9542 Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Tue, 27 Jun 2023 14:16:09 -0400
Subject: [PATCH 17/26] add test for unsupported methods

---
 tests/test_analog.py | 13 +++++++++++++
 xclim/analog.py      |  6 +++---
 2 files changed, 16 insertions(+), 3 deletions(-)

diff --git a/tests/test_analog.py b/tests/test_analog.py
index ec6be6d53..48468eea4 100644
--- a/tests/test_analog.py
+++ b/tests/test_analog.py
@@ -79,6 +79,19 @@ def test_spatial_analogs(method, open_dataset):
     np.testing.assert_allclose(diss[method], out, rtol=1e-3, atol=1e-3)
 
 
+def test_unsupported_spatial_analog_method(open_dataset):
+    method = "KonMari"
+
+    data = open_dataset("SpatialAnalogs/indicators")
+    target = data.sel(lat=46.1875, lon=-72.1875, time=slice("1970", "1990"))
+    candidates = data.sel(time=slice("1970", "1990"))
+
+    match_statement = f"Method `KonMari` is not implemented. Available methods are: {','.join(xca.metrics.keys())}"
+
+    with pytest.raises(ValueError, match=match_statement):
+        xca.spatial_analogs(target, candidates, method=method)
+
+
 def test_spatial_analogs_multi_index(open_dataset):
     # Test multi-indexes
     diss = open_dataset("SpatialAnalogs/dissimilarity")
diff --git a/xclim/analog.py b/xclim/analog.py
index ea1b9f47f..b2807a9c8 100644
--- a/xclim/analog.py
+++ b/xclim/analog.py
@@ -78,10 +78,10 @@ def spatial_analogs(
         )
 
     try:
-        metric = metrics[method]
+        metric_func = metrics[method]
     except KeyError as e:
         raise ValueError(
-            f"Method {method} is not implemented. Available methods are : {','.join(metrics.keys())}."
+            f"Method `{method}` is not implemented. Available methods are: {','.join(metrics.keys())}."
         ) from e
 
     if candidates.chunks is not None:
@@ -91,7 +91,7 @@ def spatial_analogs(
 
     # Compute dissimilarity
     diss = xr.apply_ufunc(
-        metric,
+        metric_func,
         target,
         candidates,
         input_core_dims=[(dist_dim, "_indices"), ("_dist_dim", "_indices")],

From c4535774d320ab35f3ebfa8782b125c8591fda92 Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Tue, 27 Jun 2023 15:26:09 -0400
Subject: [PATCH 18/26] add tests for uncovered functions

---
 tests/test_sdba/test_processing.py | 44 ++++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/tests/test_sdba/test_processing.py b/tests/test_sdba/test_processing.py
index 7ea8633d9..f46e9ff6f 100644
--- a/tests/test_sdba/test_processing.py
+++ b/tests/test_sdba/test_processing.py
@@ -10,6 +10,7 @@
 from xclim.sdba.base import Grouper
 from xclim.sdba.processing import (
     adapt_freq,
+    construct_moving_yearly_window,
     escore,
     from_additive_space,
     jitter,
@@ -20,6 +21,7 @@
     stack_variables,
     standardize,
     to_additive_space,
+    unpack_moving_yearly_window,
     unstack_variables,
     unstandardize,
 )
@@ -280,3 +282,45 @@ def test_stack_variables(open_dataset):
     ds1p = unstack_variables(da1)
 
     xr.testing.assert_equal(ds1, ds1p)
+
+
+@pytest.mark.parametrize(
+    "window,step,lengths",
+    [
+        (1, 1, 151),
+        (5, 5, 30),
+        (10, 10, 15),
+        (25, 25, 6),
+        (50, 50, 3),
+        (None, None, 131),
+    ],
+)
+def test_construct_moving_yearly_window(open_dataset, window, step, lengths):
+    ds = open_dataset("sdba/CanESM2_1950-2100.nc")
+
+    calls = {k: v for k, v in dict(window=window, step=step).items() if v is not None}
+    da_windowed = construct_moving_yearly_window(ds.tasmax, **calls)
+
+    assert len(da_windowed) == lengths
+
+
+def test_construct_moving_yearly_window_standard_calendar(tasmin_series):
+    tasmin = tasmin_series(np.zeros(365 * 30), start="1997-01-01", units="degC")
+
+    with pytest.raises(ValueError):
+        construct_moving_yearly_window(tasmin)
+
+
+@pytest.mark.parametrize("append_ends", [True, False])
+def test_unpack_moving_yearly_window(open_dataset, append_ends):
+    tasmax = open_dataset("sdba/ahccd_1950-2013.nc").tasmax
+
+    tasmax_windowed = construct_moving_yearly_window(tasmax)
+
+    tx_deconstructed = unpack_moving_yearly_window(
+        tasmax_windowed, append_ends=append_ends
+    )
+    if append_ends:
+        np.testing.assert_array_equal(tasmax, tx_deconstructed)
+    else:
+        assert len(tx_deconstructed.time) < len(tasmax.time)

From 068f205b8e629af4163e6970f31b17a4780a4e8c Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Tue, 27 Jun 2023 15:49:58 -0400
Subject: [PATCH 19/26] add test for unused fixture options

---
 tests/test_testing_utils.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/tests/test_testing_utils.py b/tests/test_testing_utils.py
index a1b1de217..8fa079fef 100644
--- a/tests/test_testing_utils.py
+++ b/tests/test_testing_utils.py
@@ -6,9 +6,25 @@
 
 import numpy as np
 import pytest
+from xarray import Dataset
 
 import xclim.testing.utils as utilities
 from xclim import __version__ as __xclim_version__
+from xclim.testing.helpers import test_timeseries
+
+
+class TestFixtures:
+    def test_made_up_variable(self):
+        ds = test_timeseries(
+            np.zeros(365),
+            "luminiferous_aether_flux",
+            units="W K mol A-1 m-2 s-1",
+            as_dataset=True,
+        )
+
+        assert isinstance(ds, Dataset)
+        assert ds.luminiferous_aether_flux.attrs["units"] == "W K mol A-1 m-2 s-1"
+        assert "standard_name" not in ds.luminiferous_aether_flux.attrs
 
 
 class TestFileRequests:

From d146f33cce4df2f24d1ee33ca9c4a79e7c17db0a Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:05:27 -0400
Subject: [PATCH 20/26] rename test_timeseries on import

---
 tests/test_testing_utils.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/test_testing_utils.py b/tests/test_testing_utils.py
index 8fa079fef..767a79bd3 100644
--- a/tests/test_testing_utils.py
+++ b/tests/test_testing_utils.py
@@ -10,13 +10,13 @@
 
 import xclim.testing.utils as utilities
 from xclim import __version__ as __xclim_version__
-from xclim.testing.helpers import test_timeseries
+from xclim.testing.helpers import test_timeseries as timeseries
 
 
 class TestFixtures:
-    def test_made_up_variable(self):
-        ds = test_timeseries(
-            np.zeros(365),
+    def test_timeseries_made_up_variable(self):
+        ds = timeseries(
+            np.zeros(31),
             "luminiferous_aether_flux",
             units="W K mol A-1 m-2 s-1",
             as_dataset=True,

From 8b5672d784a3e143fc6456d9a8c73a9b7632b1c8 Mon Sep 17 00:00:00 2001
From: Zeitsperre <10819524+Zeitsperre@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:30:35 -0400
Subject: [PATCH 21/26] [skip ci] update CHANGES.rst

---
 CHANGES.rst | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 65921baca..4e72f9e79 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -2,6 +2,14 @@
 Changelog
 =========
 
+v0.45.0 (unreleased)
+--------------------
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`).
+
+Internal changes
+^^^^^^^^^^^^^^^^
+* Tolerance thresholds for error in ``test_stats::test_fit`` have been relaxed to allow for more variation in the results. Previously untested ``*_moving_yearly_window`` functions are now tested. (:issue:`1400`, :pull:`1402`).
+
 v0.44.0 (2023-06-23)
 --------------------
 Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Ludwig Lierhammer (:user:`ludwiglierhammer`), David Huard (:user:`huard`).

From 94e79102359d51d46913b4e5af06881a51f4ad56 Mon Sep 17 00:00:00 2001
From: "bumpversion[bot]" <bumpversion[bot]@ouranos.ca>
Date: Wed, 28 Jun 2023 14:32:04 +0000
Subject: [PATCH 22/26] =?UTF-8?q?Bump=20version:=200.44.0=20=E2=86=92=200.?=
 =?UTF-8?q?44.1-beta?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 setup.cfg         | 2 +-
 xclim/__init__.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index 7f89011c0..66914a2bb 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.44.0
+current_version = 0.44.1-beta
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/xclim/__init__.py b/xclim/__init__.py
index c6dd869ff..7abe45cf9 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -11,7 +11,7 @@
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.44.0"
+__version__ = "0.44.1-beta"
 
 
 # Load official locales

From f4aed882dc3c6f0c1af5c1f5374e7c9796dcb225 Mon Sep 17 00:00:00 2001
From: "bumpversion[bot]" <bumpversion[bot]@ouranos.ca>
Date: Wed, 28 Jun 2023 15:47:11 +0000
Subject: [PATCH 23/26] =?UTF-8?q?Bump=20version:=200.44.1-beta=20=E2=86=92?=
 =?UTF-8?q?=200.44.2-beta?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 setup.cfg         | 2 +-
 xclim/__init__.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index 66914a2bb..162c574ae 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.44.1-beta
+current_version = 0.44.2-beta
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 7abe45cf9..51d7bd477 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -11,7 +11,7 @@
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.44.1-beta"
+__version__ = "0.44.2-beta"
 
 
 # Load official locales

From dd5e1f569cb9d799294805b6c99320c3c3ae6492 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 28 Jun 2023 14:45:07 -0400
Subject: [PATCH 24/26] remove reference to xclim-testdata hawkins-sutton
 branch

---
 tests/test_partitioning.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/test_partitioning.py b/tests/test_partitioning.py
index e5e9e7a58..ea1ce26c8 100644
--- a/tests/test_partitioning.py
+++ b/tests/test_partitioning.py
@@ -11,9 +11,7 @@ def test_hawkins_sutton_smoke(open_dataset):
     """Just a smoke test."""
     dims = {"run": "member", "scen": "scenario"}
     da = (
-        open_dataset(
-            "uncertainty_partitioning/cmip5_pr_global_mon.nc", branch="hawkins_sutton"
-        )
+        open_dataset("uncertainty_partitioning/cmip5_pr_global_mon.nc")
         .pr.sel(time=slice("1950", None))
         .rename(dims)
     )

From f3a02d02d9ca619369ed2a3008853ea5ada725a0 Mon Sep 17 00:00:00 2001
From: David Huard <huard.david@ouranos.ca>
Date: Wed, 28 Jun 2023 16:04:23 -0400
Subject: [PATCH 25/26] update xclim-testdata tag for hawkins-sutton

---
 .github/workflows/main.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 3df012b72..f4094dac7 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -20,7 +20,7 @@ on:
       - submitted
 
 env:
-  XCLIM_TESTDATA_BRANCH: v2023.6.8
+  XCLIM_TESTDATA_BRANCH: v2023.6.28
 
 jobs:
   black:

From 0f1f7b67e5c22a014c190aaede054069b4b53c1c Mon Sep 17 00:00:00 2001
From: "bumpversion[bot]" <bumpversion[bot]@ouranos.ca>
Date: Wed, 28 Jun 2023 21:33:38 +0000
Subject: [PATCH 26/26] =?UTF-8?q?Bump=20version:=200.44.2-beta=20=E2=86=92?=
 =?UTF-8?q?=200.44.3-beta?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 setup.cfg         | 2 +-
 xclim/__init__.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index 162c574ae..507b34c26 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.44.2-beta
+current_version = 0.44.3-beta
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 51d7bd477..2de116e3d 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -11,7 +11,7 @@
 
 __author__ = """Travis Logan"""
 __email__ = "logan.travis@ouranos.ca"
-__version__ = "0.44.2-beta"
+__version__ = "0.44.3-beta"
 
 
 # Load official locales