Repository: Ouranosinc/xclim

Issue:
Multivariate quantile mapping bias correction implementation
Implement https://link.springer.com/article/10.1007/s00382-017-3580-6


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index af3ab196..842e6529 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -9,9 +9,12 @@ New features and enhancements
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 * Automatic load of translations on import and possibility to pass translations for virtual modules.
 * New ``xclim.testing.list_datasets`` function listing all available test datasets in repo `xclim-testdata`.
-* New ``sdba.construct_moving_yearly_window`` and ``sdba.unpack_moving_yearly_window`` for moving window adjustments.
 * `spatial_analogs` accepts multi-indexes as the `dist_dim` parameter and will work with candidates and target arrays of different lengths.
 * `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 Cannon's version of Pitié's _N-dimensional probability density function transform_. Uses new ``sdba.utils.rand_rot_matrix``.
+* 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.
 
@@ -51,6 +54,7 @@ Internal Changes
 * Simplification of some yaml elements for virtual modules.
 * Allow injecting `freq` without the missing checks failing.
 
+
 0.27.0 (2021-05-28)
 -------------------
 
diff --git a/docs/notebooks/sdba-advanced.ipynb b/docs/notebooks/sdba-advanced.ipynb
index f5ff92d8..afd5b169 100644
--- a/docs/notebooks/sdba-advanced.ipynb
+++ b/docs/notebooks/sdba-advanced.ipynb
@@ -28,7 +28,9 @@
     "\n",
     "* Another massive bottleneck of parallelization of xarray is the thread-locking behaviour of some methods. It is quite difficult to isolate and avoid those locking, so one of the best workaround is to use Dask configurations with many _processes_ and few _threads_. The former do not share memory and thus are not impacted when the other is locked. However, this adds many memory transfer operations and, by experience, reduces dask's ability to parallelize some pipelines. \n",
     " \n",
-    "   This is done usually by creating a dask client with a large `n_workers` and a small `threads_per_worker`. The last section of this notebook goes into a bit more details by suggesting a template for large-scale bias-adjustment pipelines."
+    "   This is done usually by creating a dask client with a large `n_workers` and a small `threads_per_worker`. The last section of this notebook goes into a bit more details by suggesting a template for large-scale bias-adjustment pipelines.\n",
+    "\n",
+    "* Sometimes, datasets have auxiliary coordinates (for example : lat / lon in a rotated pole dataset). Xarray handles these variables as data variables and will **not** load them if dask is used. However, in some operations, xclim or xarray will trigger an access to those variables, triggering computations each time, since they are dask-backed. To avoid this behaviour, one can load the coordinates, or simply remove them from the inputs."
    ]
   },
   {
diff --git a/docs/notebooks/sdba.ipynb b/docs/notebooks/sdba.ipynb
index fdc716cb..8d9587f9 100644
--- a/docs/notebooks/sdba.ipynb
+++ b/docs/notebooks/sdba.ipynb
@@ -373,6 +373,230 @@
     "axT.set_title('Timeseries - Point 1')"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Fourth example : Multivariate bias-adjustment with multiple steps - Cannon 2018\n",
+    "\n",
+    "This section replicates the \"MBCn\" algorithm described by [Cannon (2018)](https://doi.org/10.1007/s00382-017-3580-6). The method relies on some univariate algorithm, an adaption of the N-pdf transform of [Pitié et al. (2005)](https://ieeexplore.ieee.org/document/1544887/) and a final reordering step.\n",
+    "\n",
+    "In the following, we use the AHCCD and CanESM2 data are reference and simulation and we correct both `pr` and `tasmax` together."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim.testing import open_dataset\n",
+    "from xclim.core.units import convert_units_to\n",
+    "\n",
+    "dref = open_dataset('sdba/ahccd_1950-2013.nc', chunks={'location': 1}, drop_variables=['lat', 'lon']).sel(time=slice('1981', '2010'))\n",
+    "dref = dref.assign(\n",
+    "    tasmax=convert_units_to(dref.tasmax, 'K'),\n",
+    "    pr=convert_units_to(dref.pr, 'kg m-2 s-1')\n",
+    ")\n",
+    "dsim = open_dataset('sdba/CanESM2_1950-2100.nc', chunks={'location': 1}, drop_variables=['lat', 'lon'])\n",
+    "\n",
+    "dhist = dsim.sel(time=slice('1981', '2010'))\n",
+    "dsim = dsim.sel(time=slice('2041', '2070'))\n",
+    "dref"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Perform an initial univariate adjustment."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# additive for tasmax\n",
+    "QDMtx = sdba.QuantileDeltaMapping(nquantiles=20, kind='+', group='time')\n",
+    "QDMtx.train(dref.tasmax, dhist.tasmax)\n",
+    "# Adjust both hist and sim, we'll feed both to the Npdf transform.\n",
+    "scenh_tx = QDMtx.adjust(dhist.tasmax)\n",
+    "scens_tx = QDMtx.adjust(dsim.tasmax)\n",
+    "\n",
+    "# remove == 0 values in pr:\n",
+    "dref['pr'] = sdba.processing.jitter_under_thresh(dref.pr, 1e-5)\n",
+    "dhist['pr'] = sdba.processing.jitter_under_thresh(dhist.pr, 1e-5)\n",
+    "dsim['pr'] = sdba.processing.jitter_under_thresh(dsim.pr, 1e-5)\n",
+    "\n",
+    "# multiplicative for pr\n",
+    "QDMpr = sdba.QuantileDeltaMapping(nquantiles=20, kind='*', group='time')\n",
+    "QDMpr.train(dref.pr, dhist.pr)\n",
+    "# Adjust both hist and sim, we'll feed both to the Npdf transform.\n",
+    "scenh_pr = QDMpr.adjust(dhist.pr)\n",
+    "scens_pr = QDMpr.adjust(dsim.pr)\n",
+    "\n",
+    "scenh = xr.Dataset(dict(tasmax=scenh_tx, pr=scenh_pr))\n",
+    "scens = xr.Dataset(dict(tasmax=scens_tx, pr=scens_pr))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Stack the variables to multivariate arrays and standardize them\n",
+    "The standardization process ensure the mean and standard deviation of each column (variable) is 0 and 1 respectively.\n",
+    "\n",
+    "`hist` and `sim` are standardized together so the two series are coherent. We keep the mean and standard deviation to be reused when we build the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Stack the variables (tasmax and pr)\n",
+    "ref = sdba.base.stack_variables(dref)\n",
+    "scenh = sdba.base.stack_variables(scenh)\n",
+    "scens = sdba.base.stack_variables(scens)\n",
+    "\n",
+    "# Standardize\n",
+    "ref, _, _ = sdba.processing.standardize(ref)\n",
+    "\n",
+    "allsim, savg, sstd = sdba.processing.standardize(xr.concat((scenh, scens), 'time'))\n",
+    "hist = allsim.sel(time=scenh.time)\n",
+    "sim = allsim.sel(time=scens.time)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Perform the N-dimensional probability density function transform\n",
+    "\n",
+    "The NpdfTransform will iteratively randomly rotate our arrays in the \"variables\" space and apply the univariate adjustment before rotating it back. In Cannon (2018) and Pitié et al. (2005), it can be seen that the source array's joint distribution converges toward the target's joint distribution when a large number of iterations is done."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim import set_options\n",
+    "\n",
+    "NpdfT = sdba.adjustment.NpdfTransform(\n",
+    "    base=sdba.QuantileDeltaMapping,  # Use QDM as the univariate adjustment.\n",
+    "    base_kws={'nquantiles': 20, 'group': 'time'},\n",
+    "    n_iter=20,  # perform 20 iteration\n",
+    "    n_escore=1000,  # only send 1000 points to the escore metric (it is realy slow)\n",
+    ")\n",
+    "\n",
+    "# See the advanced notebook for details on how this option work\n",
+    "with set_options(sdba_extra_output=True):\n",
+    "    hist, sim, extra = NpdfT.train_adjust(ref, hist, sim)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Restoring the trend\n",
+    "\n",
+    "The NpdfT has given us new \"hist\" and \"sim\" arrays with a correct rank structure. However, the trend is lost in this process. We reorder the result of the initial adjustment according to the rank structure of the NpdfT outputs to get our final bias-adjusted series.\n",
+    "\n",
+    "`sdba.processing.reordering` is one of those function that need a dataset as input, instead of taking multiple arrays. The call sequence looks a bit clumsy: 'sim' is the argument to reorder, 'ref' the argument that provides the order."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scenh = sdba.processing.reordering(scenh, hist, group='time')\n",
+    "scens = sdba.processing.reordering(scens, sim, group='time')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scenh = sdba.base.unstack_variables(scenh, 'variables')\n",
+    "scens = sdba.base.unstack_variables(scens, 'variables')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### There we are!\n",
+    "\n",
+    "Let's trigger all the computations. Here we write the data to disk and use `compute=False` in order to trigger the whole computation tree only once. There seems to be no way in xarray to do the same with a `load` call."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from dask import compute\n",
+    "from dask.diagnostics import ProgressBar\n",
+    "\n",
+    "tasks = [\n",
+    "    scenh.isel(location=2).to_netcdf('mbcn_scen_hist_loc2.nc', compute=False),\n",
+    "    scens.isel(location=2).to_netcdf('mbcn_scen_sim_loc2.nc', compute=False),\n",
+    "    extra.escores.isel(location=2).to_dataset().to_netcdf('mbcn_escores_loc2.nc', compute=False)\n",
+    "]\n",
+    "\n",
+    "with ProgressBar():\n",
+    "    compute(tasks)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's compare the series and look at the distance scores to see how well the Npdf transform has converged."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "scenh = xr.open_dataset('mbcn_scen_hist_loc2.nc')\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "\n",
+    "dref.isel(location=2).tasmax.plot(ax=ax, label='Reference')\n",
+    "scenh.tasmax.plot(ax=ax, label='Adjusted', alpha=0.65)\n",
+    "dhist.isel(location=2).tasmax.plot(ax=ax, label='Simulated')\n",
+    "\n",
+    "ax.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "escores = xr.open_dataarray('mbcn_escores_loc2.nc')\n",
+    "diff_escore = escores.differentiate('iterations')\n",
+    "diff_escore.plot()\n",
+    "plt.title('Difference of the subsequent e-scores.')\n",
+    "plt.ylabel('E-scores difference')\n",
+    "\n",
+    "assert all(diff_escore < 0.2)  # this is for testing, please ignore"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
diff --git a/xclim/sdba/_adjustment.py b/xclim/sdba/_adjustment.py
index fbdb6f7f..20c157be 100644
--- a/xclim/sdba/_adjustment.py
+++ b/xclim/sdba/_adjustment.py
@@ -9,6 +9,7 @@ import xarray as xr
 from . import nbutils as nbu
 from . import utils as u
 from .base import Grouper, map_blocks, map_groups
+from .processing import escore
 
 
 @map_groups(
@@ -176,3 +177,82 @@ def scaling_adjust(ds, *, group, interp, kind):
     af = u.broadcast(ds.af, ds.sim, group=group, interp=interp)
     scen = u.apply_correction(ds.sim, af, kind)
     return scen.rename("scen").to_dataset()
+
+
+def npdf_transform(ds, **kwargs):
+    """N-pdf transform : Iterative univariate adjustment in random rotated spaces.
+
+    Parameters
+    ----------
+    Dataset variables:
+      ref : Reference multivariate timeseries
+      hist : simulated timeseries on the reference period
+      sim : Simulated timeseries on the projected period.
+      rot_matrices : Random rotation matrices.
+
+    kwargs:
+      pts_dim : multivariate dimensionn name
+      base : Adjustment class
+      base_kws : Kwargs for initialising the adjustment object
+      adj_kws : Kwargs of the `adjust` call
+      n_escore : Number of elements to include in the e_score test (0 for all, < 0 to skip)
+
+    Returns
+    -------
+    xr.Dataset
+      Dataset with `scenh`, `scens` and `escores` DataArrays, where `scenh` and `scens` are `hist` and `sim`
+      respectively after adjustment according to `ref`. If `n_escore` is negative, `escores` will be filled with NaNs.
+    """
+    ref = ds.ref
+    hist = ds.hist
+    sim = ds.sim.rename(time_sim="time")
+    dim = kwargs["pts_dim"]
+
+    escores = []
+    for i, R in enumerate(ds.rot_matrices.transpose("iterations", ...)):
+        # @ operator stands for matrix multiplication (along named dimensions): x@R = R@x
+        # @R rotates an array defined over dimension x unto new dimension x'. x@R = x'
+        refp = ref @ R
+        histp = hist @ R
+        simp = sim @ R
+
+        # Perform univariate adjustment in rotated space (x')
+        ADJ = kwargs["base"](**kwargs["base_kws"])
+        ADJ.train(refp, histp)
+        scenhp = ADJ.adjust(histp, **kwargs["adj_kws"])
+        scensp = ADJ.adjust(simp, **kwargs["adj_kws"])
+
+        # Rotate back to original dimension x'@R = x
+        # Note that x'@R is a back rotation because the matrix multiplication is now done along x' due to xarray
+        # operating along named dimensions.
+        # In normal linear algebra, this is equivalent to taking @R.T, the back rotation.
+        hist = scenhp @ R
+        sim = scensp @ R
+
+        # Compute score
+        if kwargs["n_escore"] >= 0:
+            escores.append(
+                escore(
+                    ref,
+                    hist,
+                    dims=(dim, "time"),
+                    N=kwargs["n_escore"],
+                    scale=True,
+                ).expand_dims(iterations=[i])
+            )
+
+    if kwargs["n_escore"] >= 0:
+        escores = xr.concat(escores, "iterations")
+    else:
+        # All NaN, but with the proper shape.
+        escores = (
+            ref.isel({dim: 0, "time": 0}) * hist.isel({dim: 0, "time": 0})
+        ).expand_dims(iterations=ds.iteration) * np.NaN
+
+    return xr.Dataset(
+        data_vars={
+            "scenh": hist.transpose(*ds.hist.dims),
+            "scens": sim.rename(time="time_sim").transpose(*ds.sim.dims),
+            "escores": escores,
+        }
+    )
diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py
index 5f1c2617..0d712368 100644
--- a/xclim/sdba/adjustment.py
+++ b/xclim/sdba/adjustment.py
@@ -1,5 +1,5 @@
 """Adjustment objects."""
-from typing import Optional, Sequence, Union
+from typing import Any, Mapping, Optional, Sequence, Union
 from warnings import warn
 
 import numpy as np
@@ -8,8 +8,9 @@ from xarray.core.dataarray import DataArray
 
 from xclim.core.calendar import get_calendar
 from xclim.core.formatting import update_history
-from xclim.core.options import OPTIONS, SDBA_EXTRA_OUTPUT
+from xclim.core.options import OPTIONS, SDBA_EXTRA_OUTPUT, set_options
 from xclim.core.units import convert_units_to
+from xclim.core.utils import uses_dask
 from xclim.indices import stats
 
 from ._adjustment import (
@@ -18,12 +19,13 @@ from ._adjustment import (
     eqm_train,
     loci_adjust,
     loci_train,
+    npdf_transform,
     qdm_adjust,
     qm_adjust,
     scaling_adjust,
     scaling_train,
 )
-from .base import Grouper, ParametrizableWithDataset, parse_group
+from .base import Grouper, Parametrizable, ParametrizableWithDataset, parse_group
 from .detrending import PolyDetrend
 from .utils import (
     ADDITIVE,
@@ -32,6 +34,7 @@ from .utils import (
     get_clusters,
     get_clusters_1d,
     pc_matrix,
+    rand_rot_matrix,
 )
 
 __all__ = [
@@ -43,6 +46,7 @@ __all__ = [
     "PrincipalComponents",
     "QuantileDeltaMapping",
     "Scaling",
+    "NpdfTransform",
 ]
 
 
@@ -997,3 +1001,185 @@ class PrincipalComponents(BaseAdjustment):
         )
         scen[lbl_R] = sim.indexes[lbl_M]
         return scen.unstack(lbl_R)
+
+
+class NpdfTransform(Parametrizable):
+    """N-dimensional probability density function transform."""
+
+    def __init__(
+        self,
+        base: BaseAdjustment = QuantileDeltaMapping,
+        base_kws: Optional[Mapping[str, Any]] = None,
+        n_escore: int = 0,
+        n_iter: int = 20,
+    ):
+        r"""
+        A multivariate bias-adjustment algorithm described by [Cannon18]_, as part of the MBCn algorithm,
+        based on a color-correction algorithm described by [Pitie05]_.
+
+        This algorithm in itself, when used with QuantileDeltaMapping, is NOT trend-preserving.
+        The full MBCn algorithm includes a reordering step provided here by :py:func:`xclim.sdba.processing.reordering`.
+
+        See notes for an explanation of the algorithm.
+
+        Parameters
+        ----------
+
+        At instantiation:
+
+        base: BaseAdjustment
+          An univariate bias-adjustment class. This is untested for anything else than QuantileDeltaMapping.
+        base_kws : dict, optional
+          Arguments passed to the initialization of the univariate adjustment.
+        n_escore : int
+          The number of elements to send to the escore function. The default, 0, means all elements are included.
+          Pass -1 to skip computing the escore completely.
+          Small numbers result in less significative scores, but the execution time goes up quickly with large values.
+        n_iter : int
+          The number of iterations to perform. Defaults to 20.
+
+        In train-adjustment:
+
+        pts_dim : str
+          The name of the "multivariate" dimension. Defaults to "variables", which is the
+          normal case when using :py:func:`xclim.sdba.base.stack_variables`.
+        adj_kws : dict, optional
+          Dictionary of arguments to pass to the adjust method of the univariate adjustment.
+        rot_matrices : xr.DataArray, optional
+          The rotation matrices as a 3D array ('iterations', <pts_dim>, <anything>), with shape (n_iter, <N>, <N>).
+          If left empty, random rotation matrices will be automatically generated.
+
+        Notes
+        -----
+        The historical reference (:math:`T`, for "target"), simulated historical (:math:`H`) and simulated projected (:math:`S`)
+        datasets are constructed by stacking the timeseries of N variables together. The algoriths goes into the
+        following steps:
+
+        1. Rotate the datasets in the N-dimensional variable space with :math:`\mathbf{R}`, a random rotation NxN matrix.
+
+        ..math ::
+
+            \tilde{\mathbf{T}} = \mathbf{T}\mathbf{R} \\
+            \tilde{\mathbf{H}} = \mathbf{H}\mathbf{R} \\
+            \tilde{\mathbf{S}} = \mathbf{S}\mathbf{R}
+
+        2. An univariate bias-adjustment :math:`\mathcal{F}` is used on the rotated datasets.
+        The adjustments are made in additive mode, for each variable :math:`i`.
+
+        .. math::
+
+            \hat{\mathbf{H}}_i, \hat{\mathbf{S}}_i = \mathcal{F}\left(\tilde{\mathbf{T}}_i, \tilde{\mathbf{H}}_i, \tilde{\mathbf{S}}_i\right)
+
+        3. The bias-adjusted datasets are rotated back.
+
+        .. math::
+
+            \mathbf{H}' = \hat{\mathbf{H}}\mathbf{R} \\
+            \mathbf{S}' = \hat{\mathbf{S}}\mathbf{R}
+
+
+        These three steps are repeated a certain number of times, prescribed by argument ``n_iter``. At each
+        iteration, a new random rotation matrix is generated.
+
+        The original algorithm ([Pitie05]_), stops the iteration when some distance score converges. Following
+        [Cannon18_] and the MBCn implementation in [CannonR]_, we instead fix the number of iterations.
+
+        As done by [Cannon18]_, the distance score chosen is the "Energy distance" from [SkezelyRizzo]_
+        (see :py:func:`xclim.sdba.processing.escore`).
+
+        The random matrices are generated following a method laid out by [Mezzadri].
+
+        This is only part of the full MBCn algorithm, see :ref:`The MBCn algorithm` for an example on how
+        to replicate the full method with xclim. This includes a standardization of the simulated data beforehand,
+        an initial univariate adjustment and the reordering of those adjusted series according to the
+        rank structure of the output of this algorithm.
+
+        References
+        ----------
+        .. [Cannon18] Cannon, A. J. (2018). Multivariate quantile mapping bias correction: An N-dimensional probability density function transform for climate model simulations of multiple variables. Climate Dynamics, 50(1), 31–49. https://doi.org/10.1007/s00382-017-3580-6
+        .. [Mezzadri]Mezzadri, F. (2006). How to generate random matrices from the classical compact groups. arXiv preprint math-ph/0609050.
+        .. [Pitie05] Pitie, F., Kokaram, A. C., & Dahyot, R. (2005). N-dimensional probability density function transfer and its application to color transfer. Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, 2, 1434-1439 Vol. 2. https://doi.org/10.1109/ICCV.2005.166
+        .. [SkezelyRizzo] Szekely, G. J. and Rizzo, M. L. (2004) Testing for Equal Distributions in High Dimension, InterStat, November (5)
+        """
+        base_kws or {}
+        if "kind" in base_kws:
+            warn(
+                f'The adjustment kind cannot be controlled when using {self.__class__.__name__}, it defaults to "+".'
+            )
+        base_kws.setdefault("kind", "+")
+
+        super().__init__(
+            base=base,
+            base_kws=base_kws,
+            n_escore=n_escore,
+            n_iter=n_iter,
+        )
+
+    def train_adjust(
+        self,
+        ref: xr.DataArray,
+        hist: xr.DataArray,
+        sim: xr.DataArray,
+        *,
+        pts_dim: str = "variables",
+        adj_kws: Optional[Mapping[str, Any]] = None,
+        rot_matrices: Optional[xr.DataArray] = None,
+    ):
+        # Assuming sim has the same coords as hist
+        # We get the safest new name of the rotated dim.
+        rot_dim = xr.core.utils.get_temp_dimname(
+            set(ref.dims).union(hist.dims).union(sim.dims), pts_dim + "_prime"
+        )
+
+        # Get the rotation matrices
+        rot_matrices = rot_matrices or rand_rot_matrix(
+            ref[pts_dim], num=self.n_iter, new_dim=rot_dim
+        ).rename(matrices="iterations")
+
+        # Call a map_blocks on the iterative function
+        # Sadly, this is a bit too complicated for map_blocks, we'll do it by hand.
+        escores_tmpl = xr.broadcast(
+            ref.isel({pts_dim: 0, "time": 0}),
+            hist.isel({pts_dim: 0, "time": 0}),
+        )[0].expand_dims(iterations=rot_matrices.iterations)
+
+        template = xr.Dataset(
+            data_vars={
+                "scenh": xr.full_like(hist, np.NaN),
+                "scens": xr.full_like(sim, np.NaN).rename(time="time_sim"),
+                "escores": escores_tmpl,
+            }
+        )
+
+        # Input data, rename time dim on sim since it can't be aligned with ref or hist.
+        ds = xr.Dataset(
+            data_vars={
+                "ref": ref,
+                "hist": hist,
+                "sim": sim.rename(time="time_sim"),
+                "rot_matrices": rot_matrices,
+            }
+        )
+
+        if uses_dask(ds) and any(
+            [len(ds.chunks.get(d, [])) > 1 for d in ["time", "time_sim", pts_dim]]
+        ):
+            raise ValueError(
+                f'Inputs of {self.__class__.__name__} cannot be chunked along the main dimensions "time" and "{pts_dim}"'
+            )
+
+        kwargs = self.parameters.copy()
+        kwargs.update(pts_dim=pts_dim, adj_kws=adj_kws or {})
+
+        with set_options(sdba_extra_output=False):
+            out = ds.map_blocks(npdf_transform, template=template, kwargs=kwargs)
+
+        scenh = out.scenh
+        scens = out.scens.rename(time_sim="time")
+
+        if OPTIONS[SDBA_EXTRA_OUTPUT]:
+            extra = xr.Dataset(
+                dict(escores=out.escores, rotation_matrices=rot_matrices)
+            )
+            return scenh, scens, extra
+        return scenh, scens
diff --git a/xclim/sdba/base.py b/xclim/sdba/base.py
index eb6e1034..6922e1bd 100644
--- a/xclim/sdba/base.py
+++ b/xclim/sdba/base.py
@@ -115,12 +115,13 @@ class Grouper(Parametrizable):
           If larger than 1, a centered rolling window along the main dimension is created when grouping data.
           Units are the sampling frequency of the data along the main dimension.
         add_dims : Optional[Sequence[str]]
-          Additionnal dimensions that should be reduced in grouping operations. This behaviour is also controlled
+          Additional dimensions that should be reduced in grouping operations. This behaviour is also controlled
           by the `main_only` parameter of the `apply` method. If any of these dimensions are absent from the dataarrays,
           they will be omitted.
         interp : Union[bool, str]
           Whether to return an interpolatable index in the `get_index` method. Only effective for `month` grouping.
-          Interpolation method names are accepted for convenience, "nearest" is translated to False, all other names are translated to True.
+          Interpolation method names are accepted for convenience, "nearest" is translated to False, all other names
+          are translated to True.
           This modifies the default, but `get_index` also accepts an `interp` argument overriding the one defined here..
         """
         if "." in group:
@@ -737,3 +738,83 @@ def unpack_moving_yearly_window(da: xr.DataArray, dim: str = "movingwin"):
         out.append(slc)
 
     return xr.concat(out, "time")
+
+
+def stack_variables(ds, rechunk=True, dim="variables"):
+    """Stack different variables of a dataset into a single DataArray with a new "variables" dimension.
+
+    Variable attributes are all added as lists of attributes.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+      Input dataset.
+    rechunk : bool
+      If True (default), dask arrays are rechunked with `variables : -1`.
+    dim : str
+      Name of dimension along which variables are indexed.
+
+    Returns
+    -------
+    xr.DataArray
+      Array with variables stacked along `dim` dimension.
+    """
+    # Store original arrays' attributes
+    attrs = {}
+    nvar = len(ds.data_vars)
+    for i, var in enumerate(ds.data_vars.values()):
+        for name, attr in var.attrs.items():
+            attrs.setdefault(name, [None] * nvar)[i] = attr
+
+    # Special key used for later `unstacking`
+    attrs["is_variables"] = True
+    var_crd = xr.DataArray(
+        list(ds.data_vars.keys()), dims=(dim,), name=dim, attrs=attrs
+    )
+
+    da = xr.concat(ds.data_vars.values(), var_crd, combine_attrs="drop")
+
+    if uses_dask(da) and rechunk:
+        da = da.chunk({dim: -1})
+
+    da.attrs.update(ds.attrs)
+    return da.rename("multivariate")
+
+
+def unstack_variables(da, dim=None):
+    """Unstack a DataArray created by `stack_variables` to a dataset.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      Array holding different variables along `dim` dimension.
+    dim : str
+      Name of dimension along which the variables are stacked. If not specified (default),
+      `dim` is inferred from attributes of the coordinate.
+
+    Returns
+    -------
+    xr.Dataset
+      Dataset holding each variable in an individual DataArray.
+    """
+    if dim is None:
+        for dim, crd in da.coords.items():
+            if crd.attrs.get("is_variables"):
+                break
+        else:
+            raise ValueError("No variable coordinate found, were attributes removed?")
+
+    ds = xr.Dataset(
+        {name.item(): da.sel({dim: name.item()}, drop=True) for name in da[dim]},
+        attrs=da.attrs,
+    )
+
+    # Reset attributes
+    for name, attr_list in da.variables.attrs.items():
+        if name == "is_variables":
+            continue
+        for attr, var in zip(attr_list, da.variables):
+            if attr is not None:
+                ds[var.item()].attrs[name] = attr
+
+    return ds
diff --git a/xclim/sdba/nbutils.py b/xclim/sdba/nbutils.py
index c84a273f..11b1939e 100644
--- a/xclim/sdba/nbutils.py
+++ b/xclim/sdba/nbutils.py
@@ -1,7 +1,7 @@
 """Numba-accelerated utils."""
 import numpy as np
 import xarray as xr
-from numba import float32, float64, guvectorize, jit  # , int32, int64
+from numba import boolean, float32, float64, guvectorize, njit
 
 
 @guvectorize(
@@ -32,14 +32,13 @@ def vecquantiles(da, rnk, dim):
     return res
 
 
-@jit(
+@njit(
     [
         float32[:, :](float32[:, :], float32[:]),
         float64[:, :](float64[:, :], float64[:]),
         float32[:](float32[:], float32[:]),
         float64[:](float64[:], float64[:]),
     ],
-    nopython=True,
 )
 def _quantile(arr, q):
     if arr.ndim == 1:
@@ -94,3 +93,75 @@ def quantile(da, q, dim):
         )
 
     return res
+
+
+@njit([float32[:, :](float32[:, :]), float64[:, :](float64[:, :])])
+def remove_NaNs(x):
+    remove = np.zeros_like(x[0, :], dtype=boolean)
+    for i in range(x.shape[0]):
+        remove = remove | np.isnan(x[i, :])
+    return x[:, ~remove]
+
+
+@njit([float32(float32[:]), float64(float64[:])], fastmath=True)
+def _euclidean_norm(v):
+    """Compute the euclidean norm of vector v."""
+    return np.sqrt(np.sum(v ** 2))
+
+
+@njit(
+    [float32(float32[:, :], float32[:, :]), float64(float64[:, :], float64[:, :])],
+    fastmath=True,
+)
+def _correlation(X, Y):
+    """Compute a correlation as the mean of pairwise distances between points in X and Y.
+
+    X is KxN and Y is KxM, the result is the mean of the MxN distances.
+    Similar to scipy.spatial.distance.cdist(X, Y, 'euclidean')
+    """
+    d = 0
+    for i in range(X.shape[1]):
+        for j in range(Y.shape[1]):
+            d += _euclidean_norm(X[:, i] - Y[:, j])
+    return d / (X.shape[1] * Y.shape[1])
+
+
+@njit([float32(float32[:, :]), float64(float64[:, :])], fastmath=True)
+def _autocorrelation(X):
+    """Mean of the NxN pairwise distances of points in X of shape KxN.
+
+    Similar to scipy.spatial.distance.pdist(..., 'euclidean')
+    """
+    d = 0
+    for i in range(X.shape[1]):
+        for j in range(i + 1):
+            d += (int(i != j) + 1) * _euclidean_norm(X[:, i] - X[:, j])
+    return d / X.shape[1] ** 2
+
+
+@guvectorize(
+    [
+        (float32[:, :], float32[:, :], float32[:]),
+        (float64[:, :], float64[:, :], float64[:]),
+    ],
+    "(k, n),(k, m)->()",
+)
+def _escore(tgt, sim, out):
+    """E-score based on the Skezely-Rizzo e-distances between clusters.
+
+    tgt and sim are KxN and KxM, where dimensions are along K and observations along M and N.
+    When N > 0, only this many points of target and sim are used, taken evenly distributed in the series.
+    When std is True, X and Y are standardized according to the nanmean and nanstd (ddof = 1) of X.
+    """
+    sim = remove_NaNs(sim)
+    tgt = remove_NaNs(tgt)
+
+    n1 = sim.shape[1]
+    n2 = tgt.shape[1]
+
+    sXY = _correlation(tgt, sim)
+    sXX = _autocorrelation(tgt)
+    sYY = _autocorrelation(sim)
+
+    w = n1 * n2 / (n1 + n2)
+    out[0] = w * (sXY + sXY - sXX - sYY) / 2
diff --git a/xclim/sdba/processing.py b/xclim/sdba/processing.py
index b0f66c52..22539a1c 100644
--- a/xclim/sdba/processing.py
+++ b/xclim/sdba/processing.py
@@ -1,12 +1,14 @@
 """Pre and post processing for bias adjustment."""
+from typing import Optional, Sequence
+
 import dask.array as dsk
 import numpy as np
 import xarray as xr
 
 from xclim.core.utils import uses_dask
 
+from . import nbutils as nbu
 from .base import Grouper, map_groups
-from .nbutils import vecquantiles
 from .utils import ADDITIVE, apply_correction, ecdf, invert
 
 
@@ -29,6 +31,8 @@ def adapt_freq(
     ds : xr.Dataset
       With variables :  "ref", Target/reference data, usually observed data.
       and  "sim", Simulated data.
+    dim : str
+      Dimension name.
     group : Union[str, Grouper]
       Grouping information, see base.Grouper
     thresh : float
@@ -71,7 +75,7 @@ def adapt_freq(
     # Compute : ecdf_ref^-1( ecdf_sim( thresh ) )
     # The value in ref with the same rank as the first non zero value in sim.
     # pth is meaningless when freq. adaptation is not needed
-    pth = vecquantiles(ds.ref, P0_sim, dim).where(dP0 > 0)
+    pth = nbu.vecquantiles(ds.ref, P0_sim, dim).where(dP0 > 0)
 
     if "window" in ds.sim.dims:
         # P0_sim was computed using the window, but only the original timeseries is corrected.
@@ -137,7 +141,7 @@ def jitter_under_thresh(x: xr.DataArray, thresh: float):
         )
     else:
         jitter = np.random.uniform(low=epsilon, high=thresh, size=x.shape)
-    return x.where(~((x < thresh) & (x.notnull())), jitter)
+    return x.where(~((x < thresh) & (x.notnull())), jitter.astype(x.dtype))
 
 
 def jitter_over_thresh(x: xr.DataArray, thresh: float, upper_bnd: float):
@@ -167,7 +171,7 @@ def jitter_over_thresh(x: xr.DataArray, thresh: float, upper_bnd: float):
         )
     else:
         jitter = np.random.uniform(low=thresh, high=upper_bnd, size=x.shape)
-    return x.where(~((x > thresh) & (x.notnull())), jitter)
+    return x.where(~((x > thresh) & (x.notnull())), jitter.astype(x.dtype))
 
 
 @map_groups(reduces=[Grouper.PROP], data=[])
@@ -219,3 +223,167 @@ def uniform_noise_like(da: xr.DataArray, low: float = 1e-6, high: float = 1e-3):
     return da.copy(
         data=(high - low) * mod.random.random_sample(size=da.shape, **kw) + low
     )
+
+
+def standardize(
+    da: xr.DataArray,
+    mean: Optional[xr.DataArray] = None,
+    std: Optional[xr.DataArray] = None,
+    dim: str = "time",
+):
+    """Standardize a DataArray by centering its mean and scaling it by its standard deviation.
+
+    Either of both of mean and std can be provided if need be.
+
+    Returns the standardized data, the mean and the standard deviation.
+    """
+    if mean is None:
+        mean = da.mean(dim)
+    if std is None:
+        std = da.std(dim)
+    with xr.set_options(keep_attrs=True):
+        return (da - mean) / std, mean, std
+
+
+def unstandardize(da: xr.DataArray, mean: xr.DataArray, std: xr.DataArray):
+    """Rescale a standardized array by performing the inverse operation of `standardize`."""
+    return (std * da) + mean
+
+
+@map_groups(reordered=[Grouper.DIM], main_only=True)
+def _reordering_group(ds, *, dim):
+    """Group-wise reordering."""
+
+    def _reordering_1d(data, ordr):
+        return np.sort(data)[np.argsort(np.argsort(ordr))]
+
+    return (
+        xr.apply_ufunc(
+            _reordering_1d,
+            ds.sim,
+            ds.ref,
+            input_core_dims=[[dim], [dim]],
+            output_core_dims=[[dim]],
+            vectorize=True,
+            dask="parallelized",
+            output_dtypes=[ds.sim.dtype],
+        )
+        .rename("reordered")
+        .to_dataset()
+    )
+
+
+def reordering(sim, ref, group="time"):
+    """Reorders data in `sim` following the order of ref.
+
+    The rank structure of `ref` is used to reorder the elements of `sim` along dimension "time",
+    optionally doing the operation group-wise.
+
+    Parameters
+    ----------
+    sim : xr.DataArray
+      Array to reorder.
+    ref : xr.DataArray
+      Array whose rank order sim should replicate.
+    group : str
+      Grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
+
+    Returns
+    -------
+    xr.Dataset
+      sim reordered according to ref's rank order.
+
+    Reference
+    ---------
+    Cannon, A. J. (2018). Multivariate quantile mapping bias correction: An N-dimensional probability density function
+    transform for climate model simulations of multiple variables. Climate Dynamics, 50(1), 31–49.
+    https://doi.org/10.1007/s00382-017-3580-6
+    """
+    ds = xr.Dataset({"sim": sim, "ref": ref})
+    return _reordering_group(ds, group=group).reordered
+
+
+def escore(
+    tgt: xr.DataArray,
+    sim: xr.DataArray,
+    dims: Sequence[str] = ("variables", "time"),
+    N: int = 0,
+    scale: bool = False,
+):
+    r"""Energy score, or energy dissimilarity metric, based on [SkezelyRizzo]_ and [Cannon18]_.
+
+    Parameters
+    ----------
+    tgt: DataArray
+      Target observations.
+    sim: DataArray
+      Candidate observations. Must have the same dimensions as `tgt`.
+    dims: sequence of 2 strings
+      The name of the dimensions along which the variables and observation points are listed.
+      `tgt` and `sim` can have different length along the second one, but must be equal along the first one.
+      The result will keep all other dimensions.
+    N : int
+      If larger than 0, the number of observations to use in the score computation. The points are taken
+      evenly distributed along `obs_dim`.
+    scale: boolean
+      Whether to scale the data before computing the score. If True, both arrays as scaled according
+      to the mean and standard deviation of `tgt` along `obs_dim`. (std computed with `ddof=1` and both
+      statistics excluding NaN values.
+
+    Returns
+    -------
+    e-score
+        xr.DataArray with dimensions not in `dims`.
+
+    Notes
+    -----
+    Explanation adapted from the "energy" R package documentation.
+    The e-distance between two clusters :math:`C_i`, :math:`C_j` (tgt and sim) of size :math:`n_i,,n_j`
+    proposed by Szekely and Rizzo (2005) is defined by:
+
+    .. math::
+
+        e(C_i,C_j) = \frac{1}{2}\frac{n_i n_j}{n_i + n_j} \left[2 M_{ij} − M_{ii} − M_{jj}\right]
+
+    where
+
+    .. math::
+
+        M_{ij} = \frac{1}{n_i n_j} \sum_{p = 1}^{n_i} \sum{q = 1}^{n_j} \left\Vert X_{ip} − X{jq} \right\Vert.
+
+    :math:`\Vert\cdot\Vert` denotes Euclidean norm, :math:`X_{ip}` denotes the p-th observation in the i-th cluster.
+
+    The input scaling and the factor :math:`\frac{1}{2}` in the first equation are additions of [Cannon18]_ to
+    the metric. With that factor, the test becomes identical to the one defined by [BaringhausFranz]_.
+
+    References
+    ----------
+    .. [SkezelyRizzo] Szekely, G. J. and Rizzo, M. L. (2004) Testing for Equal Distributions in High Dimension, InterStat, November (5)
+    .. [BaringhausFranz] Baringhaus, L. and Franz, C. (2004) On a new multivariate two-sample test, Journal of Multivariate Analysis, 88(1), 190–206. https://doi.org/10.1016/s0047-259x(03)00079-4
+    """
+
+    pts_dim, obs_dim = dims
+
+    if N > 0:
+        # If N non-zero we only take around N points, evenly distributed
+        sim_step = int(np.ceil(sim[obs_dim].size / N))
+        sim = sim.isel({obs_dim: slice(None, None, sim_step)})
+        tgt_step = int(np.ceil(tgt[obs_dim].size / N))
+        tgt = tgt.isel({obs_dim: slice(None, None, tgt_step)})
+
+    if scale:
+        tgt, avg, std = standardize(tgt)
+        sim, _, _ = standardize(sim, avg, std)
+
+    # The dimension renaming is to allow different coordinates.
+    # Otherwise, apply_ufunc tries to align both obs_dim together.
+    new_dim = xr.core.utils.get_temp_dimname(tgt.dims, obs_dim)
+    sim = sim.rename({obs_dim: new_dim})
+    return xr.apply_ufunc(
+        nbu._escore,
+        tgt,
+        sim,
+        input_core_dims=[[pts_dim, obs_dim], [pts_dim, new_dim]],
+        output_dtypes=[sim.dtype],
+        dask="parallelized",
+    )
diff --git a/xclim/sdba/utils.py b/xclim/sdba/utils.py
index ad281cb6..8fe29ca5 100644
--- a/xclim/sdba/utils.py
+++ b/xclim/sdba/utils.py
@@ -685,3 +685,53 @@ def get_clusters(data: xr.DataArray, u1, u2, dim: str = "time"):
     )
 
     return ds
+
+
+def rand_rot_matrix(crd: xr.DataArray, num: int = 1, new_dim: Optional[str] = None):
+    r"""Generate random rotation matrices.
+
+    Rotation matrices are members of the SO(n) group, where n is the matrix size (`crd.size`).
+    They can be characterized as orthogonal matrices with determinant 1. A square matrix :math:`R`
+    is a rotation matrix if and only if :math:`R^t = R^{−1}` and :math:`\mathrm{det} R = 1`.
+
+    Parameters
+    ----------
+    crd: DataArray
+      1D coordinate DataArray along which the rotation occurs.
+      The output will be square with the same coordinate replicated,
+      the second renamed to `new_dim`.
+    num : int
+      If larger than 1 (default), the number of matrices to generate, stacked along a "matrices" dimension.
+    new_dim : str
+      Name of the new "prime" dimension, defaults to the same name as `crd` + "_prime".
+
+    Returns
+    -------
+    DataArray
+      float, NxN if num = 1, numxNxN otherwise, where N is the length of crd.
+
+    References
+    ----------
+    .. [Mezzadri] Mezzadri, F. (2006). How to generate random matrices from the classical compact groups. arXiv preprint math-ph/0609050.
+    """
+    if num > 1:
+        return xr.concat([rand_rot_matrix(crd, num=1) for i in range(num)], "matrices")
+
+    N = crd.size
+    dim = crd.dims[0]
+    # Rename and rebuild second coordinate : "prime" axis.
+    if new_dim is None:
+        new_dim = dim + "_prime"
+    crd2 = xr.DataArray(crd.values, dims=new_dim, name=new_dim, attrs=crd.attrs)
+
+    # Random floats from the standardized normal distribution
+    Z = np.random.standard_normal((N, N))
+
+    # QR decomposition and manipulation from Mezzadri 2006
+    Q, R = np.linalg.qr(Z)
+    num = np.diag(R)
+    denum = np.abs(num)
+    lam = np.diag(num / denum)  # "lambda"
+    return xr.DataArray(
+        Q @ lam, dims=(dim, new_dim), coords={dim: crd, new_dim: crd2}
+    ).astype("float32")


Test patch:
diff --git a/xclim/testing/tests/test_sdba/test_adjustment.py b/xclim/testing/tests/test_sdba/test_adjustment.py
index 5667a4d7..fa8fde42 100644
--- a/xclim/testing/tests/test_sdba/test_adjustment.py
+++ b/xclim/testing/tests/test_sdba/test_adjustment.py
@@ -11,12 +11,18 @@ from xclim.sdba.adjustment import (
     DetrendedQuantileMapping,
     EmpiricalQuantileMapping,
     ExtremeValues,
+    NpdfTransform,
     PrincipalComponents,
     QuantileDeltaMapping,
     Scaling,
 )
-from xclim.sdba.base import Grouper
-from xclim.sdba.processing import jitter_under_thresh, uniform_noise_like
+from xclim.sdba.base import Grouper, stack_variables, unstack_variables
+from xclim.sdba.processing import (
+    jitter_under_thresh,
+    reordering,
+    standardize,
+    uniform_noise_like,
+)
 from xclim.sdba.utils import (
     ADDITIVE,
     MULTIPLICATIVE,
diff --git a/xclim/testing/tests/test_sdba/test_processing.py b/xclim/testing/tests/test_sdba/test_processing.py
index 9c6258a6..7ddd5502 100644
--- a/xclim/testing/tests/test_sdba/test_processing.py
+++ b/xclim/testing/tests/test_sdba/test_processing.py
@@ -4,7 +4,15 @@ import pytest
 import xarray as xr
 
 from xclim.sdba.base import Grouper
-from xclim.sdba.processing import adapt_freq, jitter_over_thresh, jitter_under_thresh
+from xclim.sdba.processing import (
+    adapt_freq,
+    escore,
+    jitter_over_thresh,
+    jitter_under_thresh,
+    reordering,
+    standardize,
+    unstandardize,
+)
 
 
 def test_jitter_under_thresh():
@@ -42,7 +50,9 @@ def test_adapt_freq(use_dask):
     prsim = xr.where(pr < 20, pr / 20, pr)
     prref = xr.where(pr < 10, pr / 20, pr)
     ds_in = xr.Dataset({"sim": prsim, "ref": prref})
-    ds_ad = adapt_freq(ds_in, thresh=1, group=group)
+    ds_ad = adapt_freq(
+        ds_in, thresh=1, group=group
+    )  # Don't worry about CodeFactor complaint. dim is fed by decorator
 
     # Where the input is considered zero
     input_zeros = ds_ad.sim_ad.where(prsim <= 1)
@@ -69,3 +79,42 @@ def test_adapt_freq(use_dask):
     # Assert that Pth and dP0 are approx the good values
     np.testing.assert_allclose(ds_ad.pth, 20, rtol=0.05)
     np.testing.assert_allclose(ds_ad.dP0, 0.5, atol=0.14)
+
+
+def test_escore():
+
+    x = np.array([1, 4, 3, 6, 4, 7, 5, 8, 4, 5, 3, 7]).reshape(2, 6)
+    y = np.array([6, 6, 3, 8, 5, 7, 3, 7, 3, 6, 4, 3]).reshape(2, 6)
+
+    x = xr.DataArray(x, dims=("variables", "time"))
+    y = xr.DataArray(y, dims=("variables", "time"))
+
+    # Value taken from escore of Cannon's MBC R package.
+    np.testing.assert_allclose(escore(x, y), 1.90018550338863)
+
+
+def test_standardize():
+    x = np.random.standard_normal((2, 10000))
+    x[0, 50] = np.NaN
+    x = xr.DataArray(x, dims=("x", "y"))
+
+    xp, avg, std = standardize(x, dim="y")
+
+    np.testing.assert_allclose(avg, 0, atol=3e-2)
+    np.testing.assert_allclose(std, 1, atol=2e-2)
+
+    xp, avg, std = standardize(x, mean=avg, dim="y")
+    np.testing.assert_allclose(std, 1, atol=2e-2)
+
+    y = unstandardize(xp, 0, 1)
+
+    np.testing.assert_allclose(x, y, atol=0.1)
+
+
+def test_reordering():
+    x = xr.DataArray(np.arange(1, 11), dims=("time",))
+    y = xr.DataArray(np.arange(10, 20)[::-1], dims=("time",))
+
+    out = reordering(x, y, group="time")
+
+    np.testing.assert_array_equal(out, np.arange(1, 11)[::-1])