Repository: ouranosinc/xclim
Issue: ANUCLIM calls to lazy_indexing failing

* xclim version: '0.17.1-beta'
* Python version:Python 3.7.6
* Operating System: Linux

### Description
Calculating anuclim indices on real data failing for those relying on `lazy_indexing`,   
### What I Did
<!--Paste the command(s) you ran and the output.
If there was a crash, please include the traceback below.-->
```
from xclim import anuclim as anu
from xclim import indices as xci
import xarray as xr
from pathlib import Path
from xclim import subset
source = Path('data_monthly')
ncfiles = []
for v in ['tasmax','tasmin','pr']:
    ncfiles.append(sorted(list(source.rglob(f'*{v}*.nc'))))
ds = xr.open_mfdataset(ncfiles, chunks=dict(time=12), combine='by_coords')
ds =subset.subset_bbox(ds, lon_bnds=[-120, -60], lat_bnds=[47, 79], start_date='1955')
ds.pr.attrs['units'] = 'mm month-1'
ds['tas'] = 0.5*(ds.tasmax + ds.tasmin)
ds['tas'].attrs = ds.tasmax.attrs
ds['tas'].attrs['long_name'] = 'Monthly average temperature'
dsOut['P8_MeanTempWettestQuarter'] = anu.P8_MeanTempWettestQuarter(tas=ds.tas,
                                                                   pr=ds.pr,
                                                                   input_freq='monthly', freq='YS')
```

### What I Received
<!--Paste the output or the stack trace of the problem you experienced here.-->
```
....

Traceback (most recent call last):
  File "C:\\Program Files\\JetBrains\\PyCharm Community Edition 2019.1\\helpers\\pydev\\pydevd.py", line 1741, in <module>
    main()
  File "C:\\Program Files\\JetBrains\\PyCharm Community Edition 2019.1\\helpers\\pydev\\pydevd.py", line 1735, in main
    globals = debugger.run(setup['file'], None, None, is_module)
  File "C:\\Program Files\\JetBrains\\PyCharm Community Edition 2019.1\\helpers\\pydev\\pydevd.py", line 1135, in run
    pydev_imports.execfile(file, globals, locals)  # execute the script
  File "C:\\Program Files\\JetBrains\\PyCharm Community Edition 2019.1\\helpers\\pydev\\_pydev_imps\\_pydev_execfile.py", line 18, in execfile
    exec(compile(contents+"\\n", file, 'exec'), glob, loc)
  File "W:/Projets/Biodiversite/Humphries/Climate_indices/Calculate_climate_indices.py", line 48, in <module>
    input_freq='monthly', freq='YS')
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indicators\\__init__.py", line 179, in _wrapper
    return func(*args, **kwargs)
  File "<boltons.funcutils.FunctionBuilder-93>", line 2, in tg_mean_wetdry_quarter
  File "<boltons.funcutils.FunctionBuilder-44>", line 2, in tg_mean_wetdry_quarter
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\core\\units.py", line 361, in wrapper
    out = func(*args, **kwargs)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indices\\_anuclim.py", line 287, in tg_mean_wetdry_quarter
    out = _from_other_arg(criteria=pr_qrt, output=tas_qrt, op=xr_op, freq=freq)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indices\\_anuclim.py", line 521, in _from_other_arg
    return ds.resample(time=freq).map(get_other_op)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\resample.py", line 304, in map
    combined = self._combine(applied)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\groupby.py", line 954, in _combine
    applied_example, applied = peek_at(applied)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\utils.py", line 182, in peek_at
    peek = next(gen)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\resample.py", line 303, in <genexpr>
    applied = (func(ds, *args, **kwargs) for ds in self._iter_grouped())
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indices\\_anuclim.py", line 519, in get_other_op
    return lazy_indexing(ds.output, index=op(ds.criteria, dim=dim), dim=dim)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indices\\run_length.py", line 745, in lazy_indexing
    dask="allowed",
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\computation.py", line 1064, in apply_ufunc
    keep_attrs=keep_attrs,
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\computation.py", line 233, in apply_dataarray_vfunc
    result_var = func(*data_vars)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xarray\\core\\computation.py", line 604, in apply_variable_ufunc
    result_data = func(*input_data)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\xclim\\indices\\run_length.py", line 737, in _isel
    return np.take_along_axis(array, indx[..., None], -1)[..., 0]
  File "<__array_function__ internals>", line 6, in take_along_axis
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\array\\core.py", line 1354, in __array_function__
    return handle_nonmatching_names(func, args, kwargs)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\array\\core.py", line 1339, in handle_nonmatching_names
    args, kwargs = compute(args, kwargs)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\base.py", line 436, in compute
    results = schedule(dsk, keys, **kwargs)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\threaded.py", line 81, in get
    **kwargs
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\local.py", line 486, in get_async
    raise_exception(exc, tb)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\local.py", line 316, in reraise
    raise exc
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\local.py", line 222, in execute_task
    result = _execute_task(task, data)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\core.py", line 118, in _execute_task
    args2 = [_execute_task(a, cache) for a in args]
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\core.py", line 118, in <listcomp>
    args2 = [_execute_task(a, cache) for a in args]
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\core.py", line 119, in _execute_task
    return func(*args2)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\toolz\\functoolz.py", line 488, in __call__
    ret = f(ret)
  File "C:\\Anaconda\\envs\\xclim\\lib\\site-packages\\dask\\array\\reductions.py", line 920, in nanarg_agg
    raise ValueError("All NaN slice encountered")
ValueError: All NaN slice encountered

```


