From 2f8c293d5f83b51832faeddde92348266d4eb380 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 16 Feb 2024 17:20:27 -0500 Subject: [PATCH 01/13] _fit_start for `gamma`, `fisk` can accept negative values --- xclim/indices/stats.py | 43 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 262a6065b..01597374e 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -508,17 +508,26 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: return (chat,), {"loc": loc, "scale": scale} if dist in ["gamma"]: - x_pos = x[x > 0] + # not sure if the approximation holds for xmin < 0 + xmin = x.min() + x_pos = x - (xmin if xmin <= 0 else 0) + x_pos = x_pos[x_pos > 0] m = x_pos.mean() log_of_mean = np.log(m) mean_of_logs = np.log(x_pos).mean() a = log_of_mean - mean_of_logs alpha = (1 + np.sqrt(1 + 4 * a / 3)) / (4 * a) beta = m / alpha - return (alpha,), {"scale": beta} + kwargs = {"scale": beta} + if xmin < 0: + kwargs["loc"] = xmin + return (alpha,), kwargs if dist in ["fisk"]: - x_pos = x[x > 0] + # not sure if the approximation holds for xmin < 0 + xmin = x.min() + x_pos = x - (xmin if xmin <= 0 else 0) + x_pos = x_pos[x_pos > 0] m = x_pos.mean() v = x_pos.var() # pdf = (beta/alpha) (x/alpha)^{beta-1}/ (1+(x/alpha)^{beta})^2 @@ -528,7 +537,10 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: # In the large beta limit, f_1 -> 1 and f_1/sqrt(f_2) -> 0.56*beta - 0.25 # Solve for alpha and beta below: alpha, beta = m, (1 / 0.56) * (m / np.sqrt(v) + 1 / 4) - return (beta,), {"scale": alpha} + kwargs = {"scale": alpha} + if xmin < 0: + kwargs["loc"] = xmin + return (beta,), kwargs return (), {} @@ -554,6 +566,7 @@ def _dist_method_1D(*args, dist: str, function: str, **kwargs: Any) -> xr.DataAr array_like """ dist = get_dist(dist) + return getattr(dist, function)(*args, **kwargs) @@ -592,6 +605,27 @@ def dist_method( # Typically the data to be transformed arg = [arg] if arg is not None else [] args = arg + [fit_params.sel(dparams=dp) for dp in fit_params.dparams.values] + # print(kwargs) + # print(fit_params.dparams.values) + # print(len(args)) + # for ii in range(4): + # print(ii) + # print(args[ii]) + # import scipy + # return xr.apply_ufunc( + # scipy.stats.gamma.cdf(), + # *args, + # ) + # import scipy + # if len(args)==4: + # for ii, arg in enumerate(args): + # arg.to_netcdf(f"/home/eridup1/tmp/tmp_{ii}.nc") + # return xr.apply_ufunc( + # scipy.stats.gamma.cdf, + # *args, + # output_dtypes=[float], + # dask="parallelized", + # ) return xr.apply_ufunc( _dist_method_1D, *args, @@ -781,7 +815,6 @@ def reindex_time(da, da_ref): params, probs_of_zero = (reindex_time(dax, da) for dax in [params, probs_of_zero]) dist_probs = dist_method("cdf", params, da) probs = probs_of_zero + ((1 - probs_of_zero) * dist_probs) - params_norm = xr.DataArray( [0, 1], dims=["dparams"], From c3e4ff9339dbe2bd5893773c5623d33f041e1ce0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 16 Feb 2024 21:48:58 -0500 Subject: [PATCH 02/13] remove offset, standardize_index as general-use func --- tests/test_indices.py | 3 +- xclim/indices/_agro.py | 122 +++++++---------------------------------- xclim/indices/stats.py | 118 ++++++++++++++++++++++++++++++--------- 3 files changed, 115 insertions(+), 128 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index e815b1766..e13f2ac96 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -695,6 +695,8 @@ def test_standardized_precipitation_evapotranspiration_index( tas = tasmax - 2.5 tasmin = tasmax - 5 wb = xci.water_budget(pr, None, tasmin, tasmax, tas) + # offset like in climate indices library + wb = wb + convert_units_to("1 mm/d", wb) params = xci.stats.standardized_index_fit_params( wb.sel(time=slice("1950", "1980")), @@ -702,7 +704,6 @@ def test_standardized_precipitation_evapotranspiration_index( window=window, dist=dist, method=method, - offset="1 mm/d", ) spei = xci.standardized_precipitation_evapotranspiration_index( wb.sel(time=slice("1998", "2000")), params=params diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index cdfff68e7..0d1fc2d22 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,8 +1,6 @@ # noqa: D100 from __future__ import annotations -import warnings - import numpy as np import xarray @@ -24,11 +22,7 @@ ) from xclim.indices.generic import aggregate_between_dates, get_zones from xclim.indices.helpers import _gather_lat, day_lengths -from xclim.indices.stats import ( - preprocess_standardized_index, - standardized_index, - standardized_index_fit_params, -) +from xclim.indices.stats import standardized_index # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1094,12 +1088,10 @@ def _get_rain_season(pram): @declare_units( pr="[precipitation]", - pr_cal="[precipitation]", params="[]", ) def standardized_precipitation_index( pr: xarray.DataArray, - pr_cal: Quantified | None = None, freq: str | None = "MS", window: int = 1, dist: str = "gamma", @@ -1115,9 +1107,6 @@ def standardized_precipitation_index( ---------- pr : xarray.DataArray Daily precipitation. - pr_cal : xarray.DataArray, optional - Daily precipitation used for calibration. Usually this is a temporal subset of `pr` over some reference period. - This option will be removed in xclim >=0.47.0. Two behaviour will be possible (see below) freq : str, optional Resampling frequency. A monthly or daily frequency is expected. Option `None` assumes that desired resampling has already been applied input dataset and will skip the resampling step. @@ -1189,75 +1178,30 @@ def standardized_precipitation_index( ---------- :cite:cts:`mckee_relationship_1993` """ - if params is not None and pr_cal is None: - freq, window, indexer = ( - params.attrs[s] for s in ["freq", "window", "time_indexer"] - ) - # Unpack attrs to None and {} if needed - freq = None if freq == "" else freq - indexer = {} if indexer[0] == "" else {indexer[0]: indexer[1:]} - if cal_start or cal_end: - warnings.warn( - "Expected either `cal_{start|end}` or `params`, got both. The `params` input overrides other inputs." - "If `cal_start`, `cal_end`, `freq`, `window`, and/or `dist` were given as input, they will be ignored." + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + if dist in dist_methods.keys(): + if method not in dist_methods[dist]: + raise NotImplementedError( + f"{method} method is not implemented for {dist} distribution" ) - - if pr_cal is not None: - warnings.warn( - "Inputting a calibration array will be deprecated in xclim>=0.47.0. " - "For example, if `pr_cal` is a subset of `pr`, then instead of say:\n" - "`standardized_precipitation_index(pr=pr,pr_cal=pr.sel(time=slice(t0,t1)),...)`,\n" - "one can call:\n" - "`standardized_precipitation_index(pr=pr,cal_range=(t0,t1),...).\n" - "If for some reason `pr_cal` is not a subset of `pr`, then the following approach will still be possible:\n" - "`params = standardized_index_fit_params(da=pr_cal, freq=freq, window=window, dist=dist, method=method)`.\n" - "`spi = standardized_precipitation_index(pr=pr, params=params)`.\n" - "This approach can be used in both scenarios to break up the computations in two," - "i.e. get params, then compute standardized indices." - ) - params = standardized_index_fit_params( - pr_cal, freq=freq, window=window, dist=dist, method=method, **indexer - ) - - pr, _ = preprocess_standardized_index(pr, freq=freq, window=window, **indexer) - if params is None: - params = standardized_index_fit_params( - pr.sel(time=slice(cal_start, cal_end)), - freq=None, - window=1, - dist=dist, - method=method, - ) - - # If params only contains a subset of main dataset time grouping - # (e.g. 8/12 months, etc.), it needs to be broadcasted - template = pr.groupby(params.attrs["group"]).first() - paramsd = {k: v for k, v in params.sizes.items() if k != "dparams"} - if paramsd != template.sizes: - params = params.broadcast_like(template) - - spi = standardized_index(pr, params) - spi.attrs = params.attrs - spi.attrs["freq"] = (freq or xarray.infer_freq(spi.time)) or "undefined" - spi.attrs["window"] = window - spi.attrs["units"] = "" + else: + raise NotImplementedError(f"{dist} distribution is not implemented yet") + spi = standardized_index( + pr, freq, window, dist, method, cal_start, cal_end, params, **indexer + ) return spi @declare_units( wb="[precipitation]", - wb_cal="[precipitation]", - offset="[precipitation]", params="[]", ) def standardized_precipitation_evapotranspiration_index( wb: xarray.DataArray, - wb_cal: Quantified | None = None, freq: str | None = "MS", window: int = 1, dist: str = "gamma", method: str = "APP", - offset: Quantified = "1.000 mm/d", cal_start: DateStr | None = None, cal_end: DateStr | None = None, params: Quantified | None = None, @@ -1273,9 +1217,6 @@ def standardized_precipitation_evapotranspiration_index( ---------- wb : xarray.DataArray Daily water budget (pr - pet). - wb_cal : xarray.DataArray, optional - Daily water budget used for calibration. Usually this is a temporal subset of `wb` over some reference period. - This option will be removed in xclim >=0.47.0. Two behaviours will be possible (see below). freq : str, optional Resampling frequency. A monthly or daily frequency is expected. Option `None` assumes that desired resampling has already been applied input dataset and will skip the resampling step. @@ -1299,11 +1240,6 @@ def standardized_precipitation_evapotranspiration_index( Fit parameters. The `params` can be computed using ``xclim.indices.stats.standardized_index_fit_params`` in advance. The output can be given here as input, and it overrides other options. - offset : Quantified - For distributions bounded by zero (e.g. "gamma", "fisk"), an offset must be added to the water budget - to make sure there are no negative values. - Keep the offset as small as possible to minimize its influence on the results. - This can be given as a precipitation flux or a rate. \*\*indexer Indexing parameters to compute the indicator on a temporal subset of the data. It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`. @@ -1323,34 +1259,16 @@ def standardized_precipitation_evapotranspiration_index( See Standardized Precipitation Index (SPI) for more details on usage. """ - uses_default_offset = offset == "1.000 mm/d" - if params is not None: - params_offset = params.attrs["offset"] - if uses_default_offset is False and offset != params_offset: - warnings.warn( - "The offset in `params` differs from the input `offset`." - "Proceeding with the value given in `params`." + dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + if dist in dist_methods.keys(): + if method not in dist_methods[dist]: + raise NotImplementedError( + f"{method} method is not implemented for {dist} distribution" ) - offset = params_offset - offset = 0 if offset == "" else convert_units_to(offset, wb, context="hydro") - # Allowed distributions are constrained by the SPI function - if dist in ["gamma", "fisk"] and offset <= 0: - raise ValueError( - "The water budget must be shifted towards positive values to be used with `gamma` and `fisk` " - "distributions which are bounded by zero. A positive offset is required. Current value: " - f"{offset}{wb.attrs['units']}." - ) - # Note that the default behaviour would imply an offset for any distribution, even those distributions - # that can accommodate negative values of the water budget. This needs to be changed in future versions - # of the index. - if offset != 0: - with xarray.set_options(keep_attrs=True): - wb = wb + offset - if wb_cal is not None: - wb_cal = wb_cal + offset - - spei = standardized_precipitation_index( - wb, wb_cal, freq, window, dist, method, cal_start, cal_end, params, **indexer + else: + raise NotImplementedError(f"{dist} distribution is not implemented yet") + spei = standardized_index( + wb, freq, window, dist, method, cal_start, cal_end, params, **indexer ) return spei diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 01597374e..1f06e6e23 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -12,8 +12,7 @@ from xclim.core.calendar import compare_offsets, resample_doy, select_time from xclim.core.formatting import prefix_attrs, unprefix_attrs, update_history -from xclim.core.units import convert_units_to -from xclim.core.utils import Quantified, uses_dask +from xclim.core.utils import DateStr, Quantified, uses_dask from . import generic @@ -729,9 +728,6 @@ def standardized_index_fit_params( method : {'ML', 'APP', 'PWM'} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that doesn't involve any optimization. - offset: Quantified - Distributions bounded by zero (e.g. "gamma", "fisk") can be used for datasets with negative values - by using an offset: `da + offset`. \*\*indexer Indexing parameters to compute the indicator on a temporal subset of the data. It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`. @@ -755,11 +751,6 @@ def standardized_index_fit_params( raise NotImplementedError( f"The method `{method}` is not supported for distribution `{dist}`." ) - - if offset is not None: - with xr.set_options(keep_attrs=True): - da = da + convert_units_to(offset, da, context="hydro") - da, group = preprocess_standardized_index(da, freq, window, **indexer) params = da.groupby(group).map(fit, dist=dist, method=method) cal_range = ( @@ -774,7 +765,6 @@ def standardized_index_fit_params( "method": method, "group": group, "units": "", - "offset": offset or "", } method, args = ("", []) if indexer == {} else indexer.popitem() params.attrs["time_indexer"] = (method, *args) @@ -782,24 +772,96 @@ def standardized_index_fit_params( return params -def standardized_index(da: xr.DataArray, params: xr.DataArray): - """Compute standardized index for given fit parameters. - - This computes standardized indices which measure the deviation of variables in the dataset compared - to a reference distribution. The reference is a statistical distribution computed with fitting parameters `params` - over a given calibration period of the dataset. Those fitting parameters are obtained with - ``xclim.standardized_index_fit_params``. +def standardized_index( + da: xr.DataArray, + freq: str | None, + window: int, + dist: str, + method: str, + cal_start: DateStr | None, + cal_end: DateStr | None, + params: Quantified | None, + **indexer, +) -> xr.DataArray: + r"""Standardized Index (SI). Parameters ---------- da : xarray.DataArray - Input array. - Resampling and rolling operations should have already been applied using - ``xclim.indices.preprocess_standardized_index``. + Daily input data. + freq : str, optional + Resampling frequency. A monthly or daily frequency is expected. Option `None` assumes that desired resampling + has already been applied input dataset and will skip the resampling step. + window : int + Averaging window length relative to the resampling frequency. For example, if `freq="MS"`, + i.e. a monthly resampling, the window is an integer number of months. + dist : {"gamma", "fisk"} + Name of the univariate distribution. (see :py:mod:`scipy.stats`). + method : {'APP', 'ML'} + Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method + uses a deterministic function that doesn't involve any optimization. + cal_start : DateStr, optional + Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. + Default option `None` means that the calibration period begins at the start of the input dataset. + cal_end : DateStr, optional + End date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. + Default option `None` means that the calibration period finishes at the end of the input dataset. params : xarray.DataArray - Fit parameters computed using ``xclim.indices.stats.standardized_index_fit_params``. + Fit parameters. + The `params` can be computed using ``xclim.indices.stats.standardized_index_fit_params`` in advance. + The output can be given here as input, and it overrides other options. + \*\*indexer + Indexing parameters to compute the indicator on a temporal subset of the data. + It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`. + + Returns + ------- + xarray.DataArray, [unitless] + Standardized Precipitation Index. + + Notes + ----- + * The length `N` of the N-month SPI is determined by choosing the `window = N`. + * Supported statistical distributions are: ["gamma", "fisk"], where "fisk" is scipy's implementation of + a log-logistic distribution + * If `params` is given as input, it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. + * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in + the inversion to the normal distribution. + + References + ---------- + :cite:cts:`mckee_relationship_1993` """ + if params is not None: + freq, window, indexer = ( + params.attrs[s] for s in ["freq", "window", "time_indexer"] + ) + # Unpack attrs to None and {} if needed + freq = None if freq == "" else freq + indexer = {} if indexer[0] == "" else {indexer[0]: indexer[1:]} + if cal_start or cal_end: + warnings.warn( + "Expected either `cal_{start|end}` or `params`, got both. The `params` input overrides other inputs." + "If `cal_start`, `cal_end`, `freq`, `window`, and/or `dist` were given as input, they will be ignored." + ) + + da, _ = preprocess_standardized_index(da, freq=freq, window=window, **indexer) + if params is None: + params = standardized_index_fit_params( + da.sel(time=slice(cal_start, cal_end)), + freq=None, + window=1, + dist=dist, + method=method, + ) + + # If params only contains a subset of main dataset time grouping + # (e.g. 8/12 months, etc.), it needs to be broadcasted group = params.attrs["group"] + template = da.groupby(group).first() + paramsd = {k: v for k, v in params.sizes.items() if k != "dparams"} + if paramsd != template.sizes: + params = params.broadcast_like(template) def reindex_time(da, da_ref): if group == "time.dayofyear": @@ -809,6 +871,7 @@ def reindex_time(da, da_ref): da["time"] = da_ref.time return da if not uses_dask(da) else da.chunk({"time": -1}) + # this should be restricted to some distributions / in some context probs_of_zero = da.groupby(group).map( lambda x: (x == 0).sum("time") / x.notnull().sum("time") ) @@ -821,7 +884,7 @@ def reindex_time(da, da_ref): coords=dict(dparams=(["loc", "scale"])), attrs=dict(scipy_dist="norm"), ) - std_index = dist_method("ppf", params_norm, probs) + si = dist_method("ppf", params_norm, probs) # A cdf value of 0 or 1 gives ±np.inf when inverted to the normal distribution. # The neighbouring values 0.00...1 and 0.99...9 with a float64 give a standardized index of ± 8.21 # We use this index as reference for maximal/minimal standardized values. @@ -829,5 +892,10 @@ def reindex_time(da, da_ref): # notes that precipitation events (over a month/season) generally occur less than 100 times in the US. # Values of standardized index outside the range [-3.09, 3.09] correspond to probabilities smaller than 0.001 # or greater than 0.999, which could be considered non-statistically relevant for this sample size. - std_index = std_index.clip(-8.21, 8.21) - return std_index + si = si.clip(-8.21, 8.21) + + si.attrs = params.attrs + si.attrs["freq"] = (freq or xr.infer_freq(si.time)) or "undefined" + si.attrs["window"] = window + si.attrs["units"] = "" + return si From 2294163af9c8af5f4be518d2fa2f372773887834 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 16 Feb 2024 22:42:27 -0500 Subject: [PATCH 03/13] add context="hydro" --- tests/test_indices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index e13f2ac96..931fe6f53 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -696,7 +696,7 @@ def test_standardized_precipitation_evapotranspiration_index( tasmin = tasmax - 5 wb = xci.water_budget(pr, None, tasmin, tasmax, tas) # offset like in climate indices library - wb = wb + convert_units_to("1 mm/d", wb) + wb = wb + convert_units_to("1 mm/d", wb, context="hydro") params = xci.stats.standardized_index_fit_params( wb.sel(time=slice("1950", "1980")), From 36d3ce7275f54db1104e2e9fe3867d8e0cb32618 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 19 Feb 2024 09:25:28 -0500 Subject: [PATCH 04/13] update CHANGES --- CHANGES.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 670bebdd4..aeb136205 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -22,6 +22,7 @@ New features and enhancements * New ``xclim.core.calendar.stack_periods`` and ``unstack_periods`` for performing ``rolling(time=...).construct(..., stride=...)`` but with non-uniform temporal periods like years or months. They replace ``xclim.sdba.processing.construct_moving_yearly_window`` and ``unpack_moving_yearly_window`` which are deprecated and will be removed in a future release. * New ``as_dataset`` options for ``xclim.set_options``. When True, indicators will output Datasets instead of DataArrays. (:issue:`1257`, :pull:`1625`). * Added new option for UTCI calculation to cap low wind velocities to a minimum of 0.5 m/s following Bröde (2012) guidelines. (:issue:`1634`, :pull:`1635`). +* Distribution with negative values are directly fitted without need for an offset for distributions such as `gamma, fisk` in ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). Breaking changes ^^^^^^^^^^^^^^^^ @@ -39,6 +40,7 @@ Breaking changes * `black` formatting style has been updated to the 2024 stable conventions. `isort` has been added to the `dev` installation recipe. (:pull:`1626`). * The indice and indicator for ``winter_storm`` has been removed (deprecated since `xclim` v0.46.0 in favour of ``snd_storm_days``). (:pull:`1565`). * `xclim` has dropped support for `scipy` version below v1.9.0 and `numpy` versions below v1.20.0. (:pull:`1565`). +* Argument `offset` is removed from ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). Bug fixes ^^^^^^^^^ From 6a1809c1b6e13c933ab499203de5c7b5e9d3243a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 19 Feb 2024 09:43:06 -0500 Subject: [PATCH 05/13] Keep code compatible with offset, add deprecation warnings --- xclim/indices/_agro.py | 33 +++++++++++++++++++++++++++++++++ xclim/indices/stats.py | 9 +++++++++ 2 files changed, 42 insertions(+) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 0d1fc2d22..bbf3174cf 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,8 @@ # noqa: D100 from __future__ import annotations +import warnings + import numpy as np import xarray @@ -1194,6 +1196,7 @@ def standardized_precipitation_index( @declare_units( wb="[precipitation]", + offset="[precipitation]", params="[]", ) def standardized_precipitation_evapotranspiration_index( @@ -1202,6 +1205,7 @@ def standardized_precipitation_evapotranspiration_index( window: int = 1, dist: str = "gamma", method: str = "APP", + offset: Quantified = "1.000 mm/d", cal_start: DateStr | None = None, cal_end: DateStr | None = None, params: Quantified | None = None, @@ -1230,6 +1234,10 @@ def standardized_precipitation_evapotranspiration_index( `PWM` (probability weighted moments). The approximate method uses a deterministic function that doesn't involve any optimization. Available methods vary with the distribution: 'gamma':{'APP', 'ML', 'PWM'}, 'fisk':{'APP', 'ML'} + offset : Quantified + For distributions bounded by zero (e.g. "gamma", "fisk"), the two-parameters distributions only accept positive + values. An offset can be added to make sure this is the case. This option will be removed in xclim >=0.49.0, ``xclim`` + will rely on three-parameters distributions instead. cal_start : DateStr, optional Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. Default option `None` means that the calibration period begins at the start of the input dataset. @@ -1259,6 +1267,31 @@ def standardized_precipitation_evapotranspiration_index( See Standardized Precipitation Index (SPI) for more details on usage. """ + uses_default_offset = offset != "1.000 mm/d" + if uses_default_offset is False: + warnings.warn("Inputting an offset will be deprecated in xclim>=0.49.0. ") + if params is not None: + params_offset = params.attrs["offset"] + if uses_default_offset is False and offset != params_offset: + warnings.warn( + "The offset in `params` differs from the input `offset`." + "Proceeding with the value given in `params`." + ) + offset = params_offset + offset = 0 if offset == "" else convert_units_to(offset, wb, context="hydro") + # Allowed distributions are constrained by the SPI function + if dist in ["gamma", "fisk"] and offset <= 0: + raise ValueError( + "The water budget must be shifted towards positive values to be used with `gamma` and `fisk` " + "distributions which are bounded by zero (only when using two-parameters distributions: in xclim>=0.49.0," + "three-parameters distributions are used to accommodate negative values). Only positive offsets are accepted." + ) + # Note that the default behaviour would imply an offset for any distribution, even those distributions + # that can accommodate negative values of the water budget. This needs to be changed in future versions + # of the index. + if offset != 0: + with xarray.set_options(keep_attrs=True): + wb = wb + offset dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist in dist_methods.keys(): if method not in dist_methods[dist]: diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 1f06e6e23..2efe77736 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -12,6 +12,7 @@ from xclim.core.calendar import compare_offsets, resample_doy, select_time from xclim.core.formatting import prefix_attrs, unprefix_attrs, update_history +from xclim.core.units import convert_units_to from xclim.core.utils import DateStr, Quantified, uses_dask from . import generic @@ -743,6 +744,12 @@ def standardized_index_fit_params( * Gamma ("gamma") : "ML", "APP", "PWM" * Log-logistic ("fisk") : "ML", "APP" """ + if offset: + warnings.warn("Inputting an offset will be deprecated in xclim>=0.49.0. ") + if offset is not None: + with xr.set_options(keep_attrs=True): + da = da + convert_units_to(offset, da, context="hydro") + # "WPM" method doesn't seem to work for gamma or pearson3 dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} if dist not in dist_and_methods: @@ -769,6 +776,8 @@ def standardized_index_fit_params( method, args = ("", []) if indexer == {} else indexer.popitem() params.attrs["time_indexer"] = (method, *args) + params.attrs["offset"] = offset or "" + return params From a442fcc60643170c1777737ca4c3d7eaab6c280a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 19 Feb 2024 09:50:38 -0500 Subject: [PATCH 06/13] offset=0 by default -> produce behaviour of future xclim versions --- xclim/indices/_agro.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index bbf3174cf..70706f5f1 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1205,7 +1205,7 @@ def standardized_precipitation_evapotranspiration_index( window: int = 1, dist: str = "gamma", method: str = "APP", - offset: Quantified = "1.000 mm/d", + offset: Quantified = "", cal_start: DateStr | None = None, cal_end: DateStr | None = None, params: Quantified | None = None, @@ -1267,7 +1267,7 @@ def standardized_precipitation_evapotranspiration_index( See Standardized Precipitation Index (SPI) for more details on usage. """ - uses_default_offset = offset != "1.000 mm/d" + uses_default_offset = offset != "" if uses_default_offset is False: warnings.warn("Inputting an offset will be deprecated in xclim>=0.49.0. ") if params is not None: @@ -1286,6 +1286,7 @@ def standardized_precipitation_evapotranspiration_index( "distributions which are bounded by zero (only when using two-parameters distributions: in xclim>=0.49.0," "three-parameters distributions are used to accommodate negative values). Only positive offsets are accepted." ) + # Note that the default behaviour would imply an offset for any distribution, even those distributions # that can accommodate negative values of the water budget. This needs to be changed in future versions # of the index. From 915e9553ef54311e594a6ad86a1df0de6711b894 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Mon, 19 Feb 2024 09:53:43 -0500 Subject: [PATCH 07/13] MOre details in CHANGES, breaking --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index aeb136205..09d9617e8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -40,7 +40,7 @@ Breaking changes * `black` formatting style has been updated to the 2024 stable conventions. `isort` has been added to the `dev` installation recipe. (:pull:`1626`). * The indice and indicator for ``winter_storm`` has been removed (deprecated since `xclim` v0.46.0 in favour of ``snd_storm_days``). (:pull:`1565`). * `xclim` has dropped support for `scipy` version below v1.9.0 and `numpy` versions below v1.20.0. (:pull:`1565`). -* Argument `offset` is removed from ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). +* Argument `offset` is removed from ``xclim.indices.standardized_precipitation_evapotranspiration_index``. This will slightly change results when using SPEI with default values (:issue:`1477` :pull:`1653`). Bug fixes ^^^^^^^^^ From 2c4a684f7f598a6cd1116994bed33edef9540261 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 22 Feb 2024 17:05:29 -0500 Subject: [PATCH 08/13] update version in CHANGES --- CHANGES.rst | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index e14bcce9a..a6b93426d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,6 +2,19 @@ Changelog ========= +v0.48.2 +-------------------- +Contributors to this version: Éric Dupuis (:user:`coxipi`) + +New features and enhancements +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +* Distribution with negative values are directly fitted without need for an offset for distributions such as `gamma, fisk` in ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). + +Breaking changes +^^^^^^^^^^^^^^^^ +* Argument `offset` is removed from ``xclim.indices.standardized_precipitation_evapotranspiration_index``. This will slightly change results when using SPEI with default values (:issue:`1477` :pull:`1653`). + + v0.48.1 (2024-02-20) -------------------- Contributors to this version: Trevor James Smith (:user:`Zeitsperre`). From be11a475015c1f51d3ef8e5535069cf47e859867 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 22 Feb 2024 17:06:35 -0500 Subject: [PATCH 09/13] remove entries in old versions, CHANGES --- CHANGES.rst | 2 -- 1 file changed, 2 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index a6b93426d..b925b8f62 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -50,7 +50,6 @@ New features and enhancements * Distribution instances can now be passed to the ``dist`` argument of most statistical indices. (:pull:`1644`). * Added a new ``xclim.indices.generic.select_rolling_resample_op`` function to allow for computing rolling statistics. (:issue:`1480`, :pull:`1643`). * Add the possibility to use a group with a window in ``xc.sdba.processing.reordering``. (:pull:`1566`). -* Distribution with negative values are directly fitted without need for an offset for distributions such as `gamma, fisk` in ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). Breaking changes ^^^^^^^^^^^^^^^^ @@ -70,7 +69,6 @@ Breaking changes * `xclim` has dropped support for `scipy` versions below v1.9.0 and `numpy` versions below v1.20.0. (:pull:`1565`). * For generic function ``select_resample_op`` and ``core.units.to_agg_units``, operation "sum" will now return the same units as the input, and not implicitly be translated to an "integral". (:issue:`1645`, :pull:`1649`). * `lmoments3` was removed as a dependency of `xclim` due to incompatible licensing (GPLv3 vs `xclim`'s Apache 2.0). Depending on the outcome of efforts to modify the licensing of `lmoments3`, this change may eventually be reverted. See `Ouranosinc/lmoments3#12 `_. See also the "frequency analysis" notebook for an example on how to continue using the probability weighted moments method for fitting distributions. (:issue:`1620`, :pull:`1644`). -* Argument `offset` is removed from ``xclim.indices.standardized_precipitation_evapotranspiration_index``. This will slightly change results when using SPEI with default values (:issue:`1477` :pull:`1653`). Bug fixes ^^^^^^^^^ From b989dbff00e245aa0424ccdc22419ad407acdb20 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Wed, 20 Mar 2024 16:12:12 -0400 Subject: [PATCH 10/13] better estimators using method of moments for 2-p fisk --- xclim/indices/stats.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 2c8ebbabc..336a292a4 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -511,17 +511,17 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: x_pos = x - (xmin if xmin <= 0 else 0) x_pos = x_pos[x_pos > 0] m = x_pos.mean() - v = x_pos.var() - # pdf = (beta/alpha) (x/alpha)^{beta-1}/ (1+(x/alpha)^{beta})^2 - # Compute f_1 and f_2 which only depend on beta: - # f_1 := mean/alpha = /alpha - # f_2 := variance/alpha^2 = ( - ^2)/alpha^2 - # In the large beta limit, f_1 -> 1 and f_1/sqrt(f_2) -> 0.56*beta - 0.25 - # Solve for alpha and beta below: - alpha, beta = m, (1 / 0.56) * (m / np.sqrt(v) + 1 / 4) + m2 = (x_pos**2).mean() + # method of moments: + # LHS is computed analytically with the two-parameters log-logistic distribution + # and depends on alpha,beta + # RHS is from the sample + # = m + # / ^2 = m2/m**2 + # solving these equations yields + alpha = 2 * m**3 / (m2 + m**2) + beta = np.pi * m / np.sqrt(3) / np.sqrt(m2 - m**2) kwargs = {"scale": alpha} - if xmin < 0: - kwargs["loc"] = xmin return (beta,), kwargs return (), {} From 52bff862c3237270b97f6f8aa9f3496c37c4b2b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 22 Mar 2024 17:50:28 -0400 Subject: [PATCH 11/13] add `loc` estimates in `_fit_start` for 'gamma' and 'fisk' --- xclim/indices/_agro.py | 10 ++++-- xclim/indices/stats.py | 73 ++++++++++++++++++++++++++++++------------ 2 files changed, 61 insertions(+), 22 deletions(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 70706f5f1..5653aaabd 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1098,6 +1098,7 @@ def standardized_precipitation_index( window: int = 1, dist: str = "gamma", method: str = "APP", + fitkwargs: dict = {}, cal_start: DateStr | None = None, cal_end: DateStr | None = None, params: Quantified | None = None, @@ -1120,6 +1121,8 @@ def standardized_precipitation_index( method : {'APP', 'ML'} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that doesn't involve any optimization. + fitkwargs : dict + Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). cal_start : DateStr, optional Start date of the calibration period. A `DateStr` is expected, that is a `str` in format `"YYYY-MM-DD"`. Default option `None` means that the calibration period begins at the start of the input dataset. @@ -1189,7 +1192,7 @@ def standardized_precipitation_index( else: raise NotImplementedError(f"{dist} distribution is not implemented yet") spi = standardized_index( - pr, freq, window, dist, method, cal_start, cal_end, params, **indexer + pr, freq, window, dist, method, fitkwargs, cal_start, cal_end, params, **indexer ) return spi @@ -1205,6 +1208,7 @@ def standardized_precipitation_evapotranspiration_index( window: int = 1, dist: str = "gamma", method: str = "APP", + fitkwargs: dict = {}, offset: Quantified = "", cal_start: DateStr | None = None, cal_end: DateStr | None = None, @@ -1234,6 +1238,8 @@ def standardized_precipitation_evapotranspiration_index( `PWM` (probability weighted moments). The approximate method uses a deterministic function that doesn't involve any optimization. Available methods vary with the distribution: 'gamma':{'APP', 'ML', 'PWM'}, 'fisk':{'APP', 'ML'} + fitkwargs : dict + Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). offset : Quantified For distributions bounded by zero (e.g. "gamma", "fisk"), the two-parameters distributions only accept positive values. An offset can be added to make sure this is the case. This option will be removed in xclim >=0.49.0, ``xclim`` @@ -1302,7 +1308,7 @@ def standardized_precipitation_evapotranspiration_index( else: raise NotImplementedError(f"{dist} distribution is not implemented yet") spei = standardized_index( - wb, freq, window, dist, method, cal_start, cal_end, params, **indexer + wb, freq, window, dist, method, fitkwargs, cal_start, cal_end, params, **indexer ) return spei diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 336a292a4..892aa603b 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -490,26 +490,36 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: return (chat,), {"loc": loc, "scale": scale} if dist in ["gamma"]: - # not sure if the approximation holds for xmin < 0 - xmin = x.min() - x_pos = x - (xmin if xmin <= 0 else 0) - x_pos = x_pos[x_pos > 0] + if "floc" in fitkwargs: + loc0 = fitkwargs["floc"] + else: + xs = sorted(x) + x1, x2, xn = xs[0], xs[1], xs[-1] + n = len(x) + cv = x.std() / x.mean() + p = 100 * (0.48265 + 0.32967 * cv) * n ** (-0.2984 * cv) + xp = np.percentile(x, p) + loc0 = (x1 * xn - xp**2) / (x1 + xn - 2 * xp) + loc0 = loc0 if loc0 < x1 else 0.9999 * x1 + x_pos = x - loc0 m = x_pos.mean() log_of_mean = np.log(m) mean_of_logs = np.log(x_pos).mean() - a = log_of_mean - mean_of_logs - alpha = (1 + np.sqrt(1 + 4 * a / 3)) / (4 * a) - beta = m / alpha - kwargs = {"scale": beta} - if xmin < 0: - kwargs["loc"] = xmin - return (alpha,), kwargs + A = log_of_mean - mean_of_logs + a0 = (1 + np.sqrt(1 + 4 * A / 3)) / (4 * A) + scale0 = m / a0 + kwargs = {"scale": scale0, "loc": loc0} + return (a0,), kwargs if dist in ["fisk"]: - # not sure if the approximation holds for xmin < 0 - xmin = x.min() - x_pos = x - (xmin if xmin <= 0 else 0) - x_pos = x_pos[x_pos > 0] + if "floc" in fitkwargs: + loc0 = fitkwargs["floc"] + else: + xs = sorted(x) + x1, x2, xn = xs[0], xs[1], xs[-1] + loc0 = (x1 * xn - x2**2) / (x1 + xn - 2 * x2) + loc0 = loc0 if loc0 < x1 else 0.9999 * x1 + x_pos = x - loc0 m = x_pos.mean() m2 = (x_pos**2).mean() # method of moments: @@ -519,10 +529,10 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: # = m # / ^2 = m2/m**2 # solving these equations yields - alpha = 2 * m**3 / (m2 + m**2) - beta = np.pi * m / np.sqrt(3) / np.sqrt(m2 - m**2) - kwargs = {"scale": alpha} - return (beta,), kwargs + scale0 = 2 * m**3 / (m2 + m**2) + c0 = np.pi * m / np.sqrt(3) / np.sqrt(m2 - m**2) + kwargs = {"scale": scale0, "loc": loc0} + return (c0,), kwargs return (), {} @@ -673,6 +683,7 @@ def standardized_index_fit_params( window: int, dist: str | scipy.stats.rv_continuous, method: str, + fitkwargs: dict = {}, offset: Quantified | None = None, **indexer, ) -> xr.DataArray: @@ -698,6 +709,8 @@ def standardized_index_fit_params( method : {'ML', 'APP', 'PWM'} Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method uses a deterministic function that doesn't involve any optimization. + fitkwargs : dict + Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). \*\*indexer Indexing parameters to compute the indicator on a temporal subset of the data. It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`. @@ -729,7 +742,24 @@ def standardized_index_fit_params( f"The method `{method}` is not supported for distribution `{dist.name}`." ) da, group = preprocess_standardized_index(da, freq, window, **indexer) - params = da.groupby(group).map(fit, dist=dist, method=method) + + # convert floc units if needed + # this should be passed to scipy eventually, so it should not support strings + # with units. Perhaps using `fitkwargs` in this context is misleading? + for fpar in ["floc", "fscale"]: + if fpar in fitkwargs.keys(): + if np.isscalar(fitkwargs[fpar]) is False: + fitkwargs[fpar] = convert_units_to(fitkwargs[fpar], da, context="hydro") + + # Use zero inflated distributions + # Idea: Perhaps having specific zero-inflated distributions or an option passed + # to the standardized index (SPI/SPEI) to state we want to use zero-inflated + # would be a better way to organize this + if da.min() == 0 and dist.name in ["gamma", "fisk"]: + da0 = da.where(da > 0).copy() + else: + da0 = da + params = da0.groupby(group).map(fit, dist=dist, method=method, **fitkwargs) cal_range = ( da.time.min().dt.strftime("%Y-%m-%d").item(), da.time.max().dt.strftime("%Y-%m-%d").item(), @@ -757,6 +787,7 @@ def standardized_index( window: int, dist: str | scipy.stats.rv_continuous | None, method: str, + fitkwargs: dict, cal_start: DateStr | None, cal_end: DateStr | None, params: Quantified | None, @@ -837,6 +868,7 @@ def standardized_index( window=1, dist=dist, method=method, + fitkwargs=fitkwargs, ) # If params only contains a subset of main dataset time grouping @@ -860,6 +892,7 @@ def reindex_time(da, da_ref): lambda x: (x == 0).sum("time") / x.notnull().sum("time") ) params, probs_of_zero = (reindex_time(dax, da) for dax in [params, probs_of_zero]) + # should `da` below exclude zeros? dist_probs = dist_method("cdf", params, da, dist=dist) probs = probs_of_zero + ((1 - probs_of_zero) * dist_probs) params_norm = xr.DataArray( From b20d7a9477eb32aa293b469fc830ced7536ee137 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 22 Mar 2024 18:17:11 -0400 Subject: [PATCH 12/13] APP method can only be used with floc=0 --- xclim/indices/stats.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 892aa603b..b9a7f3123 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -725,6 +725,8 @@ def standardized_index_fit_params( Supported combinations of `dist` and `method` are: * Gamma ("gamma") : "ML", "APP", "PWM" * Log-logistic ("fisk") : "ML", "APP" + + "APP" method only supports two-parameter distributions. Parameter `loc` will be set to 0 (setting `floc=0` in `fitkwargs`). """ if offset: warnings.warn("Inputting an offset will be deprecated in xclim>=0.49.0. ") @@ -741,6 +743,19 @@ def standardized_index_fit_params( raise NotImplementedError( f"The method `{method}` is not supported for distribution `{dist.name}`." ) + if method == "APP": + if "floc" in fitkwargs.keys(): + if fitkwargs["floc"] != 0: + raise ValueError( + "The APP method is only supported for two-parameter distributions with `gamma` or `fisk`. `loc` parameter must be set 0." + "Pass `floc=0` in `fitkwargs`." + ) + else: + fitkwargs["floc"] = 0 + warnings.warn( + "The APP method is only supported for two-parameter distributions with `gamma` or `fisk`." + "Location parameter `loc` will be set to 0. To avoid this warning, set `floc=0` manually in `fitkwargs`." + ) da, group = preprocess_standardized_index(da, freq, window, **indexer) # convert floc units if needed @@ -775,7 +790,6 @@ def standardized_index_fit_params( } method, args = ("", []) if indexer == {} else indexer.popitem() params.attrs["time_indexer"] = (method, *args) - params.attrs["offset"] = offset or "" return params From 30ab01835b5137abc707f982a5e180edf4643786 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Tue, 16 Apr 2024 15:57:13 -0400 Subject: [PATCH 13/13] Update CHANGES.rst missing 's' Co-authored-by: Vasco Schiavo <115561717+VascoSch92@users.noreply.github.com> --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index ecc0072f9..6d56a2b1c 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,7 +12,7 @@ Announcements New features and enhancements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -* Distribution with negative values are directly fitted without need for an offset for distributions such as `gamma, fisk` in ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). +* Distributions with negative values are directly fitted without need for an offset for distributions such as `gamma, fisk` in ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1477` :pull:`1653`). Breaking changes ^^^^^^^^^^^^^^^^