From 9811ebe0ba4ce642f68c8048f8cf01ac0aaf8990 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 12:28:33 +0200 Subject: [PATCH 01/11] Branch off lazy mmstats PR --- esmvalcore/_recipe_checks.py | 5 +- esmvalcore/preprocessor/_multimodel.py | 526 +++++++++--------- .../multimodel_statistics/test_multimodel.py | 73 ++- .../_multimodel/test_multimodel.py | 231 ++++++-- 4 files changed, 512 insertions(+), 323 deletions(-) diff --git a/esmvalcore/_recipe_checks.py b/esmvalcore/_recipe_checks.py index 21cff6a4a9..1133a36639 100644 --- a/esmvalcore/_recipe_checks.py +++ b/esmvalcore/_recipe_checks.py @@ -9,7 +9,8 @@ import yamale from ._data_finder import get_start_end_year -from .preprocessor import PreprocessingTask, TIME_PREPROCESSORS +from .preprocessor import TIME_PREPROCESSORS, PreprocessingTask +from .preprocessor._multimodel import STATISTIC_MAPPING logger = logging.getLogger(__name__) @@ -191,7 +192,7 @@ def extract_shape(settings): def valid_multimodel_statistic(statistic): """Check that `statistic` is a valid argument for multimodel stats.""" - valid_names = ["mean", "median", "std", "min", "max"] + valid_names = ['std'] + list(STATISTIC_MAPPING.keys()) valid_patterns = [r"^(p\d{1,2})(\.\d*)?$"] if not (statistic in valid_names or re.match(r'|'.join(valid_patterns), statistic)): diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index ad010ae097..483019c6ec 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -3,184 +3,65 @@ import logging import re from datetime import datetime -from functools import partial, reduce +from functools import reduce import cf_units +import dask.array as da import iris import numpy as np -import scipy +from iris.util import equalise_attributes logger = logging.getLogger(__name__) +STATISTIC_MAPPING = { + 'gmean': iris.analysis.GMEAN, # not lazy in iris + 'hmean': iris.analysis.HMEAN, # not lazy in iris + 'max': iris.analysis.MAX, + 'median': iris.analysis.MEDIAN, # not lazy in iris + 'min': iris.analysis.MIN, + 'rms': iris.analysis.RMS, + 'sum': iris.analysis.SUM, + 'mean': iris.analysis.MEAN, + 'std_dev': iris.analysis.STD_DEV, + 'variance': iris.analysis.VARIANCE, + # The following require extra kwargs, + # atm this is only supported for percentiles via e.g. `pXX` + 'count': iris.analysis.COUNT, + 'peak': iris.analysis.PEAK, + 'percentile': iris.analysis.PERCENTILE, # not lazy in iris + 'proportion': iris.analysis.PROPORTION, # not lazy in iris + 'wpercentile': iris.analysis.WPERCENTILE, # not lazy in iris +} + +CONCAT_DIM = 'multi-model' + + +def _resolve_operator(statistic: str): + """Find the operator corresponding to the statistic.""" + statistic = statistic.lower() + kwargs = {} + + # special cases + if statistic == 'std': + logger.warning( + "Multicube statistics is aligning its behaviour with iris.analysis" + ". Please consider replacing 'std' with 'std_dev' in your code.") + statistic = 'std_dev' + + elif re.match(r"^(p\d{1,2})(\.\d*)?$", statistic): + # percentiles between p0 and p99.99999... + percentile = float(statistic[1:]) + kwargs['percent'] = percentile + statistic = 'percentile' -def _plev_fix(dataset, pl_idx): - """Extract valid plev data. - - This function takes care of situations in which certain plevs are - completely masked due to unavailable interpolation boundaries. - """ - if np.ma.is_masked(dataset): - # keep only the valid plevs - if not np.all(dataset.mask[pl_idx]): - statj = np.ma.array(dataset[pl_idx], mask=dataset.mask[pl_idx]) - else: - logger.debug('All vals in plev are masked, ignoring.') - statj = None - else: - mask = np.zeros_like(dataset[pl_idx], bool) - statj = np.ma.array(dataset[pl_idx], mask=mask) - - return statj + try: + operator = STATISTIC_MAPPING[statistic] + except KeyError as err: + raise ValueError( + f'Statistic `{statistic}` not supported by multicube statistics. ' + f'Must be one of {tuple(STATISTIC_MAPPING.keys())}.') from err - -def _quantile(data, axis, quantile): - """Calculate quantile. - - Workaround for calling scipy's mquantiles with arrays of >2 - dimensions Similar to iris' _percentiles function, see their - discussion: https://github.com/SciTools/iris/pull/625 - """ - # Ensure that the target axis is the last dimension. - data = np.rollaxis(data, axis, start=data.ndim) - shape = data.shape[:-1] - # Flatten any leading dimensions. - if shape: - data = data.reshape([np.prod(shape), data.shape[-1]]) - # Perform the quantile calculation. - result = scipy.stats.mstats.mquantiles(data, - quantile, - axis=-1, - alphap=1, - betap=1) - # Ensure to unflatten any leading dimensions. - if shape: - result = result.reshape(shape) - # Check whether to reduce to a scalar result - if result.shape == (1, ): - result = result[0] - - return result - - -def _compute_statistic(data, statistic_name): - """Compute multimodel statistic.""" - data = np.ma.array(data) - statistic = data[0] - - if statistic_name == 'median': - statistic_function = np.ma.median - elif statistic_name == 'mean': - statistic_function = np.ma.mean - elif statistic_name == 'std': - statistic_function = np.ma.std - elif statistic_name == 'max': - statistic_function = np.ma.max - elif statistic_name == 'min': - statistic_function = np.ma.min - elif re.match(r"^(p\d{1,2})(\.\d*)?$", statistic_name): - # percentiles between p0 and p99.99999... - quantile = float(statistic_name[1:]) / 100 - statistic_function = partial(_quantile, quantile=quantile) - else: - raise ValueError(f'No such statistic: `{statistic_name}`') - - # no plevs - if len(data[0].shape) < 3: - # get all NOT fully masked data - u_data - # data is per time point - # so we can safely NOT compute stats for single points - if data.ndim == 1: - u_datas = data - else: - u_datas = [d for d in data if not np.all(d.mask)] - if len(u_datas) > 1: - statistic = statistic_function(data, axis=0) - else: - statistic.mask = True - return statistic - - # plevs - for j in range(statistic.shape[0]): - plev_check = [] - for cdata in data: - fixed_data = _plev_fix(cdata, j) - if fixed_data is not None: - plev_check.append(fixed_data) - - # check for nr datasets - if len(plev_check) > 1: - plev_check = np.ma.array(plev_check) - statistic[j] = statistic_function(plev_check, axis=0) - else: - statistic.mask[j] = True - - return statistic - - -def _put_in_cube(template_cube, cube_data, statistic, t_axis): - """Quick cube building and saving.""" - tunits = template_cube.coord('time').units - times = iris.coords.DimCoord(t_axis, - standard_name='time', - units=tunits, - var_name='time', - long_name='time') - times.bounds = None - times.guess_bounds() - - coord_names = [c.long_name for c in template_cube.coords()] - coord_names.extend([c.standard_name for c in template_cube.coords()]) - if 'latitude' in coord_names: - lats = template_cube.coord('latitude') - else: - lats = None - if 'longitude' in coord_names: - lons = template_cube.coord('longitude') - else: - lons = None - - # no plevs - if len(template_cube.shape) == 3: - cspec = [(times, 0), (lats, 1), (lons, 2)] - # plevs - elif len(template_cube.shape) == 4: - if 'air_pressure' in coord_names: - plev = template_cube.coord('air_pressure') - elif 'altitude' in coord_names: - plev = template_cube.coord('altitude') - else: - raise ValueError(f"4D cube with dimensions {coord_names} not " - f"supported at the moment") - cspec = [(times, 0), (plev, 1), (lats, 2), (lons, 3)] - elif len(template_cube.shape) == 1: - cspec = [ - (times, 0), - ] - elif len(template_cube.shape) == 2: - # If you're going to hardwire air_pressure into this, - # might as well have depth here too. - plev = template_cube.coord('depth') - cspec = [ - (times, 0), - (plev, 1), - ] - - # correct dspec if necessary - fixed_dspec = np.ma.fix_invalid(cube_data, copy=False, fill_value=1e+20) - # put in cube - stats_cube = iris.cube.Cube(fixed_dspec, - dim_coords_and_dims=cspec, - long_name=statistic) - coord_names = [coord.name() for coord in template_cube.coords()] - if 'air_pressure' in coord_names: - if len(template_cube.shape) == 3: - stats_cube.add_aux_coord(template_cube.coord('air_pressure')) - - stats_cube.var_name = template_cube.var_name - stats_cube.long_name = template_cube.long_name - stats_cube.standard_name = template_cube.standard_name - stats_cube.units = template_cube.units - return stats_cube + return operator, kwargs def _get_consistent_time_unit(cubes): @@ -218,7 +99,6 @@ def _unify_time_coordinates(cubes): if 0 not in np.diff(years): # yearly data dates = [datetime(year, 7, 1, 0, 0, 0) for year in years] - elif 0 not in np.diff(months): # monthly data dates = [ @@ -248,115 +128,219 @@ def _unify_time_coordinates(cubes): cube.coord('time').guess_bounds() -def _get_time_slice(cubes, time): - """Fill time slice array with cubes' data if time in cube, else mask.""" - time_slice = [] - for cube in cubes: - cube_time = cube.coord('time').points - if time in cube_time: - idx = int(np.argwhere(cube_time == time)) - subset = cube.data[idx] - else: - subset = np.ma.empty(list(cube.shape[1:])) - subset.mask = True - time_slice.append(subset) - return time_slice +def _time_coords_are_aligned(cubes): + """Return `True` if time coords are aligned.""" + first_time_array = cubes[0].coord('time').points + for cube in cubes[1:]: + other_time_array = cube.coord('time').points + if not np.array_equal(first_time_array, other_time_array): + return False -def _assemble_data(cubes, statistic, span='overlap'): - """Get statistical data in iris cubes.""" - # New time array representing the union or intersection of all cubes - time_spans = [cube.coord('time').points for cube in cubes] - if span == 'overlap': - new_times = reduce(np.intersect1d, time_spans) - elif span == 'full': - new_times = reduce(np.union1d, time_spans) - n_times = len(new_times) + return True - # Target array to populate with computed statistics - new_shape = [n_times] + list(cubes[0].shape[1:]) - stats_data = np.ma.zeros(new_shape, dtype=np.dtype('float32')) - # Realize all cubes at once instead of separately for each time slice - _ = [cube.data for cube in cubes] +def _subset(cube, time_points): + """Subset cube to given time range.""" + begin = cube.coord('time').units.num2date(time_points[0]) + end = cube.coord('time').units.num2date(time_points[-1]) + constraint = iris.Constraint(time=lambda cell: begin <= cell.point <= end) + return cube.extract(constraint) - # Make time slices and compute stats - for i, time in enumerate(new_times): - time_data = _get_time_slice(cubes, time) - stats_data[i] = _compute_statistic(time_data, statistic) - template = cubes[0] - stats_cube = _put_in_cube(template, stats_data, statistic, new_times) - return stats_cube +def _extend(cube, time_points): + """Extend cube to a specified time range. + If time points are missing before the start/after the end of the + time range, cubes for each missing time pointwith masked data will + be added to pad the time range to match `time_points`. This method + supports lazy operation. + """ + cube.coord('time').bounds = None + cube_points = cube.coord('time').points -def _multicube_statistics(cubes, statistics, span): - """Compute statistics over multiple cubes. + begin = cube_points[0] + end = cube_points[-1] - Can be used e.g. for ensemble or multi-model statistics. + pad_begin = time_points[time_points < begin] + pad_end = time_points[time_points > end] - This function was designed to work on (max) four-dimensional data: - time, vertical axis, two horizontal axes. + if (len(pad_begin)) == 0 and (len(pad_end) == 0): + return cube - Apart from the time coordinate, cubes must have consistent shapes. There - are two options to combine time coordinates of different lengths, see - the `span` argument. + template_cube = cube[:1].copy() + template_cube.data = np.ma.array(da.zeros_like(template_cube.data), + mask=True, + dtype=template_cube.data.dtype) - Parameters - ---------- - cubes: list - list of cubes over which the statistics will be computed; - statistics: list - statistical metrics to be computed. Available options: mean, median, - max, min, std, or pXX.YY (for percentile XX.YY; decimal part optional). - span: str - overlap or full; if overlap, statitsticss are computed on common time- - span; if full, statistics are computed on full time spans, ignoring - missing data. + cube_list = [] - Returns - ------- - dict - dictionary of statistics cubes with statistics' names as keys. + for time_point in pad_begin: + new_slice = template_cube.copy() + new_slice.coord('time').points = float(time_point) + cube_list.append(new_slice) - Raises - ------ - ValueError - If span is neither overlap nor full. - """ - if len(cubes) < 2: - logger.info('Found only 1 cube; no statistics computed for %r', - list(cubes)[0]) - return {statistic: cubes[0] for statistic in statistics} + cube_list.append(cube) + + for time_point in pad_end: + new_slice = template_cube.copy() + new_slice.coord('time').points = float(time_point) + cube_list.append(new_slice) + + cube_list = iris.cube.CubeList(cube_list) + + new_cube = cube_list.concatenate_cube() - logger.debug('Multicube statistics: computing: %s', statistics) + return new_cube - # Reset time coordinates and make cubes share the same calendar + +def _align(cubes, span): + """Expand or subset cubes so they share a common time span.""" _unify_time_coordinates(cubes) - # Check whether input is valid + if _time_coords_are_aligned(cubes): + return cubes + + all_time_arrays = [cube.coord('time').points for cube in cubes] + if span == 'overlap': - # check if we have any time overlap - times = [cube.coord('time').points for cube in cubes] - overlap = reduce(np.intersect1d, times) - if len(overlap) <= 1: - logger.info("Time overlap between cubes is none or a single point." - "check datasets: will not compute statistics.") - return cubes - logger.debug("Using common time overlap between " - "datasets to compute statistics.") + common_time_points = reduce(np.intersect1d, all_time_arrays) + new_cubes = [_subset(cube, common_time_points) for cube in cubes] elif span == 'full': - logger.debug("Using full time spans to compute statistics.") + all_time_points = reduce(np.union1d, all_time_arrays) + new_cubes = [_extend(cube, all_time_points) for cube in cubes] else: - raise ValueError( - "Unexpected value for span {}, choose from 'overlap', 'full'". - format(span)) + raise ValueError(f"Invalid argument for span: {span!r}" + "Must be one of 'overlap', 'full'.") + + for cube in new_cubes: + # Make sure bounds exist and are consistent + cube.coord('time').bounds = None + cube.coord('time').guess_bounds() + + return new_cubes + + +def _combine(cubes): + """Merge iris cubes into a single big cube with new dimension. + + This assumes that all input cubes have the same shape. + """ + equalise_attributes(cubes) # in-place + + for i, cube in enumerate(cubes): + concat_dim = iris.coords.AuxCoord(i, var_name=CONCAT_DIM) + + cube.add_aux_coord(concat_dim) + + # Clear some metadata that can cause merge to fail + # https://scitools-iris.readthedocs.io/en/stable/userguide/ + # merge_and_concat.html#common-issues-with-merge-and-concatenate + + cube.cell_methods = None + + for coord in cube.coords(): + coord.long_name = None + coord.attributes = None + + cubes = iris.cube.CubeList(cubes) + + merged_cube = cubes.merge_cube() + + return merged_cube + + +def rechunk(cube): + """Rechunk the cube to speed up out-of-memory computation.""" + new_chunks = {0: -1} # don't chunk along the multimodel dimension + if cube.ndim > 1: + new_chunks[1] = 'auto' # do chunk along the first subsequent dimension + + cube.data = cube.lazy_data().rechunk(new_chunks) + + logger.debug("Total data size: %s MB", cube.lazy_data().nbytes * 1e-6) + logger.debug("New chunk block size: %s MB", + cube.lazy_data().nbytes / cube.lazy_data().npartitions * 1e-6) + logger.debug("New chunk configuration: %s", cube.lazy_data()) + + +def _compute_eager(cubes: list, *, operator: iris.analysis.Aggregator, + **kwargs): + """Loop over slices of a cube if iris has no lazy aggregator.""" + _ = [cube.data for cube in cubes] # make sure the cubes' data are realized + + result_slices = [] + for i in range(cubes[0].shape[0]): + single_model_slices = [cube[i] for cube in cubes + ] # maybe filter the iris warning here? + combined_slice = _combine(single_model_slices) + collapsed_slice = combined_slice.collapsed(CONCAT_DIM, operator, + **kwargs) + result_slices.append(collapsed_slice) + + result_cube = iris.cube.CubeList(result_slices).merge_cube() + + # For consistency with lazy procedure + result_cube.data = np.ma.array(result_cube.data) + + result_cube.remove_coord(CONCAT_DIM) + + return result_cube + + +def _compute_lazy(cubes: list, *, operator: iris.analysis.Aggregator, + **kwargs): + """Compute statistics using lazy iris function.""" + cube = _combine( + cubes) # this is now done for each statistic, can we avoid that? + rechunk(cube) + + # This will always return a masked array + result_cube = cube.collapsed(CONCAT_DIM, operator, **kwargs) + result_cube.remove_coord(CONCAT_DIM) + + for cube in cubes: + cube.remove_coord(CONCAT_DIM) + + return result_cube + + +def _multicube_statistics(cubes, statistics, span): + """Compute statistics over multiple cubes. + + Can be used e.g. for ensemble or multi-model statistics. + + Cubes are merged and subsequently collapsed along a new auxiliary + coordinate. Inconsistent attributes will be removed. + """ + if len(cubes) == 1: + raise ValueError('Cannot perform multicube statistics ' + 'for a single cube.') + + lazy_input = bool(all(cube.has_lazy_data() for cube in cubes)) + + copied_cubes = [cube.copy() for cube in cubes] # avoid modifying inputs + aligned_cubes = _align(copied_cubes, span=span) - # Compute statistics statistics_cubes = {} for statistic in statistics: - statistic_cube = _assemble_data(cubes, statistic, span) - statistics_cubes[statistic] = statistic_cube + logger.debug('Multicube statistics: computing: %s', statistic) + operator, kwargs = _resolve_operator(statistic) + + if operator.lazy_func is None: + result_cube = _compute_eager(aligned_cubes, + operator=operator, + **kwargs) + else: + result_cube = _compute_lazy(aligned_cubes, + operator=operator, + **kwargs) + + # lazy input --> lazy output + result_cube.data = result_cube.lazy_data( + ) if lazy_input else result_cube.data + + statistics_cubes[statistic] = result_cube return statistics_cubes @@ -399,23 +383,39 @@ def multi_model_statistics(products, keep_input_datasets=True): """Compute multi-model statistics. - This function computes multi-model statistics on cubes or products. - Products (or: preprocessorfiles) are used internally by ESMValCore to store + This function computes multi-model statistics on a list of ``products``, + which can be instances of :py:class:`~iris.cube.Cube` or + :py:class:`~esmvalcore.preprocessor.PreprocessorFile`. + The latter is used internally by ESMValCore to store workflow and provenance information, and this option should typically be ignored. - This function was designed to work on (max) four-dimensional data: time, - vertical axis, two horizontal axes. Apart from the time coordinate, cubes - must have consistent shapes. There are two options to combine time - coordinates of different lengths, see the `span` argument. + Apart from the time coordinate, cubes must have consistent shapes. There + are two options to combine time coordinates of different lengths, see the + ``span`` argument. + + Uses the statistical operators in :py:mod:`iris.analysis`, including + ``mean``, ``median``, ``min``, ``max``, and ``std``. Percentiles are also + supported and can be specified like ``pXX.YY`` (for percentile ``XX.YY``; + decimal part optional). + + Notes + ----- + Some of the operators in :py:mod:`iris.analysis` require additional + arguments. Except for percentiles, these operators are currently not + supported. + + Lazy operation is supported for all statistics, except + ``median``, ``percentile``, ``gmean`` and ``hmean``. Parameters ---------- products: list Cubes (or products) over which the statistics will be computed. statistics: list - Statistical metrics to be computed. Available options: mean, median, - max, min, std, or pXX.YY (for percentile XX.YY; decimal part optional). + Statistical metrics to be computed, e.g. [``mean``, ``max``]. Choose + from the operators listed in the iris.analysis package. Percentiles can + be specified like ``pXX.YY``. span: str Overlap or full; if overlap, statitstics are computed on common time- span; if full, statistics are computed on full time spans, ignoring diff --git a/tests/sample_data/multimodel_statistics/test_multimodel.py b/tests/sample_data/multimodel_statistics/test_multimodel.py index a62498d672..485328d85b 100644 --- a/tests/sample_data/multimodel_statistics/test_multimodel.py +++ b/tests/sample_data/multimodel_statistics/test_multimodel.py @@ -1,7 +1,7 @@ """Test using sample data for :func:`esmvalcore.preprocessor._multimodel`.""" import pickle -import sys +import platform from itertools import groupby from pathlib import Path @@ -24,7 +24,11 @@ reason='Cannot calculate statistics with single cube in list')), '365_day', 'gregorian', - 'proleptic_gregorian', + pytest.param( + 'proleptic_gregorian', + marks=pytest.mark.xfail( + raises=iris.exceptions.MergeError, + reason='https://github.com/ESMValGroup/ESMValCore/issues/956')), pytest.param( 'julian', marks=pytest.mark.skip( @@ -39,7 +43,30 @@ def assert_array_almost_equal(this, other): if np.ma.isMaskedArray(this) or np.ma.isMaskedArray(other): np.testing.assert_array_equal(this.mask, other.mask) - np.testing.assert_array_almost_equal(this, other) + np.testing.assert_allclose(this, other) + + +def assert_coords_equal(this: list, other: list): + """Assert coords list `this` equals coords list `other`.""" + for this_coord, other_coord in zip(this, other): + np.testing.assert_equal(this_coord.points, other_coord.points) + assert this_coord.var_name == other_coord.var_name + assert this_coord.standard_name == other_coord.standard_name + assert this_coord.units == other_coord.units + + +def assert_metadata_equal(this, other): + """Assert metadata `this` are equal to metadata `other`.""" + assert this.standard_name == other.standard_name + assert this.long_name == other.long_name + assert this.var_name == other.var_name + assert this.units == other.units + + +def fix_metadata(cubes): + """Fix metadata.""" + for cube in cubes: + cube.coord('air_pressure').bounds = None def preprocess_data(cubes, time_slice: dict = None): @@ -67,13 +94,10 @@ def get_cache_key(value): If this doesn't avoid problems with unpickling the cached data, manually clean the pytest cache with the command `pytest --cache-clear`. """ - return ' '.join([ - str(value), - iris.__version__, - np.__version__, - sys.version, - f"rev-{TEST_REVISION}", - ]) + py_version = platform.python_version() + return (f'{value}_iris-{iris.__version__}_' + f'numpy-{np.__version__}_python-{py_version}' + f'rev-{TEST_REVISION}') @pytest.fixture(scope="module") @@ -102,6 +126,8 @@ def timeseries_cubes_month(request): request.config.cache.set(cache_key, pickle.dumps(cubes).decode('latin1')) + fix_metadata(cubes) + return cubes @@ -132,6 +158,8 @@ def timeseries_cubes_day(request): request.config.cache.set(cache_key, pickle.dumps(cubes).decode('latin1')) + fix_metadata(cubes) + def calendar(cube): return cube.coord('time').units.calendar @@ -172,17 +200,10 @@ def multimodel_regression_test(cubes, span, name): filename = Path(__file__).with_name(f'{name}-{span}-{statistic}.nc') if filename.exists(): reference_cube = iris.load_cube(str(filename)) - assert_array_almost_equal(result_cube.data, reference_cube.data) - - # Compare coords - for this_coord, other_coord in zip(result_cube.coords(), - reference_cube.coords()): - assert this_coord == other_coord - # remove Conventions which are added by Iris on save - reference_cube.attributes.pop('Conventions', None) - - assert reference_cube.metadata == result_cube.metadata + assert_array_almost_equal(result_cube.data, reference_cube.data) + assert_metadata_equal(result_cube.metadata, reference_cube.metadata) + assert_coords_equal(result_cube.coords(), reference_cube.coords()) else: # The test will fail if no regression data are available. @@ -190,6 +211,9 @@ def multimodel_regression_test(cubes, span, name): raise RuntimeError(f'Wrote reference data to {filename.absolute()}') +@pytest.mark.xfail( + raises=iris.exceptions.MergeError, + reason='https://github.com/ESMValGroup/ESMValCore/issues/956') @pytest.mark.use_sample_data @pytest.mark.parametrize('span', SPAN_PARAMS) def test_multimodel_regression_month(timeseries_cubes_month, span): @@ -228,8 +252,11 @@ def test_multimodel_no_vertical_dimension(timeseries_cubes_month): @pytest.mark.use_sample_data @pytest.mark.xfail( - 'iris.exceptions.CoordinateNotFoundError', - reason='https://github.com/ESMValGroup/ESMValCore/issues/891') + raises=iris.exceptions.MergeError, + reason='https://github.com/ESMValGroup/ESMValCore/issues/956') +# @pytest.mark.xfail( +# raises=iris.exceptions.CoordinateNotFoundError, +# reason='https://github.com/ESMValGroup/ESMValCore/issues/891') def test_multimodel_no_horizontal_dimension(timeseries_cubes_month): """Test statistic without horizontal dimension using monthly data.""" span = 'full' @@ -252,7 +279,7 @@ def test_multimodel_only_time_dimension(timeseries_cubes_month): @pytest.mark.use_sample_data @pytest.mark.xfail( - 'ValueError', + raises=ValueError, reason='https://github.com/ESMValGroup/ESMValCore/issues/890') def test_multimodel_no_time_dimension(timeseries_cubes_month): """Test statistic without time dimension using monthly data.""" diff --git a/tests/unit/preprocessor/_multimodel/test_multimodel.py b/tests/unit/preprocessor/_multimodel/test_multimodel.py index c63d3bb39d..dc9092eb2f 100644 --- a/tests/unit/preprocessor/_multimodel/test_multimodel.py +++ b/tests/unit/preprocessor/_multimodel/test_multimodel.py @@ -2,6 +2,7 @@ from datetime import datetime +import dask.array as da import iris import numpy as np import pytest @@ -56,6 +57,7 @@ def generate_cube_from_dates( fill_val=1, len_data=3, var_name=None, + lazy=False, ): """Generate test cube from list of dates / frequency specification. @@ -80,61 +82,63 @@ def generate_cube_from_dates( if isinstance(dates, str): time = timecoord(dates, calendar, offset=offset, num=len_data) else: + len_data = len(dates) unit = Unit(offset, calendar=calendar) time = iris.coords.DimCoord(unit.date2num(dates), standard_name='time', units=unit) - return Cube((fill_val, ) * len_data, - dim_coords_and_dims=[(time, 0)], - var_name=var_name) + data = np.array((fill_val, ) * len_data) + if lazy: + data = da.from_array(data) -def get_cubes_for_validation_test(frequency): + return Cube(data, dim_coords_and_dims=[(time, 0)], var_name=var_name) + + +def get_cubes_for_validation_test(frequency, lazy=False): """Set up cubes used for testing multimodel statistics.""" # Simple 1d cube with standard time cord - cube1 = generate_cube_from_dates(frequency) + cube1 = generate_cube_from_dates(frequency, lazy=lazy) # Cube with masked data cube2 = cube1.copy() - cube2.data = np.ma.array([5, 5, 5], mask=[True, False, False]) + data2 = np.ma.array([5, 5, 5], mask=[True, False, False]) + if lazy: + data2 = da.from_array(data2) + cube2.data = data2 # Cube with deviating time coord cube3 = generate_cube_from_dates(frequency, calendar='360_day', offset='days since 1950-01-01', len_data=2, - fill_val=9) + fill_val=9, + lazy=lazy) return [cube1, cube2, cube3] VALIDATION_DATA_SUCCESS = ( ('full', 'mean', (5, 5, 3)), - pytest.param( - 'full', - 'std', (5.656854249492381, 4, 2.8284271247461903), - marks=pytest.mark.xfail( - raises=AssertionError, - reason='https://github.com/ESMValGroup/ESMValCore/issues/1024')), + ('full', 'std_dev', (5.656854249492381, 4, 2.8284271247461903)), + ('full', 'std', (5.656854249492381, 4, 2.8284271247461903)), ('full', 'min', (1, 1, 1)), ('full', 'max', (9, 9, 5)), ('full', 'median', (5, 5, 3)), ('full', 'p50', (5, 5, 3)), ('full', 'p99.5', (8.96, 8.96, 4.98)), + ('full', 'peak', ([9], [9], [5])), ('overlap', 'mean', (5, 5)), - pytest.param( - 'full', - 'std', (5.656854249492381, 4), - marks=pytest.mark.xfail( - raises=AssertionError, - reason='https://github.com/ESMValGroup/ESMValCore/issues/1024')), + ('overlap', 'std_dev', (5.656854249492381, 4)), + ('overlap', 'std', (5.656854249492381, 4)), ('overlap', 'min', (1, 1)), ('overlap', 'max', (9, 9)), ('overlap', 'median', (5, 5)), ('overlap', 'p50', (5, 5)), ('overlap', 'p99.5', (8.96, 8.96)), + ('overlap', 'peak', ([9], [9])), # test multiple statistics ('overlap', ('min', 'max'), ((1, 1), (9, 9))), ('full', ('min', 'max'), ((1, 1, 1), (9, 9, 5))), @@ -144,14 +148,7 @@ def get_cubes_for_validation_test(frequency): @pytest.mark.parametrize('frequency', FREQUENCY_OPTIONS) @pytest.mark.parametrize('span, statistics, expected', VALIDATION_DATA_SUCCESS) def test_multimodel_statistics(frequency, span, statistics, expected): - """High level test for multicube statistics function. - - - Should work for multiple data frequencies - - Should be able to deal with multiple statistics - - Should work for both span arguments - - Should deal correctly with different mask options - - Return type should be a dict with all requested statistics as keys - """ + """High level test for multicube statistics function.""" cubes = get_cubes_for_validation_test(frequency) if isinstance(statistics, str): @@ -165,10 +162,83 @@ def test_multimodel_statistics(frequency, span, statistics, expected): for i, statistic in enumerate(statistics): result_cube = result[statistic] + # make sure that temporary coord has been removed + with pytest.raises(iris.exceptions.CoordinateNotFoundError): + result_cube.coord('multi-model') + # test that real data in => real data out + assert result_cube.has_lazy_data() is False expected_data = np.ma.array(expected[i], mask=False) assert_array_allclose(result_cube.data, expected_data) +@pytest.mark.parametrize('span', SPAN_OPTIONS) +def test_lazy_data_consistent_times(span): + """Test laziness of multimodel statistics with consistent time axis.""" + cubes = ( + generate_cube_from_dates('monthly', fill_val=1, lazy=True), + generate_cube_from_dates('monthly', fill_val=3, lazy=True), + generate_cube_from_dates('monthly', fill_val=6, lazy=True), + ) + + for cube in cubes: + assert cube.has_lazy_data() + + statistic = 'sum' + statistics = (statistic, ) + + result = mm._multicube_statistics(cubes, span=span, statistics=statistics) + + result_cube = result[statistic] + assert result_cube.has_lazy_data() + + +@pytest.mark.parametrize('span', SPAN_OPTIONS) +def test_lazy_data_inconsistent_times(span): + """Test laziness of multimodel statistics with inconsistent time axis. + + This hits `_align`, which adds additional computations which must be + lazy. + """ + + cubes = ( + generate_cube_from_dates( + [datetime(1850, i, 15, 0, 0, 0) for i in range(1, 10)], lazy=True), + generate_cube_from_dates( + [datetime(1850, i, 15, 0, 0, 0) for i in range(3, 8)], lazy=True), + generate_cube_from_dates( + [datetime(1850, i, 15, 0, 0, 0) for i in range(2, 9)], lazy=True), + ) + + for cube in cubes: + assert cube.has_lazy_data() + + statistic = 'sum' + statistics = (statistic, ) + + result = mm._multicube_statistics(cubes, span=span, statistics=statistics) + + result_cube = result[statistic] + assert result_cube.has_lazy_data() + + +VALIDATION_DATA_FAIL = ( + ('percentile', ValueError), + ('wpercentile', ValueError), + ('count', TypeError), + ('proportion', TypeError), +) + + +@pytest.mark.parametrize('statistic, error', VALIDATION_DATA_FAIL) +def test_unsupported_statistics_fail(statistic, error): + """Check that unsupported statistics raise an exception.""" + cubes = get_cubes_for_validation_test('monthly') + span = 'overlap' + statistics = (statistic, ) + with pytest.raises(error): + _ = multi_model_statistics(cubes, span, statistics) + + @pytest.mark.parametrize('calendar1, calendar2, expected', ( ('360_day', '360_day', '360_day'), ('365_day', '365_day', '365_day'), @@ -193,6 +263,91 @@ def test_get_consistent_time_unit(calendar1, calendar2, expected): assert result.calendar == expected +@pytest.mark.parametrize('span', SPAN_OPTIONS) +def test_align(span): + """Test _align function.""" + + # TODO --> check that if a cube is extended, + # the extended points are masked (not NaN!) + + len_data = 3 + + cubes = [] + + for calendar in CALENDAR_OPTIONS: + cube = generate_cube_from_dates('monthly', + calendar=calendar, + len_data=3) + cubes.append(cube) + + result_cubes = mm._align(cubes, span) + + calendars = set(cube.coord('time').units.calendar for cube in result_cubes) + + assert len(calendars) == 1 + assert list(calendars)[0] == 'gregorian' + + shapes = set(cube.shape for cube in result_cubes) + + assert len(shapes) == 1 + assert tuple(shapes)[0] == (len_data, ) + + +@pytest.mark.parametrize('span', SPAN_OPTIONS) +def test_combine_same_shape(span): + """Test _combine with same shape of cubes.""" + len_data = 3 + num_cubes = 5 + cubes = [] + + for i in range(num_cubes): + cube = generate_cube_from_dates('monthly', + '360_day', + fill_val=i, + len_data=len_data) + cubes.append(cube) + + result_cube = mm._combine(cubes) + + dim_coord = result_cube.coord(mm.CONCAT_DIM) + assert dim_coord.var_name == mm.CONCAT_DIM + assert result_cube.shape == (num_cubes, len_data) + + desired = np.linspace((0, ) * len_data, + num_cubes - 1, + num=num_cubes, + dtype=int) + np.testing.assert_equal(result_cube.data, desired) + + +def test_combine_different_shape_fail(): + """Test _combine with inconsistent data.""" + num_cubes = 5 + cubes = [] + + for num in range(1, num_cubes + 1): + cube = generate_cube_from_dates('monthly', '360_day', len_data=num) + cubes.append(cube) + + with pytest.raises(iris.exceptions.MergeError): + _ = mm._combine(cubes) + + +def test_combine_inconsistent_var_names_fail(): + """Test _combine with inconsistent var names.""" + num_cubes = 5 + cubes = [] + + for num in range(num_cubes): + cube = generate_cube_from_dates('monthly', + '360_day', + var_name=f'test_var_{num}') + cubes.append(cube) + + with pytest.raises(iris.exceptions.MergeError): + _ = mm._combine(cubes) + + @pytest.mark.parametrize('span', SPAN_OPTIONS) def test_edge_case_different_time_offsets(span): cubes = ( @@ -219,10 +374,6 @@ def test_edge_case_different_time_offsets(span): desired = np.array((14., 45., 73.)) np.testing.assert_array_equal(time_coord.points, desired) - # input cubes are updated in-place - for cube in cubes: - np.testing.assert_array_equal(cube.coord('time').points, desired) - def generate_cubes_with_non_overlapping_timecoords(): """Generate sample data where time coords do not overlap.""" @@ -296,10 +447,6 @@ def test_edge_case_time_not_in_middle_of_months(span): desired = np.array((14., 45., 73.)) np.testing.assert_array_equal(time_coord.points, desired) - # input cubes are updated in-place - for cube in cubes: - np.testing.assert_array_equal(cube.coord('time').points, desired) - @pytest.mark.parametrize('span', SPAN_OPTIONS) def test_edge_case_sub_daily_data_fail(span): @@ -314,6 +461,19 @@ def test_edge_case_sub_daily_data_fail(span): _ = multi_model_statistics(cubes, span, statistics) +@pytest.mark.parametrize('span', SPAN_OPTIONS) +def test_edge_case_single_cube_fail(span): + """Test that an error is raised when a single cube is passed.""" + cube = generate_cube_from_dates('monthly') + cubes = (cube, ) + + statistic = 'min' + statistics = (statistic, ) + + with pytest.raises(ValueError): + _ = multi_model_statistics(cubes, span, statistics) + + def test_unify_time_coordinates(): """Test set common calendar.""" cube1 = generate_cube_from_dates('monthly', @@ -360,6 +520,7 @@ def test_return_products(): result1 = mm._multiproduct_statistics(products, keep_input_datasets=True, **kwargs) + result2 = mm._multiproduct_statistics(products, keep_input_datasets=False, **kwargs) From 81032a7f4159b9feb1ee349ff8bd8d39ab28d3c3 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 13:12:25 +0200 Subject: [PATCH 02/11] Remove code related to lazy evaluation --- esmvalcore/preprocessor/_multimodel.py | 57 +++----------------------- 1 file changed, 6 insertions(+), 51 deletions(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 483019c6ec..f83fa66e37 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -250,57 +250,23 @@ def _combine(cubes): return merged_cube -def rechunk(cube): - """Rechunk the cube to speed up out-of-memory computation.""" - new_chunks = {0: -1} # don't chunk along the multimodel dimension - if cube.ndim > 1: - new_chunks[1] = 'auto' # do chunk along the first subsequent dimension - - cube.data = cube.lazy_data().rechunk(new_chunks) - - logger.debug("Total data size: %s MB", cube.lazy_data().nbytes * 1e-6) - logger.debug("New chunk block size: %s MB", - cube.lazy_data().nbytes / cube.lazy_data().npartitions * 1e-6) - logger.debug("New chunk configuration: %s", cube.lazy_data()) - - def _compute_eager(cubes: list, *, operator: iris.analysis.Aggregator, **kwargs): - """Loop over slices of a cube if iris has no lazy aggregator.""" + """Compute statistics one slice at a time.""" _ = [cube.data for cube in cubes] # make sure the cubes' data are realized result_slices = [] for i in range(cubes[0].shape[0]): - single_model_slices = [cube[i] for cube in cubes - ] # maybe filter the iris warning here? + single_model_slices = [cube[i] for cube in cubes] combined_slice = _combine(single_model_slices) collapsed_slice = combined_slice.collapsed(CONCAT_DIM, operator, **kwargs) result_slices.append(collapsed_slice) result_cube = iris.cube.CubeList(result_slices).merge_cube() - - # For consistency with lazy procedure - result_cube.data = np.ma.array(result_cube.data) - - result_cube.remove_coord(CONCAT_DIM) - - return result_cube - - -def _compute_lazy(cubes: list, *, operator: iris.analysis.Aggregator, - **kwargs): - """Compute statistics using lazy iris function.""" - cube = _combine( - cubes) # this is now done for each statistic, can we avoid that? - rechunk(cube) - - # This will always return a masked array - result_cube = cube.collapsed(CONCAT_DIM, operator, **kwargs) result_cube.remove_coord(CONCAT_DIM) - for cube in cubes: - cube.remove_coord(CONCAT_DIM) + result_cube.data = np.ma.array(result_cube.data) return result_cube @@ -317,8 +283,6 @@ def _multicube_statistics(cubes, statistics, span): raise ValueError('Cannot perform multicube statistics ' 'for a single cube.') - lazy_input = bool(all(cube.has_lazy_data() for cube in cubes)) - copied_cubes = [cube.copy() for cube in cubes] # avoid modifying inputs aligned_cubes = _align(copied_cubes, span=span) @@ -327,18 +291,9 @@ def _multicube_statistics(cubes, statistics, span): logger.debug('Multicube statistics: computing: %s', statistic) operator, kwargs = _resolve_operator(statistic) - if operator.lazy_func is None: - result_cube = _compute_eager(aligned_cubes, - operator=operator, - **kwargs) - else: - result_cube = _compute_lazy(aligned_cubes, - operator=operator, - **kwargs) - - # lazy input --> lazy output - result_cube.data = result_cube.lazy_data( - ) if lazy_input else result_cube.data + result_cube = _compute_eager(aligned_cubes, + operator=operator, + **kwargs) statistics_cubes[statistic] = result_cube From 0f47ed677c7a107700e7604affc63b1d9cf95233 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 15:00:11 +0200 Subject: [PATCH 03/11] Turn off tests for lazy data --- tests/unit/preprocessor/_multimodel/test_multimodel.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit/preprocessor/_multimodel/test_multimodel.py b/tests/unit/preprocessor/_multimodel/test_multimodel.py index dc9092eb2f..570549273a 100644 --- a/tests/unit/preprocessor/_multimodel/test_multimodel.py +++ b/tests/unit/preprocessor/_multimodel/test_multimodel.py @@ -171,6 +171,7 @@ def test_multimodel_statistics(frequency, span, statistics, expected): assert_array_allclose(result_cube.data, expected_data) +@pytest.mark.xfail(reason='Lazy data not (yet) supported.') @pytest.mark.parametrize('span', SPAN_OPTIONS) def test_lazy_data_consistent_times(span): """Test laziness of multimodel statistics with consistent time axis.""" @@ -192,6 +193,7 @@ def test_lazy_data_consistent_times(span): assert result_cube.has_lazy_data() +@pytest.mark.xfail(reason='Lazy data not (yet) supported.') @pytest.mark.parametrize('span', SPAN_OPTIONS) def test_lazy_data_inconsistent_times(span): """Test laziness of multimodel statistics with inconsistent time axis. From 2e214bfc65e07bce2bd6681bac3ce01ecde5765a Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 15:01:43 +0200 Subject: [PATCH 04/11] Enforce dtype float32 in test cubes --- tests/unit/preprocessor/_multimodel/test_multimodel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit/preprocessor/_multimodel/test_multimodel.py b/tests/unit/preprocessor/_multimodel/test_multimodel.py index 570549273a..cca3f52e79 100644 --- a/tests/unit/preprocessor/_multimodel/test_multimodel.py +++ b/tests/unit/preprocessor/_multimodel/test_multimodel.py @@ -88,7 +88,7 @@ def generate_cube_from_dates( standard_name='time', units=unit) - data = np.array((fill_val, ) * len_data) + data = np.array((fill_val, ) * len_data, dtype=np.float32) if lazy: data = da.from_array(data) @@ -104,7 +104,7 @@ def get_cubes_for_validation_test(frequency, lazy=False): # Cube with masked data cube2 = cube1.copy() - data2 = np.ma.array([5, 5, 5], mask=[True, False, False]) + data2 = np.ma.array([5, 5, 5], mask=[True, False, False], dtype=np.float32) if lazy: data2 = da.from_array(data2) cube2.data = data2 From a5e0e9f50e3c76ee823b7306e475e455cdc6a295 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 16:15:57 +0200 Subject: [PATCH 05/11] Fix dtype issue --- esmvalcore/preprocessor/_multimodel.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index f83fa66e37..4ad2619e69 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -261,6 +261,10 @@ def _compute_eager(cubes: list, *, operator: iris.analysis.Aggregator, combined_slice = _combine(single_model_slices) collapsed_slice = combined_slice.collapsed(CONCAT_DIM, operator, **kwargs) + + # iris collapse seems to change dtype on some cubes, fix this + collapsed_slice.data = collapsed_slice.data.astype(np.float32) + result_slices.append(collapsed_slice) result_cube = iris.cube.CubeList(result_slices).merge_cube() From 223acd7841e8c3e2c2e8c3af0de57c50e80ec02b Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Fri, 28 May 2021 17:24:41 +0200 Subject: [PATCH 06/11] Better note on dtype fix --- esmvalcore/preprocessor/_multimodel.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 4ad2619e69..90c0c96b1a 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -262,7 +262,8 @@ def _compute_eager(cubes: list, *, operator: iris.analysis.Aggregator, collapsed_slice = combined_slice.collapsed(CONCAT_DIM, operator, **kwargs) - # iris collapse seems to change dtype on some cubes, fix this + # some iris aggregators modify dtype, see e.g. + # https://numpy.org/doc/stable/reference/generated/numpy.ma.average.html collapsed_slice.data = collapsed_slice.data.astype(np.float32) result_slices.append(collapsed_slice) From b1b12482947923d8c95f88c0f28610a30b53050a Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Mon, 7 Jun 2021 15:00:05 +0200 Subject: [PATCH 07/11] improve warning about std_dev --- esmvalcore/preprocessor/_multimodel.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 90c0c96b1a..2792fc5e88 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -44,8 +44,10 @@ def _resolve_operator(statistic: str): # special cases if statistic == 'std': logger.warning( - "Multicube statistics is aligning its behaviour with iris.analysis" - ". Please consider replacing 'std' with 'std_dev' in your code.") + "Changing statistics from specified `std` to `std_dev`, " + "since multimodel statistics is now using the iris.analysis module" + ", which also uses `std_dev`. Please consider replacing 'std' " + " with 'std_dev' in your recipe or code.") statistic = 'std_dev' elif re.match(r"^(p\d{1,2})(\.\d*)?$", statistic): From 6bfc754e9722a3075d2cf0df4cbfbef8c84bb76c Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Mon, 7 Jun 2021 15:42:12 +0200 Subject: [PATCH 08/11] Catch errors in cube alignment with nice error message --- esmvalcore/preprocessor/_multimodel.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 2792fc5e88..7e3837200a 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -147,7 +147,13 @@ def _subset(cube, time_points): begin = cube.coord('time').units.num2date(time_points[0]) end = cube.coord('time').units.num2date(time_points[-1]) constraint = iris.Constraint(time=lambda cell: begin <= cell.point <= end) - return cube.extract(constraint) + try: + return cube.extract(constraint) + except Exception as excinfo: + raise ValueError( + "Tried to align cubes in multi-model statistics, but failed for" + f" cube {cube} and time points {time_points}. Encountered the " + f"following exception: {excinfo}") def _extend(cube, time_points): @@ -191,7 +197,13 @@ def _extend(cube, time_points): cube_list = iris.cube.CubeList(cube_list) - new_cube = cube_list.concatenate_cube() + try: + new_cube = cube_list.concatenate_cube() + except Exception as excinfo: + raise ValueError( + "Tried to align cubes in multi-model statistics, but failed for" + f" cube {cube} and time points {time_points}. Encountered the " + f"following exception: {excinfo}") return new_cube From 58bd700bae3bec1f486108218c58637274c8e0df Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Mon, 7 Jun 2021 15:59:51 +0200 Subject: [PATCH 09/11] Catch concat errors in compute_eager --- esmvalcore/preprocessor/_multimodel.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 7e3837200a..6ee155d3ae 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -282,11 +282,18 @@ def _compute_eager(cubes: list, *, operator: iris.analysis.Aggregator, result_slices.append(collapsed_slice) - result_cube = iris.cube.CubeList(result_slices).merge_cube() - result_cube.remove_coord(CONCAT_DIM) + try: + result_cube = iris.cube.CubeList(result_slices).merge_cube() + except Exception as excinfo: + raise ValueError( + "Multi-model statistics failed to concatenate results into a" + f" single array. This happened for operator {operator}" + f" with computed statistics {result_slices}." + "This can happen e.g. if the calculation results in inconsistent" + f" dtypes. Encountered the following exception: {excinfo}") result_cube.data = np.ma.array(result_cube.data) - + result_cube.remove_coord(CONCAT_DIM) return result_cube From a0a4d1feace3b19207c3060f26d388f4555dd1a5 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Mon, 7 Jun 2021 16:00:15 +0200 Subject: [PATCH 10/11] Remove doc statement about lazy evaluation --- esmvalcore/preprocessor/_multimodel.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/esmvalcore/preprocessor/_multimodel.py b/esmvalcore/preprocessor/_multimodel.py index 6ee155d3ae..8469f9177b 100644 --- a/esmvalcore/preprocessor/_multimodel.py +++ b/esmvalcore/preprocessor/_multimodel.py @@ -386,9 +386,6 @@ def multi_model_statistics(products, arguments. Except for percentiles, these operators are currently not supported. - Lazy operation is supported for all statistics, except - ``median``, ``percentile``, ``gmean`` and ``hmean``. - Parameters ---------- products: list From 68a53b62a17198245efbd94e5e02a94fe9705d03 Mon Sep 17 00:00:00 2001 From: Peter Kalverla Date: Mon, 7 Jun 2021 16:21:21 +0200 Subject: [PATCH 11/11] Regenerate sample data for failing regression tests --- .../timeseries_daily_365_day-full-mean.nc | Bin 25378 -> 23299 bytes .../timeseries_daily_365_day-overlap-mean.nc | Bin 25378 -> 23299 bytes 2 files changed, 0 insertions(+), 0 deletions(-) diff --git a/tests/sample_data/multimodel_statistics/timeseries_daily_365_day-full-mean.nc b/tests/sample_data/multimodel_statistics/timeseries_daily_365_day-full-mean.nc index c473d78f49b5ca3320a423179772e5ee52a0b855..81e9e60a627758ee5e14304c4a5434366843d204 100644 GIT binary patch literal 23299 zcmeHP349bq)~_J}f(D2}I0Ru-kUJ)va)^O{2sdOoMbRXX0cH)8km;G8?w;(%0~J9; z1rhMWtm3txqN3<>7!U=IRXkYG1vmPkh#W=)5k>ZU)m<;i3<)aV>ZiXi1uy@3)vu0z z?^Shm9X)4AT6+C+E;vWX5Emz6MS}g4|NL8pm}{IE=J|RS|ELB>4@nzSE55YW#rE}> z`bzAh9?p;ao(|rt@SIx0C1Vm>QcZ-2X)0QeawO7QI+&K~ptwyUSBFttuCD$BP-R3( za!wMHRIMtDB-r|SCdP?tt1ew|-rCKhF=Tk$`7$ELm9J`Be6Ax9e~}Qi>DKCqIQkdx zc{M*>XhbA$*S=-K;qko%iMT|4TWqQ)H#c*pN1KtCBc*Ne^D2qY4CHFwE;I9TJh^?` zGdV3kQxHl14?;5<=#e(|)76 z(R!kuEkkkr_cybp>PabSHj+}xt&m!kbYaW%Awvh2TggrxGttQ;W~_siRVsgPf9keg zp3Et}9RGdY9EmzjB)Vd`^qAaiZKqOkzfp&~*RO$M+U*r(gqP z#DKB#<{J^oJlMc#Lq`r8HFoIeQJG^gV8kF1PyXYgf3`0Wtltx84o>5`7#peTLenK2{ z;0y|Dq7-J^hoo6y_5RZeYg3)VloWLu3XH%fNM>z5MJC^CO8V^Tb|C)uRy+P0gL7a)Ij#H2^D7*)mg z^!P41SdV$B|&-CRLXfvcs#nGh(b^(4Z+n1AFkdx`lp6Q_^vOZA- zd9^@}=U9uZt+FSoogj0qRjZauN;zTq;P@Y7kIO$ zQ*vctr}LSun#TL8Ty7y}=u)ndqAt;xJh)l82O}@9Ajj*=);zXb zp6tW7EYm36Y~D0uHuW%#uK7*t!()Ylxt?q$b)4w(_=YNjx>505A5PE6`S4?m z<2fH^lvPpgEL9OTNGLCt^NEitm-CTbh76oP=R*$zJ8*nLdxvir70>z0wk`w9bN+Hr zAOpv9{`TOeQatA``yv@wp5q%k(dF?yyE+kxQSmHqj{+*?!}XEgi(V+ z=p#aC78Ane&vj(&*YQTJLX59nh!HP@(c^DnT>O_XrV}mSqhrfQI!?Z-t`$G*^^ zAJ9>DPZ+Hhg>mBEFitHFQ09Sp-d9LDiu zVQl>-j2#JuC?hH`OiT-#c<4eCpR_iy=rIGg5RKUxM#io%mfaObG11V4VGJf(^kxX9 zZ_<*OBnOj{C}VxWCeIe@Y=b^(e&q%fcwTKaAalVf6k}7;k(YLRwh}_5Uk`Ll1g^VacUUke z@9u9}h+S%7#S#lw=aT;XI?^ZTc=l?tb3+}65(pEh-c7<-^0tM&?^x)u$if1m53_9$_9hD2 znW)#n#47q-@6U8?wvO9Mf_Q9G5DnJ`(TwQvKUx@^OS<{aL>ATS=vyXc6D__lguB{> zkY1$Y(i?P?WaxN=XuID+{45J=-DIQI78)Nharajy{`8EA_3KQ`>}bO5Vxr9jCQddt zG3i|sS+)-D+YJ(KQcQX!oFnb*6>ui9AIXdi~kLfNOMI;nC4=j*f@t>3HZ~9d+*0k@JLw zu1{MS)PwTwPX0jDs+)yPJ$0-dpyTG+I`-7najr|pokVx|=t|VDO&A44v+ISihG^ao79J-$NO^ta)-m{* zAU<6a#K}j35RXw=9ShNwX!a@UxBeYQhbDzMl30j+jf0q!5XAgOLEKEV^{OB`5S0us z#O@J==uFg~=;**gOegAdT>xI9ttSF#`&|H&Is|bm(TmApWD)I+qyCVn?^hPy{L;eL zYXZ1)T>$SiB_D1c#OPK*B(@IX#-$d9F1B#vP^$9~3;Wwx7)Nw8OUKITIy7nnf4@bC zDA4iOd>so9n%Hs3#762{A1Crr``J(JrxUgHdOF!OO-BQ=UH*+a&ZG90xY$Ik`%N64 zZKC;YCXP=xQ8t70|4k6-F*@FUD~KNN1TiHP#4SXz@0nQk4{Gl|6DRVh9jB824lpss zMg3u21J5D^q=&G$eF$2I5KiTgKTNYQ*G2ZKW#QdY6Y(FL*iO%vk3ZFM`5_%g4(eF` zKK1GEk$vY;dm$RsDvaxhy6+@=6YYFWN9;BobIE6JAsV=u+HeW+b}$fIt&Wnvl0!Zn>DG{y2je%33K7eTbUJbYI(=hH{4c&;AY|-#NQLLc;?UVuE^9G8ZGjQEs3```N z`e6XC><-|8>jRiQGk^_Lmm7!{?l#easMj15BZ!8R4Z0E?*c8C|+o(P%0bJKLfS!E< zsM9xq9=W7zq6O4{v^q|G?q3bmrRP*(n*hG;5WsCBfTCJd)=>?qL>uY(-LIR8MCxl+ z5RKbm;BlhyFKd{#O~WIbH0-7Oja{Lk-$NRvP@ni{OA~L#oA{ssl|y~lDWYvLeq2`1 zkCX&Ivd;J8{-7T_iNb5iPyS{gm7bf|%r_vWICOHqAJY!{amP#zi9QW4taRx8j}0{6 zYan-#fxC&O-$uG6y6y)J=_fUey(xfxL~~XIu%Aey=j@U01_m`W(17UOQz4urdixI= zo*;^wui;gq#rJ7&FVWC7k#yhKfagORAAJx){`L^Y6OE<&b-$wF3z8d3H0D|jmlJui zLwJzr)HNY&nH)kO7=l5xbv4<3m4?>_Xt+96!+@QB^dg!*C4@&k&NJ}INg=G55kjJu z#$aFhvF)fI4chzhbw@uwZ{)}1L_c~j4WW!ENekg4B5QaE=Z*;BvxWhjY8pUa^3h9& z2auE;fSah-qXs5DZs56>4XoZ`pz}=z1{2*xa^DfX@|}i9wa zY(~_zhvrY<(Q1C%W}{>g!*$&@b7-Cn@yopl9l4qV>mvSSEDL zsu#rZbA$M6{UC}OQlIr&80XV7Z({8*>eiv>W2uE}KDHpJe}4aC>Z5w;@buA9OwXn; z_0Nk6!g!0wb72^3iN+nb@FdY4_XKctaRBGhy>=QIz#BwiYMU3+*do6^Jr6Y#D}vNt z5w;`Dni9mzGlJ-y5k%jyK{OZ{MCs@t=DuxW-F6d`H_&*FXx%Oo9X>G8kbL2>t2C@_ z96+BY)VACKTzOFd%jX!_al3&NUR!lI^kEqce4I`-ypWl=Ea-!ms0sQfY08$4S z*f_|*+K=>iZ@3v@@9zBcUo1yr>dt~{VP-N61)N895G3~ z1b?Hhw3P3N6sq3BB20c6aipsJTXj-7#}F&dQE#;WP${e2C{wj)WLZC_kkfYJ{FDMJ z6<3Y8vj?P-+%Nxg*O(K}$x_^n4x~ye zFhSWbISMlr&zP@x#s!LJT&;Mz#l)7ggAJ!>fG~vY)#>mZGn{2q||jSGfpD zhlA|Er-jNJZe=5oL_J6BNB@$&%2PU7yGhfnM{vTjqj%-mD@C=XoZNGn)Bd8t&ee)D zMUBH34YQ@{NGZ9_N>W3krJhL2`nE;8f#XJ`(;{ngS?py_MgO+qJGXd1ieKbxbiZ}O zz+S@MS)A--vSr`B{lp>>)$gD|56{Rtjd$z2SCo zy_B#WMF>xP|9iyK@=`x}u7iSi`O{**{9&(03en`~p>@5QHe4;EYQit10o(rdBCH!> zw+Nd>ctwQ0BMcX5*2;bo;b+g^HMi3*bXPUm{!|TA^t@N7z5PWhY5WV$4vydEwgB02 z$iT>p^USwCvEq*OA!!54r^amQOM9(3!6(Rm>IuWD#yPT!x} zR8W3g`PX{8k$3pZyxp{+oTaJw?M5E;I^lL7KZ7^f1J>VLC(J$)Or;yF6nZ@M`h0n> ze3GB%3gkFn2Dx2O?(F=1JaYIE!>yC=mU(@|f^|+?8XUD^$NFLX&&>>6ry5OyqtNoU z^M7>IKhJlTZ`rbU`Dt(2SN+7J{&_xoJmXAfZBRWVuQJBd%VyOd0EnlGutE6)0G@2W z7M~(tPtP;s zCW>v@A$J)+hynC7QOsMmU*=j9ep?!_?O#WFgB=`ODNh#459;t$#{um5AQ>e+y(P6_ z6^aR)R?Mb3F2dggP<|w5#q8iXo9r9x;CN2q!-iW$)WEN!fin&!b~&dUew%{{9zHYv z*kFRsBJi_l)WNg{_9O&D54#lVRU9sGt>PhK{4@bG^Q!7r0m0wBX4+^+qxmm^W2YAYdPPFFX zggC0Eof+S_im!nR4V?KD)2B+y{57iblTK_Dg?tF4BHyZik*lz`s1cDo#WYY2I9#G* z{~pQC>+5;SsiVgtoNx3grs(grwN-+Rf?pm~iwF-&}a+j4tm)TPi**Gp)|Jf$- zs;YCp13pO3=KYmOU$q4$MhwMAv z+kEt!9M5+9ZixN0LghsSW$p;eMwloy* zVfF}@UZd8E^=!*%jyqI z$@GNiL>xVOd|Iac^*C7_*?^`0OAnX9m**>=nbqVc%_>A~S=j(tmS$yL_Gjth1b+aH zrDe9VNwBDv%`bqqB-J1{dx|HQB&4KlS1ehO;!VjSbT5`BS5T|GyfDZ2wXZ3v)hM$<6L8PsGm1nbzZql%$m8q@J_|!tP1ZMyHL< z9CX#VF$0pu`*P_+l$s|elRil4d>irf9>rPfkzcQj6WLQWZ_uk5mAi?`y$?nX&EOhG z_g1o_kmC<-__VW0%Pa8G@{x-=5t+Luq!%fXW(7RHsk~4nt9x2!nP`6l(M#XW^=9XK z&1~7Nm6!T6ia_k;R7330cGm$^T@O!8F%XRnjEas>#zaTFs zFjce5im#+P`%A=K@&a0ZK$DM8Ssr=J>Gt~Bp5lBqS>`2YSfqzkHd(edGtbXjjC@|7 N-m2ttn$lMa{{uKC=1~9u delta 6075 zcmb_gd0bW1w?F4{FET2E%qRlNpvWX(2qf@j(8LL~BsERJBoT-^oRM(MOjB{IX(mpT z4Q3i@^(-~Bw4XkQ%A70uIrUS^%p%L|t#i-4ioD|+j z`(k9Gw^8ThxK*_~cd@}dyk;Ic;yG+!eu=@GUJi;aiLZ^*d2aAJ^$+8;w_H5A(qMh2 zHd&`DliBguFe5jS<5*5yP{SADCa<2fRm1UU9Qci6d18keV?Yxu|qxN_-f^0x7wb~ZPSt+Smia_j?a^pN#o zQ^(ihPjRBlt#~q(ei~A!?P;Q9hcld~aMBA$zPjtMFZ;$NgAI=ji7!diw+PY+=QObg z_e8orb;CiIyROv^Wq}BxB(k>nsxE1>}`WC{$+=FVziI@;>!l< zGc68fbK`XJ%j_Nu2rxO`xKB(zp?Lt{}xZWohRD#saEXYaIdb6`uW-9w>{gYx@% zF3avAK5{@IJLI~PZAb}auXMO*+>!Xcck^7gF}z_Z&%SkkLL&n_0z#HLBCLyJJ2wa2 zlgz3;0yM699)%h?N>6FA?6Urp=xKi{@=UlVp?fs1 zAY_jtoTJ+<)5hUm0U>Fz`n-d1PA%!}VYKW#M6Y&llM3~;?A$xP-ExhIy+g}24vSNU z+J{p(Er073VAzw~-0HQ?lP_W`gFZ5j&Yk->3m@8zU&zGZ7~Z4dcno`SXh2|ersHjH z{G_2T>fyQY!SD=&BcPmN>7=vf*JNojJt&^-!`H}W29KK@4X4Z{Pc_VGpCi!&^ws+b zmpGzeZ~L^2)>4edeJ;7%f;wG4ik7zNGkB}3&g<06N$Y;7?P{0ns%2;=Vk(OrK3eC! zFZb&11XY*ocTa9mn2*C{Av!zfq{}+4|H)pPZkQ^x$}y*0IN0NWgx#M?u$~oh;Jk=y zTSS~Ed_SWO#|ZoT*I_xKtWO?)%rT*mU(B|J|kIZ5(PsaV~hVjZDqnTi6!io+sWj*#Z9LOOi@6tL=sfbgFMypmjp z-f4CCVzr7NUr}+cPDQ_GRUEn^BL9X6?=>Q>uN9HU)_3(A@#-uU+&C2*2pNeo5>jQ9 z^^kFauqH~ze8QscDwY#UK2vbG1*m67A7 zAm3fVcQ1%|cA<#39u+ZYj0mqZ5nTxfddslJ6Ynn)?%a~F+$7^+R~ais0o!K__@z|9 zBNN!r@BnX1j12cU8E?hPaAVJe_wCGLL=?^t5Is}CzTpCz*pBe1`s!XX)IKtzLuA|u zmho6Kh13=mO<82COcmDO6~q(r`YRX=1utaFaLto(ja2CKtAxvkRSZ6&V&Q+N==PzC z-D_3!Sf`?Ag^Z7;%g8Q}kzOR@+w&4m68t|T9}-@jO17_1>anFz!6HJ-Z{)r|B)mCC z#xF)0N%LgvBb07b5VJwS2eN`nLfH&*M!kZZu>!Ko1!T?_@ajSV8&(J~t)%RDTS3Y; z1qEpe-by8>5yJW^*qbh3XRd%bCIRO=2=L$qEavM4Ed4>j?Mn(Kou!S=DCoIe!Kh{h zn+FPT9VEcqMZoYt0S$X>7_!%fPZudjB2k70TZ#;sPkk{p)bCiw&hd zHZ1Y8VGiNI7#m^;dk5FyY;hgp2|0vc^Xf2-(EkZ5rV$R@wxZ`BR+PrsFpsb+NyS9M z1vjc9gv_569R5kc&Fxk!-etw{K$|Px2qr7B5s`lGKMquIxrc%hHa#-1eshh0W!Ge! zxGrNim7Gn4s)aHxFOboja>-pFcjOD`L{_hvCBT=$+5aUOhF4_Vm@6anc^R!w(f_AX znflF!LKgu?kJylQ+=epIhUW-|FJ-L$TE_M&8Mmt`O>$&pi=Mi%D3D~<=zyDnHT z@3IA>R#=ciShe4RKM6)o#^}3L^L9#jc87!~)=PMdF!_uXpPjQ}%`_|KOt<30N3EDi zSawcE7$IZ6jAFuIvRew_+Q(M-9<<_8vK8Y~tVr)~g(=gD)Jhr0D&=}Cp?nQ%m2me> z3GQ!6sOw?HA8}SZ&sj0UK+$u{f;_@*`Yp=rE5o01yNNLBgoI6m$3C%O%0UZW``Cg@ z`$+L73m)ER!9*%TuZGDu>?Px5M;S|~Xxt$jbTK2@-Hg6IW=!;D-FgP_<5*hH7(SGZ z?-?3AVzw1og!xTY{6MHpmvA#xLOyHg=~w^G2n*a6S@0R*rIi+RUqx2)mr&D1LirgH z#ivD79TPE%P(q4@e`>*1;u}OL8E?Tugr}Yqv5vrv6S1FQwTX}jpKY~Z!Db6S&$VD| zjs>|V&FDucn+NO897#3X`h};X4JT zEeiNk0={e!@J2rYW&H)Lqe?4J6|lTk#SubzR~7FON?H|cA-u?}7`M&}Khn6*5GxK5 zWD4CL^lMh*O}}Us8BI2s2{GRF`=;8EG2DjCkv4Q1V#9aCY?yyk#s|k_jQdE&0>Z9S zGGb24@FB;&HpYVYx>%7BKq1@RiiabtXq+$M#0wJsBsV`sXnIsa2Eq5d1w$z{eAB6J z5!T!6%uXrl?B=VnxbS~F zh8trk6B4u@kmJ1iYJWNIb)9hfY}Oe4@seYf3u?#>3!l^f+T(VU_H?v{W6vL{<*Nz`~@Eyt^pb+Z`2xcZMalcS2TG5w8e8eu-BIM$;evG<+$lBxPUS2}fAr{ne$UUghk@YDlY>?f>$ zpO3dUIn;I7M`)`3}zV`%D^otis)7xQfQLzW;1-+dlpo%&Gg6U1TR z0}2P`=iR$&p%sqX8@+gB^VZA-Jinb~WlOAEa0u&P>ec*C_8vaap}-MMsU(i;Ktxl3 zok)-Jd0CsJo?UCM4vL8-EVy4(x~NS?YDT3|(YTGa)8-!Rc+@LtF2at5(0eUl7NC$isjuCcX2>28{zuJAs8 zvv1O(jatEB$5Z{;(V-pK+Th-7Sh`=w$Wr@O^v5f}RqvK5dF-R7ixx;7BPTdH_U+tT z$@uJE+0Kx6DjHp+vpAopbWHPbxXw#cfn0Ylb{g5q^oWXc8Y@j}$JLz2?a?H}>W)jf zv{B;hWr%JyG{y-FX1>5Zt#N2~G(^h19%kG-+EO1qD(t?#tZ^LFBSQ0a%Qm%$Tx;W} zIOIGagKN~-BPdZd+ID}gt)K_DX{$Ha_G}rqmB2O;;JoY7M00Dt>ez=OAfJ71^kLr` zQ(0qLcfOS6m}1%NtgUpjsgLL)mYp5RW|@L1B9gRd=+ZVCil4MpSgd8`w99_YrWh zZ>j&^-06H-v$r<>}i|t=@X%InKuKIkC7+_owuCiAj3JIz{5&RemzzA*f*{mHN_v6GKgR={Q*!qHGEe8EH`#C*dTwV2)GYCK7 F{|m5cfZqTB diff --git a/tests/sample_data/multimodel_statistics/timeseries_daily_365_day-overlap-mean.nc b/tests/sample_data/multimodel_statistics/timeseries_daily_365_day-overlap-mean.nc index c473d78f49b5ca3320a423179772e5ee52a0b855..81e9e60a627758ee5e14304c4a5434366843d204 100644 GIT binary patch literal 23299 zcmeHP349bq)~_J}f(D2}I0Ru-kUJ)va)^O{2sdOoMbRXX0cH)8km;G8?w;(%0~J9; z1rhMWtm3txqN3<>7!U=IRXkYG1vmPkh#W=)5k>ZU)m<;i3<)aV>ZiXi1uy@3)vu0z z?^Shm9X)4AT6+C+E;vWX5Emz6MS}g4|NL8pm}{IE=J|RS|ELB>4@nzSE55YW#rE}> z`bzAh9?p;ao(|rt@SIx0C1Vm>QcZ-2X)0QeawO7QI+&K~ptwyUSBFttuCD$BP-R3( za!wMHRIMtDB-r|SCdP?tt1ew|-rCKhF=Tk$`7$ELm9J`Be6Ax9e~}Qi>DKCqIQkdx zc{M*>XhbA$*S=-K;qko%iMT|4TWqQ)H#c*pN1KtCBc*Ne^D2qY4CHFwE;I9TJh^?` zGdV3kQxHl14?;5<=#e(|)76 z(R!kuEkkkr_cybp>PabSHj+}xt&m!kbYaW%Awvh2TggrxGttQ;W~_siRVsgPf9keg zp3Et}9RGdY9EmzjB)Vd`^qAaiZKqOkzfp&~*RO$M+U*r(gqP z#DKB#<{J^oJlMc#Lq`r8HFoIeQJG^gV8kF1PyXYgf3`0Wtltx84o>5`7#peTLenK2{ z;0y|Dq7-J^hoo6y_5RZeYg3)VloWLu3XH%fNM>z5MJC^CO8V^Tb|C)uRy+P0gL7a)Ij#H2^D7*)mg z^!P41SdV$B|&-CRLXfvcs#nGh(b^(4Z+n1AFkdx`lp6Q_^vOZA- zd9^@}=U9uZt+FSoogj0qRjZauN;zTq;P@Y7kIO$ zQ*vctr}LSun#TL8Ty7y}=u)ndqAt;xJh)l82O}@9Ajj*=);zXb zp6tW7EYm36Y~D0uHuW%#uK7*t!()Ylxt?q$b)4w(_=YNjx>505A5PE6`S4?m z<2fH^lvPpgEL9OTNGLCt^NEitm-CTbh76oP=R*$zJ8*nLdxvir70>z0wk`w9bN+Hr zAOpv9{`TOeQatA``yv@wp5q%k(dF?yyE+kxQSmHqj{+*?!}XEgi(V+ z=p#aC78Ane&vj(&*YQTJLX59nh!HP@(c^DnT>O_XrV}mSqhrfQI!?Z-t`$G*^^ zAJ9>DPZ+Hhg>mBEFitHFQ09Sp-d9LDiu zVQl>-j2#JuC?hH`OiT-#c<4eCpR_iy=rIGg5RKUxM#io%mfaObG11V4VGJf(^kxX9 zZ_<*OBnOj{C}VxWCeIe@Y=b^(e&q%fcwTKaAalVf6k}7;k(YLRwh}_5Uk`Ll1g^VacUUke z@9u9}h+S%7#S#lw=aT;XI?^ZTc=l?tb3+}65(pEh-c7<-^0tM&?^x)u$if1m53_9$_9hD2 znW)#n#47q-@6U8?wvO9Mf_Q9G5DnJ`(TwQvKUx@^OS<{aL>ATS=vyXc6D__lguB{> zkY1$Y(i?P?WaxN=XuID+{45J=-DIQI78)Nharajy{`8EA_3KQ`>}bO5Vxr9jCQddt zG3i|sS+)-D+YJ(KQcQX!oFnb*6>ui9AIXdi~kLfNOMI;nC4=j*f@t>3HZ~9d+*0k@JLw zu1{MS)PwTwPX0jDs+)yPJ$0-dpyTG+I`-7najr|pokVx|=t|VDO&A44v+ISihG^ao79J-$NO^ta)-m{* zAU<6a#K}j35RXw=9ShNwX!a@UxBeYQhbDzMl30j+jf0q!5XAgOLEKEV^{OB`5S0us z#O@J==uFg~=;**gOegAdT>xI9ttSF#`&|H&Is|bm(TmApWD)I+qyCVn?^hPy{L;eL zYXZ1)T>$SiB_D1c#OPK*B(@IX#-$d9F1B#vP^$9~3;Wwx7)Nw8OUKITIy7nnf4@bC zDA4iOd>so9n%Hs3#762{A1Crr``J(JrxUgHdOF!OO-BQ=UH*+a&ZG90xY$Ik`%N64 zZKC;YCXP=xQ8t70|4k6-F*@FUD~KNN1TiHP#4SXz@0nQk4{Gl|6DRVh9jB824lpss zMg3u21J5D^q=&G$eF$2I5KiTgKTNYQ*G2ZKW#QdY6Y(FL*iO%vk3ZFM`5_%g4(eF` zKK1GEk$vY;dm$RsDvaxhy6+@=6YYFWN9;BobIE6JAsV=u+HeW+b}$fIt&Wnvl0!Zn>DG{y2je%33K7eTbUJbYI(=hH{4c&;AY|-#NQLLc;?UVuE^9G8ZGjQEs3```N z`e6XC><-|8>jRiQGk^_Lmm7!{?l#easMj15BZ!8R4Z0E?*c8C|+o(P%0bJKLfS!E< zsM9xq9=W7zq6O4{v^q|G?q3bmrRP*(n*hG;5WsCBfTCJd)=>?qL>uY(-LIR8MCxl+ z5RKbm;BlhyFKd{#O~WIbH0-7Oja{Lk-$NRvP@ni{OA~L#oA{ssl|y~lDWYvLeq2`1 zkCX&Ivd;J8{-7T_iNb5iPyS{gm7bf|%r_vWICOHqAJY!{amP#zi9QW4taRx8j}0{6 zYan-#fxC&O-$uG6y6y)J=_fUey(xfxL~~XIu%Aey=j@U01_m`W(17UOQz4urdixI= zo*;^wui;gq#rJ7&FVWC7k#yhKfagORAAJx){`L^Y6OE<&b-$wF3z8d3H0D|jmlJui zLwJzr)HNY&nH)kO7=l5xbv4<3m4?>_Xt+96!+@QB^dg!*C4@&k&NJ}INg=G55kjJu z#$aFhvF)fI4chzhbw@uwZ{)}1L_c~j4WW!ENekg4B5QaE=Z*;BvxWhjY8pUa^3h9& z2auE;fSah-qXs5DZs56>4XoZ`pz}=z1{2*xa^DfX@|}i9wa zY(~_zhvrY<(Q1C%W}{>g!*$&@b7-Cn@yopl9l4qV>mvSSEDL zsu#rZbA$M6{UC}OQlIr&80XV7Z({8*>eiv>W2uE}KDHpJe}4aC>Z5w;@buA9OwXn; z_0Nk6!g!0wb72^3iN+nb@FdY4_XKctaRBGhy>=QIz#BwiYMU3+*do6^Jr6Y#D}vNt z5w;`Dni9mzGlJ-y5k%jyK{OZ{MCs@t=DuxW-F6d`H_&*FXx%Oo9X>G8kbL2>t2C@_ z96+BY)VACKTzOFd%jX!_al3&NUR!lI^kEqce4I`-ypWl=Ea-!ms0sQfY08$4S z*f_|*+K=>iZ@3v@@9zBcUo1yr>dt~{VP-N61)N895G3~ z1b?Hhw3P3N6sq3BB20c6aipsJTXj-7#}F&dQE#;WP${e2C{wj)WLZC_kkfYJ{FDMJ z6<3Y8vj?P-+%Nxg*O(K}$x_^n4x~ye zFhSWbISMlr&zP@x#s!LJT&;Mz#l)7ggAJ!>fG~vY)#>mZGn{2q||jSGfpD zhlA|Er-jNJZe=5oL_J6BNB@$&%2PU7yGhfnM{vTjqj%-mD@C=XoZNGn)Bd8t&ee)D zMUBH34YQ@{NGZ9_N>W3krJhL2`nE;8f#XJ`(;{ngS?py_MgO+qJGXd1ieKbxbiZ}O zz+S@MS)A--vSr`B{lp>>)$gD|56{Rtjd$z2SCo zy_B#WMF>xP|9iyK@=`x}u7iSi`O{**{9&(03en`~p>@5QHe4;EYQit10o(rdBCH!> zw+Nd>ctwQ0BMcX5*2;bo;b+g^HMi3*bXPUm{!|TA^t@N7z5PWhY5WV$4vydEwgB02 z$iT>p^USwCvEq*OA!!54r^amQOM9(3!6(Rm>IuWD#yPT!x} zR8W3g`PX{8k$3pZyxp{+oTaJw?M5E;I^lL7KZ7^f1J>VLC(J$)Or;yF6nZ@M`h0n> ze3GB%3gkFn2Dx2O?(F=1JaYIE!>yC=mU(@|f^|+?8XUD^$NFLX&&>>6ry5OyqtNoU z^M7>IKhJlTZ`rbU`Dt(2SN+7J{&_xoJmXAfZBRWVuQJBd%VyOd0EnlGutE6)0G@2W z7M~(tPtP;s zCW>v@A$J)+hynC7QOsMmU*=j9ep?!_?O#WFgB=`ODNh#459;t$#{um5AQ>e+y(P6_ z6^aR)R?Mb3F2dggP<|w5#q8iXo9r9x;CN2q!-iW$)WEN!fin&!b~&dUew%{{9zHYv z*kFRsBJi_l)WNg{_9O&D54#lVRU9sGt>PhK{4@bG^Q!7r0m0wBX4+^+qxmm^W2YAYdPPFFX zggC0Eof+S_im!nR4V?KD)2B+y{57iblTK_Dg?tF4BHyZik*lz`s1cDo#WYY2I9#G* z{~pQC>+5;SsiVgtoNx3grs(grwN-+Rf?pm~iwF-&}a+j4tm)TPi**Gp)|Jf$- zs;YCp13pO3=KYmOU$q4$MhwMAv z+kEt!9M5+9ZixN0LghsSW$p;eMwloy* zVfF}@UZd8E^=!*%jyqI z$@GNiL>xVOd|Iac^*C7_*?^`0OAnX9m**>=nbqVc%_>A~S=j(tmS$yL_Gjth1b+aH zrDe9VNwBDv%`bqqB-J1{dx|HQB&4KlS1ehO;!VjSbT5`BS5T|GyfDZ2wXZ3v)hM$<6L8PsGm1nbzZql%$m8q@J_|!tP1ZMyHL< z9CX#VF$0pu`*P_+l$s|elRil4d>irf9>rPfkzcQj6WLQWZ_uk5mAi?`y$?nX&EOhG z_g1o_kmC<-__VW0%Pa8G@{x-=5t+Luq!%fXW(7RHsk~4nt9x2!nP`6l(M#XW^=9XK z&1~7Nm6!T6ia_k;R7330cGm$^T@O!8F%XRnjEas>#zaTFs zFjce5im#+P`%A=K@&a0ZK$DM8Ssr=J>Gt~Bp5lBqS>`2YSfqzkHd(edGtbXjjC@|7 N-m2ttn$lMa{{uKC=1~9u delta 6075 zcmb_gd0bW1w?F4{FET2E%qRlNpvWX(2qf@j(8LL~BsERJBoT-^oRM(MOjB{IX(mpT z4Q3i@^(-~Bw4XkQ%A70uIrUS^%p%L|t#i-4ioD|+j z`(k9Gw^8ThxK*_~cd@}dyk;Ic;yG+!eu=@GUJi;aiLZ^*d2aAJ^$+8;w_H5A(qMh2 zHd&`DliBguFe5jS<5*5yP{SADCa<2fRm1UU9Qci6d18keV?Yxu|qxN_-f^0x7wb~ZPSt+Smia_j?a^pN#o zQ^(ihPjRBlt#~q(ei~A!?P;Q9hcld~aMBA$zPjtMFZ;$NgAI=ji7!diw+PY+=QObg z_e8orb;CiIyROv^Wq}BxB(k>nsxE1>}`WC{$+=FVziI@;>!l< zGc68fbK`XJ%j_Nu2rxO`xKB(zp?Lt{}xZWohRD#saEXYaIdb6`uW-9w>{gYx@% zF3avAK5{@IJLI~PZAb}auXMO*+>!Xcck^7gF}z_Z&%SkkLL&n_0z#HLBCLyJJ2wa2 zlgz3;0yM699)%h?N>6FA?6Urp=xKi{@=UlVp?fs1 zAY_jtoTJ+<)5hUm0U>Fz`n-d1PA%!}VYKW#M6Y&llM3~;?A$xP-ExhIy+g}24vSNU z+J{p(Er073VAzw~-0HQ?lP_W`gFZ5j&Yk->3m@8zU&zGZ7~Z4dcno`SXh2|ersHjH z{G_2T>fyQY!SD=&BcPmN>7=vf*JNojJt&^-!`H}W29KK@4X4Z{Pc_VGpCi!&^ws+b zmpGzeZ~L^2)>4edeJ;7%f;wG4ik7zNGkB}3&g<06N$Y;7?P{0ns%2;=Vk(OrK3eC! zFZb&11XY*ocTa9mn2*C{Av!zfq{}+4|H)pPZkQ^x$}y*0IN0NWgx#M?u$~oh;Jk=y zTSS~Ed_SWO#|ZoT*I_xKtWO?)%rT*mU(B|J|kIZ5(PsaV~hVjZDqnTi6!io+sWj*#Z9LOOi@6tL=sfbgFMypmjp z-f4CCVzr7NUr}+cPDQ_GRUEn^BL9X6?=>Q>uN9HU)_3(A@#-uU+&C2*2pNeo5>jQ9 z^^kFauqH~ze8QscDwY#UK2vbG1*m67A7 zAm3fVcQ1%|cA<#39u+ZYj0mqZ5nTxfddslJ6Ynn)?%a~F+$7^+R~ais0o!K__@z|9 zBNN!r@BnX1j12cU8E?hPaAVJe_wCGLL=?^t5Is}CzTpCz*pBe1`s!XX)IKtzLuA|u zmho6Kh13=mO<82COcmDO6~q(r`YRX=1utaFaLto(ja2CKtAxvkRSZ6&V&Q+N==PzC z-D_3!Sf`?Ag^Z7;%g8Q}kzOR@+w&4m68t|T9}-@jO17_1>anFz!6HJ-Z{)r|B)mCC z#xF)0N%LgvBb07b5VJwS2eN`nLfH&*M!kZZu>!Ko1!T?_@ajSV8&(J~t)%RDTS3Y; z1qEpe-by8>5yJW^*qbh3XRd%bCIRO=2=L$qEavM4Ed4>j?Mn(Kou!S=DCoIe!Kh{h zn+FPT9VEcqMZoYt0S$X>7_!%fPZudjB2k70TZ#;sPkk{p)bCiw&hd zHZ1Y8VGiNI7#m^;dk5FyY;hgp2|0vc^Xf2-(EkZ5rV$R@wxZ`BR+PrsFpsb+NyS9M z1vjc9gv_569R5kc&Fxk!-etw{K$|Px2qr7B5s`lGKMquIxrc%hHa#-1eshh0W!Ge! zxGrNim7Gn4s)aHxFOboja>-pFcjOD`L{_hvCBT=$+5aUOhF4_Vm@6anc^R!w(f_AX znflF!LKgu?kJylQ+=epIhUW-|FJ-L$TE_M&8Mmt`O>$&pi=Mi%D3D~<=zyDnHT z@3IA>R#=ciShe4RKM6)o#^}3L^L9#jc87!~)=PMdF!_uXpPjQ}%`_|KOt<30N3EDi zSawcE7$IZ6jAFuIvRew_+Q(M-9<<_8vK8Y~tVr)~g(=gD)Jhr0D&=}Cp?nQ%m2me> z3GQ!6sOw?HA8}SZ&sj0UK+$u{f;_@*`Yp=rE5o01yNNLBgoI6m$3C%O%0UZW``Cg@ z`$+L73m)ER!9*%TuZGDu>?Px5M;S|~Xxt$jbTK2@-Hg6IW=!;D-FgP_<5*hH7(SGZ z?-?3AVzw1og!xTY{6MHpmvA#xLOyHg=~w^G2n*a6S@0R*rIi+RUqx2)mr&D1LirgH z#ivD79TPE%P(q4@e`>*1;u}OL8E?Tugr}Yqv5vrv6S1FQwTX}jpKY~Z!Db6S&$VD| zjs>|V&FDucn+NO897#3X`h};X4JT zEeiNk0={e!@J2rYW&H)Lqe?4J6|lTk#SubzR~7FON?H|cA-u?}7`M&}Khn6*5GxK5 zWD4CL^lMh*O}}Us8BI2s2{GRF`=;8EG2DjCkv4Q1V#9aCY?yyk#s|k_jQdE&0>Z9S zGGb24@FB;&HpYVYx>%7BKq1@RiiabtXq+$M#0wJsBsV`sXnIsa2Eq5d1w$z{eAB6J z5!T!6%uXrl?B=VnxbS~F zh8trk6B4u@kmJ1iYJWNIb)9hfY}Oe4@seYf3u?#>3!l^f+T(VU_H?v{W6vL{<*Nz`~@Eyt^pb+Z`2xcZMalcS2TG5w8e8eu-BIM$;evG<+$lBxPUS2}fAr{ne$UUghk@YDlY>?f>$ zpO3dUIn;I7M`)`3}zV`%D^otis)7xQfQLzW;1-+dlpo%&Gg6U1TR z0}2P`=iR$&p%sqX8@+gB^VZA-Jinb~WlOAEa0u&P>ec*C_8vaap}-MMsU(i;Ktxl3 zok)-Jd0CsJo?UCM4vL8-EVy4(x~NS?YDT3|(YTGa)8-!Rc+@LtF2at5(0eUl7NC$isjuCcX2>28{zuJAs8 zvv1O(jatEB$5Z{;(V-pK+Th-7Sh`=w$Wr@O^v5f}RqvK5dF-R7ixx;7BPTdH_U+tT z$@uJE+0Kx6DjHp+vpAopbWHPbxXw#cfn0Ylb{g5q^oWXc8Y@j}$JLz2?a?H}>W)jf zv{B;hWr%JyG{y-FX1>5Zt#N2~G(^h19%kG-+EO1qD(t?#tZ^LFBSQ0a%Qm%$Tx;W} zIOIGagKN~-BPdZd+ID}gt)K_DX{$Ha_G}rqmB2O;;JoY7M00Dt>ez=OAfJ71^kLr` zQ(0qLcfOS6m}1%NtgUpjsgLL)mYp5RW|@L1B9gRd=+ZVCil4MpSgd8`w99_YrWh zZ>j&^-06H-v$r<>}i|t=@X%InKuKIkC7+_owuCiAj3JIz{5&RemzzA*f*{mHN_v6GKgR={Q*!qHGEe8EH`#C*dTwV2)GYCK7 F{|m5cfZqTB