Changed files:
tests/test_utils.py, xclim/indices/run_length.py, xclim/indices/_anuclim.py, tests/test_run_length.py, xclim/sdba/utils.py, xclim/core/calendar.py, xclim/core/utils.py

diff --git a/tests/test_run_length.py b/tests/test_run_length.py
index 391cb08e..90e2a0ae 100644
--- a/tests/test_run_length.py
+++ b/tests/test_run_length.py
@@ -350,3 +350,21 @@ class TestRunsWithDates:
 
         out = runs.resample(time="MS").map(func, window=1, date="07-01", dim="time")
         assert out.isnull().all()
+
+
+@pytest.mark.parametrize("use_dask", [True, False])
+def test_lazy_indexing_special_cases(use_dask):
+    a = xr.DataArray(np.random.rand(10, 10, 10), dims=("x", "y", "z"))
+    b = xr.DataArray(np.random.rand(10, 10, 10), dims=("x", "y", "z"))
+
+    if use_dask:
+        a = a.chunk({"y": 5, "z": 5})
+        b = b.chunk({"y": 5, "z": 1})
+
+    with pytest.raises(ValueError):
+        rl.lazy_indexing(a, b)
+
+    b = b.argmin("x").argmin("y")
+
+    with pytest.raises(ValueError, match="more than one dimension more than index"):
+        rl.lazy_indexing(a, b)
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 828af424..2c7184d0 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -4,8 +4,10 @@
 from inspect import signature
 
 import numpy as np
+import xarray as xr
 
 from xclim.core.indicator import Indicator
+from xclim.core.utils import ensure_chunk_size
 from xclim.core.utils import walk_map
 from xclim.core.utils import wrapped_partial
 
@@ -60,3 +62,19 @@ def test_wrapped_indicator(tas_series):
 
     assert ind2(tas, tas2) == 366
     assert ind1(tas, thresh=1111) == 366
+
+
+def test_ensure_chunk_size():
+    da = xr.DataArray(np.zeros((20, 21, 20)), dims=("x", "y", "z"))
+
+    out = ensure_chunk_size(da, x=10, y=-1)
+
+    assert da is out
+
+    dac = da.chunk({"x": (1,) * 20, "y": (10, 10, 1), "z": (10, 10)})
+
+    out = ensure_chunk_size(dac, x=3, y=5, z=-1)
+
+    assert out.chunks[0] == (3, 3, 3, 3, 3, 5)
+    assert out.chunks[1] == (10, 11)
+    assert out.chunks[2] == (20,)
diff --git a/xclim/core/calendar.py b/xclim/core/calendar.py
index ffe6e1c8..41a9ee27 100644
--- a/xclim/core/calendar.py
+++ b/xclim/core/calendar.py
@@ -52,13 +52,15 @@ max_doy = {
 }
 
 
-def get_calendar(arr: Union[xr.DataArray, xr.Dataset]) -> str:
+def get_calendar(arr: Union[xr.DataArray, xr.Dataset], dim: str = "time") -> str:
     """Return the calendar of the time coord of the DataArray
 
     Parameters
     ----------
-    arr : xr.DataArray
-      Array with `time` coordinate. Values must either be of datetime64 dtype or have a dt.calendar attribute.
+    arr : Union[xr.DataArray, xr.Dataset]
+      Array/dataset with a datetime coordinate. Values must either be of datetime64 dtype or have a dt.calendar attribute.
+    dim : str
+      Name of the coordinate to check.
 
     Raises
     ------
@@ -70,14 +72,14 @@ def get_calendar(arr: Union[xr.DataArray, xr.Dataset]) -> str:
     str
       The cftime calendar name or "default" when the data is using numpy's datetime type (numpy.datetime64.
     """
-    if arr.time.dtype == "O":  # Assume cftime, if it fails, not our fault
-        non_na_item = arr.time.where(arr.time.notnull(), drop=True)[0].item()
+    if arr[dim].dtype == "O":  # Assume cftime, if it fails, not our fault
+        non_na_item = arr[dim].where(arr[dim].notnull(), drop=True)[0].item()
         cal = non_na_item.calendar
-    elif "datetime64" in arr.time.dtype.name:
+    elif "datetime64" in arr[dim].dtype.name:
         cal = "default"
     else:
         raise ValueError(
-            f"Cannot infer calendars from timeseries of type {arr.time[0].dtype}"
+            f"Cannot infer calendars from timeseries of type {arr[dim][0].dtype}"
         )
     return cal
 
@@ -86,6 +88,7 @@ def convert_calendar(
     source: Union[xr.DataArray, xr.Dataset],
     target: Union[xr.DataArray, str],
     align_on: Optional[str] = None,
+    dim: str = "time",
 ) -> xr.DataArray:
     """Convert a DataArray/Dataset to another calendar using the specified method.
     Only converts the individual timestamps, does not modify any data except in dropping invalid/surplus dates.
@@ -108,6 +111,8 @@ def convert_calendar(
       that are missing in the converted `source` are filled by NaNs.
     align_on : {None, 'date', 'year'}
       Must be specified when either source or target is a `360_day` calendar, ignored otherwise. See Notes.
+    dim : str
+      Name of the time coordinate.
 
     Returns
     -------
@@ -147,12 +152,12 @@ def convert_calendar(
 
         This option is best used with data on a frequency coarser than daily.
     """
-    cal_src = get_calendar(source)
+    cal_src = get_calendar(source, dim=dim)
 
     if isinstance(target, str):
         cal_tgt = target
     else:
-        cal_tgt = get_calendar(target)
+        cal_tgt = get_calendar(target, dim=dim)
 
     if cal_src == cal_tgt:
         return source
@@ -182,36 +187,36 @@ def convert_calendar(
                 / days_in_year(yr, cal_src)
             ).astype(int)
 
-        new_doy = source.time.groupby("time.year").map(_yearly_interp_doy)
+        new_doy = source.time.groupby(f"{dim}.year").map(_yearly_interp_doy)
 
         # Convert the source datetimes, but override the doy with our new doys
-        out["time"] = xr.DataArray(
+        out[dim] = xr.DataArray(
             [
                 _convert_datetime(datetime, new_doy=doy, calendar=cal_tgt)
-                for datetime, doy in zip(source.time.indexes["time"], new_doy)
+                for datetime, doy in zip(source[dim].indexes[dim], new_doy)
             ],
-            dims=("time",),
-            name="time",
+            dims=(dim,),
+            name=dim,
         )
         # Remove duplicate timestamps, happens when reducing the number of days
-        out = out.isel(time=np.unique(out.time, return_index=True)[1])
+        out = out.isel({dim: np.unique(out[dim], return_index=True)[1]})
     else:
-        time_idx = source.time.indexes["time"]
-        out["time"] = xr.DataArray(
+        time_idx = source[dim].indexes[dim]
+        out[dim] = xr.DataArray(
             [_convert_datetime(time, calendar=cal_tgt) for time in time_idx],
-            dims=("time",),
-            name="time",
+            dims=(dim,),
+            name=dim,
         )
         # Remove NaN that where put on invalid dates in target calendar
-        out = out.where(out.time.notnull(), drop=True)
+        out = out.where(out[dim].notnull(), drop=True)
 
     if isinstance(target, xr.DataArray):
-        out = out.reindex(time=target)
+        out = out.reindex({dim: target})
     return out
 
 
 def interp_calendar(
-    source: Union[xr.DataArray, xr.Dataset], target: xr.DataArray,
+    source: Union[xr.DataArray, xr.Dataset], target: xr.DataArray, dim: str = "time",
 ) -> xr.DataArray:
     """Interpolates a DataArray/Dataset to another calendar based on decimal year measure.
 
@@ -228,22 +233,22 @@ def interp_calendar(
       The source data to interpolate, must have a time coordinate of a valid dtype (np.datetime64 or cftime objects)
     target: xr.DataArray
       The target time coordinate of a valid dtype (np.datetime64 or cftime objects)
+    dim : str
+      The time coordinate name.
 
     Return
     ------
     Union[xr.DataArray, xr.Dataset]
       The source interpolated on the decimal years of target,
     """
-    cal_src = get_calendar(source)
-    cal_tgt = get_calendar(target)
+    cal_src = get_calendar(source, dim=dim)
+    cal_tgt = get_calendar(target, dim=dim)
 
     out = source.copy()
-    out["time"] = datetime_to_decimal_year(source.time, calendar=cal_src).drop_vars(
-        "time"
-    )
+    out[dim] = datetime_to_decimal_year(source[dim], calendar=cal_src).drop_vars(dim)
     target_idx = datetime_to_decimal_year(target, calendar=cal_tgt)
     out = out.interp(time=target_idx)
-    out["time"] = target
+    out[dim] = target
     return out
 
 
diff --git a/xclim/core/utils.py b/xclim/core/utils.py
index fefa3b90..4189bc27 100644
--- a/xclim/core/utils.py
+++ b/xclim/core/utils.py
@@ -10,6 +10,9 @@ from collections import defaultdict
 from functools import partial
 from types import FunctionType
 
+import dask.array as dsk
+import numpy as np
+import xarray as xr
 from boltons.funcutils import update_wrapper
 
 
@@ -82,3 +85,51 @@ class ValidationError(ValueError):
     @property
     def msg(self):
         return self.args[0]
+
+
+def ensure_chunk_size(da: xr.DataArray, max_iter=10, **minchunks):
+    """Ensure that the input dataarray has chunks of at least the given size.
+
+    If only one chunk is too small, it is merged with an adjacent chunk.
+    If many chunks are too small, they are grouped together by merging adjacent chunks.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+      The input dataarray, with or without the dask backend. Does nothing when passed a non-dask array.
+    **minchunks : int
+      A kwarg mapping from dimension name to minimum chunk size.
+      Pass -1 to force a single chunk along that dimension.
+    """
+    if not isinstance(da.data, dsk.Array):
+        return da
+
+    all_chunks = dict(zip(da.dims, da.chunks))
+    chunking = {}
+    for dim, minchunk in minchunks.items():
+        chunks = all_chunks[dim]
+        if minchunk == -1 and len(chunks) > 1:
+            # Rechunk to single chunk only if it's not already one
+            chunking[dim] = -1
+
+        toosmall = np.array(chunks) < minchunk  # Chunks that are too small
+        if toosmall.sum() > 1:
+            # Many chunks are too small, merge them by groups
+            fac = np.ceil(minchunk / min(chunks)).astype(int)
+            chunking[dim] = tuple(
+                sum(chunks[i : i + fac]) for i in range(0, len(chunks), fac)
+            )
+            # Reset counter is case the last chunks are still too small
+            chunks = chunking[dim]
+            toosmall = np.array(chunks) < minchunk
+        if toosmall.sum() == 1:
+            # Only one, merge it with adjacent chunk
+            ind = np.where(toosmall)[0][0]
+            new_chunks = list(chunks)
+            sml = new_chunks.pop(ind)
+            new_chunks[max(ind - 1, 0)] += sml
+            chunking[dim] = tuple(new_chunks)
+
+    if chunking:
+        return da.chunk(chunks=chunking)
+    return da
diff --git a/xclim/indices/_anuclim.py b/xclim/indices/_anuclim.py
index eaf302b2..adb809f6 100644
--- a/xclim/indices/_anuclim.py
+++ b/xclim/indices/_anuclim.py
@@ -1,3 +1,4 @@
+import numpy as np
 import xarray
 
 from ._multivariate import daily_temperature_range
@@ -11,6 +12,7 @@ from xclim.core.units import declare_units
 from xclim.core.units import pint_multiply
 from xclim.core.units import units
 from xclim.core.units import units2pint
+from xclim.core.utils import ensure_chunk_size
 
 
 xarray.set_options(enable_cftimeindex=True)  # Set xarray to use cftimeindex
@@ -516,7 +518,9 @@ def _from_other_arg(criteria, output, op, freq):
     dim = "time"
 
     def get_other_op(ds):
-        return lazy_indexing(ds.output, index=op(ds.criteria, dim=dim), dim=dim)
+        all_nans = ds.criteria.isnull().all(dim=dim)
+        index = op(ds.criteria.where(~all_nans, 0), dim=dim)
+        return lazy_indexing(ds.output, index=index, dim=dim).where(~all_nans)
 
     return ds.resample(time=freq).map(get_other_op)
 
@@ -547,6 +551,11 @@ def _to_quarter(freq, pr=None, tas=None):
             f'Unknown input time frequency "{freq}": must be one of "daily", "weekly" or "monthly".'
         )
 
+    if tas is not None:
+        tas = ensure_chunk_size(tas, time=np.ceil(window / 2))
+    if pr is not None:
+        pr = ensure_chunk_size(pr, time=np.ceil(window / 2))
+
     with xarray.set_options(keep_attrs=True):
         if pr is not None:
             pr = pint_multiply(pr, 1 * u, "mm")
@@ -560,4 +569,5 @@ def _to_quarter(freq, pr=None, tas=None):
             )
             out.attrs = tas.attrs
 
+    out = ensure_chunk_size(out, time=-1)
     return out
diff --git a/xclim/indices/run_length.py b/xclim/indices/run_length.py
index 85dfc193..95ba91c8 100644
--- a/xclim/indices/run_length.py
+++ b/xclim/indices/run_length.py
@@ -18,7 +18,6 @@ import dask.array as dsk
 import numpy as np
 import xarray as xr
 
-
 logging.captureWarnings(True)
 npts_opt = 9000
 
@@ -696,7 +695,7 @@ def first_run_ufunc(x: xr.DataArray, window: int, dim: str = "time",) -> xr.appl
     return ind
 
 
-def lazy_indexing(da: xr.DataArray, index: xr.DataArray, dim=None, drop=True):
+def lazy_indexing(da: xr.DataArray, index: xr.DataArray, dim=None):
     """Get values of `da` at indices `index` in a NaN-aware and lazy manner.
 
     The algorithm differs whether da is 1D or not.
@@ -704,9 +703,9 @@ def lazy_indexing(da: xr.DataArray, index: xr.DataArray, dim=None, drop=True):
     Parameters
     ----------
     da : xr.DataArray
-      Input array. If not 1D, `dim` must be given and the corresponding dimension shouldn't be chunked.
+      Input array. If not 1D, `dim` must be given and must not appear in index.
     index : xr.DataArray
-      N-d integer indices
+      N-d integer indices, all dimensions of index must be in da
     dim : Dimension along which to index,
           unused if `da` is 1D, should not be present in `index`.
 
@@ -732,15 +731,27 @@ def lazy_indexing(da: xr.DataArray, index: xr.DataArray, dim=None, drop=True):
             out = out.drop_vars(da.dims[0])
         return out
 
-    # else:
-    def _isel(array, indx):
-        return np.take_along_axis(array, indx[..., None], -1)[..., 0]
+    if dim is None:
+        diff_dims = set(da.dims) - set(index.dims)
+        if len(diff_dims) == 0:
+            raise ValueError(
+                "da must have at least one dimension more than index for lazy_indexing."
+            )
+        if len(diff_dims) > 1:
+            raise ValueError(
+                "If da has more than one dimension more than index, the indexing dim must be given through `dim`"
+            )
+        dim = diff_dims.pop()
+
+    def _index_from_nd_array(array, indices):
+        return np.take_along_axis(array, indices[..., np.newaxis], axis=-1)[..., 0]
 
     return xr.apply_ufunc(
-        _isel,
+        _index_from_nd_array,
         da,
         index,
         input_core_dims=[[dim], []],
         output_core_dims=[[]],
-        dask="allowed",
+        dask="parallelized",
+        output_dtypes=[da.dtype],
     )
diff --git a/xclim/sdba/utils.py b/xclim/sdba/utils.py
index 5bdc3ab7..24e51e1e 100644
--- a/xclim/sdba/utils.py
+++ b/xclim/sdba/utils.py
@@ -14,7 +14,7 @@ from scipy.interpolate import interp1d
 from .base import Grouper
 from .base import parse_group
 from xclim.core.calendar import _interpolate_doy_calendar
-
+from xclim.core.utils import ensure_chunk_size
 
 MULTIPLICATIVE = "*"
 ADDITIVE = "+"
@@ -278,10 +278,7 @@ def add_cyclic_bounds(da: xr.DataArray, att: str, cyclic_coords: bool = True):
         vals[-1] = vals[-2] + diff[-1]
         qmf = qmf.assign_coords({att: vals})
         qmf[att].attrs.update(da.coords[att].attrs)
-    if da.chunks is not None and len(da.chunks[da.get_axis_num(att)]) == 1:
-        # Downstream methos do not expect this method to change the number of chunks
-        return qmf.chunk({att: -1})
-    return qmf
+    return ensure_chunk_size(qmf, **{att: -1})
 
 
 def extrapolate_qm(qf: xr.DataArray, xq: xr.DataArray, method: str = "constant"):
@@ -365,10 +362,7 @@ def add_endpoints(
         elems.append(y)
     l, r = elems  # pylint: disable=unbalanced-tuple-unpacking
     out = xr.concat((l, da, r), dim=dim)
-    if da.chunks is not None and len(da.chunks[da.get_axis_num(dim)]) == 1:
-        # Downstream methos do not expect this method to change the number of chunks
-        return out.chunk({dim: -1})
-    return out
+    return ensure_chunk_size(out, **{dim: -1})
 
 
 @parse_group