From 68267b01a4fe63038d8e11e146d802b48d96ab14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 17:20:22 +0200 Subject: [PATCH 01/49] transform: add field, store, process kw to routines; respect normalizer; add wrapper; better checks --- docs/source/transform.rst | 2 - gstools/transform/__init__.py | 48 +- gstools/transform/field.py | 834 ++++++++++++++++++++++++++-------- 3 files changed, 696 insertions(+), 188 deletions(-) diff --git a/docs/source/transform.rst b/docs/source/transform.rst index 1d41d590a..a456cb121 100644 --- a/docs/source/transform.rst +++ b/docs/source/transform.rst @@ -2,8 +2,6 @@ gstools.transform ================= .. automodule:: gstools.transform - :members: - :undoc-members: .. raw:: latex diff --git a/gstools/transform/__init__.py b/gstools/transform/__init__.py index 0c52fd6b2..3a2a50061 100644 --- a/gstools/transform/__init__.py +++ b/gstools/transform/__init__.py @@ -4,10 +4,20 @@ .. currentmodule:: gstools.transform -Field-Transformations +Wrapper +^^^^^^^ + +.. autosummary:: + :toctree: generated + + apply + +Field Transformations ^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + binary discrete boxcox @@ -17,11 +27,29 @@ normal_to_uniform normal_to_arcsin normal_to_uquad + apply_function + +Array Transformations +^^^^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: generated + + array_discrete + array_boxcox + array_zinnharvey + array_force_moments + array_to_lognormal + array_to_uniform + array_to_arcsin + array_to_uquad ---- """ from gstools.transform.field import ( + apply, + apply_function, binary, discrete, boxcox, @@ -31,9 +59,19 @@ normal_to_uniform, normal_to_arcsin, normal_to_uquad, + array_discrete, + array_boxcox, + array_zinnharvey, + array_force_moments, + array_to_lognormal, + array_to_uniform, + array_to_arcsin, + array_to_uquad, ) __all__ = [ + "apply", + "apply_function", "binary", "discrete", "boxcox", @@ -43,4 +81,12 @@ "normal_to_uniform", "normal_to_arcsin", "normal_to_uquad", + "array_discrete", + "array_boxcox", + "array_zinnharvey", + "array_force_moments", + "array_to_lognormal", + "array_to_uniform", + "array_to_arcsin", + "array_to_uquad", ] diff --git a/gstools/transform/field.py b/gstools/transform/field.py index bc0e201ee..efd7387c9 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -6,7 +6,17 @@ The following functions are provided +Wrapper +^^^^^^^ + .. autosummary:: + apply + +Transformations +^^^^^^^^^^^^^^^ + +.. autosummary:: + apply_function binary discrete boxcox @@ -16,14 +26,33 @@ normal_to_uniform normal_to_arcsin normal_to_uquad + +Low-Level Routines +^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + array_discrete + array_boxcox + array_zinnharvey + array_force_moments + array_to_lognormal + array_to_uniform + array_to_arcsin + array_to_uquad """ -# pylint: disable=C0103 +# pylint: disable=C0103, C0123, R0911 from warnings import warn import numpy as np from scipy.special import erf, erfinv - +from gstools.normalizer import ( + Normalizer, + remove_trend_norm_mean, + apply_mean_norm_trend, +) __all__ = [ + "apply", + "apply_function", "binary", "discrete", "boxcox", @@ -33,10 +62,203 @@ "normal_to_uniform", "normal_to_arcsin", "normal_to_uquad", + "array_discrete", + "array_boxcox", + "array_zinnharvey", + "array_force_moments", + "array_to_lognormal", + "array_to_uniform", + "array_to_arcsin", + "array_to_uquad", ] -def binary(fld, divide=None, upper=None, lower=None): +def _get_field_io(field, store): + if not isinstance(field, str): + raise ValueError( + f"transform: given field name '{field}' is not a string." + ) + if isinstance(store, str): + if not store.isidentifier(): + raise ValueError( + f"transform: given store name '{store}' is not valid" + ) + save = True + output_field = store + else: + save = bool(store) + output_field = field + return output_field, save + + +def _check_input_field(fld, field): + if not hasattr(fld, field) or getattr(fld, field) is None: + raise ValueError(f"transform: selected field '{field}' not present") + return getattr(fld, field) + + +def _pre_process(fld, data, keep_mean): + return remove_trend_norm_mean( + pos=fld.pos, + field=data, + mean=None if keep_mean else fld.mean, + normalizer=fld.normalizer, + trend=fld.trend, + mesh_type=fld.mesh_type, + value_type=fld.value_type, + check_shape=False, + ) + + +def _post_process(fld, data, keep_mean): + return apply_mean_norm_trend( + pos=fld.pos, + field=data, + mean=None if keep_mean else fld.mean, + normalizer=fld.normalizer, + trend=fld.trend, + mesh_type=fld.mesh_type, + value_type=fld.value_type, + check_shape=False, + ) + + +def _check_for_default_normal(fld): + if not type(fld.normalizer) == Normalizer: + raise ValueError( + "transform: need a normal field but there is a normalizer defined" + ) + if fld.trend is not None: + raise ValueError( + "transform: need a normal field but there is a trend defined" + ) + if callable(fld.mean) or fld.mean is None: + raise ValueError( + "transform: need a normal field but mean is not constant" + ) + + +def apply(fld, method, field="field", store=True, process=False, **kwargs): + """ + Apply field transformation. + + Parameters + ---------- + fld : :any:`Field` + Field class containing a generated field. + method : :class:`str` + Method to use. See :any:`transform` for available transformations. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + **kwargs + Keyword arguments forwarded to selected method. + + Raises + ------ + ValueError + When method is unknown. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + kwargs["field"] = field + kwargs["store"] = store + kwargs["process"] = process + method = str(method) # ensure method is a string + if method == "binary": + return binary(fld, **kwargs) + if method == "discrete": + return discrete(fld, **kwargs) + if method == "boxcox": + return boxcox(fld, **kwargs) + if method == "zinnharvey": + return zinnharvey(fld, **kwargs) + if method.endswith("force_moments"): + return normal_force_moments(fld, **kwargs) + if method.endswith("lognormal"): + return normal_to_lognormal(fld, **kwargs) + if method.endswith("uniform"): + return normal_to_uniform(fld, **kwargs) + if method.endswith("arcsin"): + return normal_to_arcsin(fld, **kwargs) + if method.endswith("uquad"): + return normal_to_uquad(fld, **kwargs) + if method.endswith("function"): + return apply_function(fld, **kwargs) + raise ValueError(f"transform.apply: unknown method '{method}'.") + + +def apply_function( + fld, + function, + field="field", + store=True, + process=False, + keep_mean=False, + **kwargs, +): + """ + Apply function as field transformation. + + Parameters + ---------- + fld : :any:`Field` + Field class containing a generated field. + function : :any:`callable` + Function to use. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is False. + **kwargs + Keyword arguments forwarded to given function. + + Raises + ------ + ValueError + When function is not callable. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not callable(function): + raise ValueError("apply_function: function not callable.") + output_field, save = _get_field_io(field, store) + data = _check_input_field(fld, field) + if process: + data = _pre_process(fld, data, keep_mean=keep_mean) + data = function(data, **kwargs) + if process: + data = _post_process(fld, data, keep_mean=keep_mean) + return fld.post_field(data, name=output_field, process=False, save=save) + + +def binary( + fld, + divide=None, + upper=None, + lower=None, + field="field", + store=True, + process=False, +): """ Binary transformation. @@ -45,8 +267,7 @@ def binary(fld, divide=None, upper=None, lower=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. divide : :class:`float`, optional The dividing value. Default: ``fld.mean`` @@ -56,17 +277,50 @@ def binary(fld, divide=None, upper=None, lower=None): lower : :class:`float`, optional The resulting lower value of the field. Default: ``mean - sqrt(fld.model.sill)`` - """ - if fld.field is None: - print("binary: no field stored in SRF class.") - else: - divide = fld.mean if divide is None else divide - upper = fld.mean + np.sqrt(fld.model.sill) if upper is None else upper - lower = fld.mean - np.sqrt(fld.model.sill) if lower is None else lower - discrete(fld, [lower, upper], thresholds=[divide]) + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. - -def discrete(fld, values, thresholds="arithmetic"): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process and divide is None: + _check_for_default_normal(fld) + if divide is None: + mean = 0.0 if process else fld.mean + divide = mean if divide is None else divide + upper = mean + np.sqrt(fld.model.sill) if upper is None else upper + lower = mean - np.sqrt(fld.model.sill) if lower is None else lower + kw = dict( + values=[lower, upper], + thresholds=[divide], + ) + return apply_function( + fld=fld, + function=array_discrete, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) + + +def discrete( + fld, + values, + thresholds="arithmetic", + field="field", + store=True, + process=False, +): """ Discrete transformation. @@ -76,8 +330,7 @@ def discrete(fld, values, thresholds="arithmetic"): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. values : :any:`numpy.ndarray` The discrete values the field will take thresholds : :class:`str` or :any:`numpy.ndarray`, optional @@ -87,50 +340,40 @@ def discrete(fld, values, thresholds="arithmetic"): * "equal": devide the field into equal parts * an array of explicitly given thresholds Default: "arithmetic" - """ - if fld.field is None: - print("discrete: no field stored in SRF class.") - else: - if thresholds == "arithmetic": - # just in case, sort the values - values = np.sort(values) - thresholds = (values[1:] + values[:-1]) / 2 - elif thresholds == "equal": - values = np.array(values) - n = len(values) - p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] - rescale = np.sqrt(fld.model.sill * 2) - # use quantile of the normal distribution to get equal ratios - thresholds = fld.mean + rescale * erfinv(2 * p - 1) - else: - if len(values) != len(thresholds) + 1: - raise ValueError( - "discrete transformation: " - "len(values) != len(thresholds) + 1" - ) - values = np.array(values) - thresholds = np.array(thresholds) - # check thresholds - if not np.all(thresholds[:-1] < thresholds[1:]): - raise ValueError( - "discrete transformation: thresholds need to be ascending." - ) - # use a separate result so the intermediate results are not affected - result = np.empty_like(fld.field) - # handle edge cases - result[fld.field <= thresholds[0]] = values[0] - result[fld.field > thresholds[-1]] = values[-1] - for i, value in enumerate(values[1:-1]): - result[ - np.logical_and( - thresholds[i] < fld.field, fld.field <= thresholds[i + 1] - ) - ] = value - # overwrite the field - fld.field = result - - -def boxcox(fld, lmbda=1, shift=0): + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process and thresholds == "equal": + _check_for_default_normal(fld) + kw = dict( + values=values, + thresholds=thresholds, + mean=0.0 if process else fld.mean, + var=fld.model.sill, + ) + return apply_function( + fld=fld, + function=array_discrete, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) + + +def boxcox(fld, lmbda=1, shift=0, field="field", store=True, process=False): """ (Inverse) Box-Cox transformation to denormalize data. @@ -142,8 +385,7 @@ def boxcox(fld, lmbda=1, shift=0): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. lmbda : :class:`float`, optional The lambda parameter of the Box-Cox transformation. For ``lmbda=0`` one obtains the log-normal transformation. @@ -152,20 +394,33 @@ def boxcox(fld, lmbda=1, shift=0): The shift parameter from the two-parametric Box-Cox transformation. The field will be shifted by that value before transformation. Default: ``0`` + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. """ - if fld.field is None: - print("Box-Cox: no field stored in SRF class.") - else: - fld.mean += shift - fld.field += shift - if np.isclose(lmbda, 0): - normal_to_lognormal(fld) - if np.min(fld.field) < -1 / lmbda: - warn("Box-Cox: Some values will be cut off!") - fld.field = (np.maximum(lmbda * fld.field + 1, 0)) ** (1 / lmbda) + kw = dict(lmbda=lmbda, shift=shift) + return apply_function( + fld=fld, + function=array_boxcox, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) -def zinnharvey(fld, conn="high"): +def zinnharvey(fld, conn="high", field="field", store=True, process=False): """ Zinn and Harvey transformation to connect low or high values. @@ -174,19 +429,39 @@ def zinnharvey(fld, conn="high"): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. conn : :class:`str`, optional Desired connectivity. Either "low" or "high". Default: "high" + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. """ - if fld.field is None: - print("zinnharvey: no field stored in SRF class.") - else: - fld.field = _zinnharvey(fld.field, conn, fld.mean, fld.model.sill) + if not process: + _check_for_default_normal(fld) + kw = dict(conn=conn, mean=0.0 if process else fld.mean, var=fld.model.sill) + return apply_function( + fld=fld, + function=array_zinnharvey, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) -def normal_force_moments(fld): +def normal_force_moments(fld, field="field", store=True, process=False): """ Force moments of a normal distributed field. @@ -195,16 +470,36 @@ def normal_force_moments(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. """ - if fld.field is None: - print("normal_force_moments: no field stored in SRF class.") - else: - fld.field = _normal_force_moments(fld.field, fld.mean, fld.model.sill) + if not process: + _check_for_default_normal(fld) + kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill) + return apply_function( + fld=fld, + function=array_force_moments, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) -def normal_to_lognormal(fld): +def normal_to_lognormal(fld, field="field", store=True, process=False): """ Transform normal distribution to log-normal distribution. @@ -213,16 +508,32 @@ def normal_to_lognormal(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. """ - if fld.field is None: - print("normal_to_lognormal: no field stored in SRF class.") - else: - fld.field = _normal_to_lognormal(fld.field) + return apply_function( + fld=fld, + function=array_to_lognormal, + field=field, + store=store, + process=process, + keep_mean=True, # apply to normal field including mean + ) -def normal_to_uniform(fld): +def normal_to_uniform(fld, field="field", store=True, process=False): """ Transform normal distribution to uniform distribution on [0, 1]. @@ -231,16 +542,25 @@ def normal_to_uniform(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. - """ - if fld.field is None: - print("normal_to_uniform: no field stored in SRF class.") - else: - fld.field = _normal_to_uniform(fld.field, fld.mean, fld.model.sill) - - -def normal_to_arcsin(fld, a=None, b=None): + Field class containing a generated field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill) + return apply_function( + fld=fld, + function=array_to_uniform, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) + + +def normal_to_arcsin( + fld, a=None, b=None, field="field", store=True, process=False +): """ Transform normal distribution to the bimodal arcsin distribution. @@ -251,27 +571,44 @@ def normal_to_arcsin(fld, a=None, b=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. a : :class:`float`, optional Parameter a of the arcsin distribution (lower bound). Default: keep mean and variance b : :class:`float`, optional Parameter b of the arcsin distribution (upper bound). Default: keep mean and variance - """ - if fld.field is None: - print("normal_to_arcsin: no field stored in SRF class.") - else: - a = fld.mean - np.sqrt(2.0 * fld.model.sill) if a is None else a - b = fld.mean + np.sqrt(2.0 * fld.model.sill) if b is None else b - fld.field = _normal_to_arcsin( - fld.field, fld.mean, fld.model.sill, a, b - ) - fld.mean = (a + b) / 2.0 + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. - -def normal_to_uquad(fld, a=None, b=None): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill, a=a, b=b) + return apply_function( + fld=fld, + function=array_to_arcsin, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) + + +def normal_to_uquad( + fld, a=None, b=None, field="field", store=True, process=False +): """ Transform normal distribution to U-quadratic distribution. @@ -282,36 +619,154 @@ def normal_to_uquad(fld, a=None, b=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. a : :class:`float`, optional Parameter a of the U-quadratic distribution (lower bound). Default: keep mean and variance b : :class:`float`, optional Parameter b of the U-quadratic distribution (upper bound). Default: keep mean and variance + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill, a=a, b=b) + return apply_function( + fld=fld, + function=array_to_uquad, + field=field, + store=store, + process=process, + keep_mean=False, + **kw, + ) + + +# low level functions + + +def array_discrete( + field, values, thresholds="arithmetic", mean=None, var=None +): """ - if fld.field is None: - print("normal_to_uquad: no field stored in SRF class.") + Discrete transformation. + + After this transformation, the field has only `len(values)` discrete + values. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + values : :any:`numpy.ndarray` + The discrete values the field will take + thresholds : :class:`str` or :any:`numpy.ndarray`, optional + the thresholds, where the value classes are separated + possible values are: + * "arithmetic": the mean of the 2 neighbouring values + * "equal": devide the field into equal parts + * an array of explicitly given thresholds + Default: "arithmetic" + mean : :class:`float`or :any:`None` + Mean of the field for "equal" thresholds. Default: np.mean(field) + var : :class:`float`or :any:`None` + Variance of the field for "equal" thresholds. Default: np.var(field) + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + if thresholds == "arithmetic": + # just in case, sort the values + values = np.sort(values) + thresholds = (values[1:] + values[:-1]) / 2 + elif thresholds == "equal": + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + values = np.array(values) + n = len(values) + p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] + rescale = np.sqrt(var * 2) + # use quantile of the normal distribution to get equal ratios + thresholds = mean + rescale * erfinv(2 * p - 1) else: - a = fld.mean - np.sqrt(5.0 / 3.0 * fld.model.sill) if a is None else a - b = fld.mean + np.sqrt(5.0 / 3.0 * fld.model.sill) if b is None else b - fld.field = _normal_to_uquad(fld.field, fld.mean, fld.model.sill, a, b) - fld.mean = (a + b) / 2.0 + if len(values) != len(thresholds) + 1: + raise ValueError( + "discrete transformation: " + "len(values) != len(thresholds) + 1" + ) + values = np.array(values) + thresholds = np.array(thresholds) + # check thresholds + if not np.all(thresholds[:-1] < thresholds[1:]): + raise ValueError( + "discrete transformation: thresholds need to be ascending." + ) + # use a separate result so the intermediate results are not affected + result = np.empty_like(field) + # handle edge cases + result[field <= thresholds[0]] = values[0] + result[field > thresholds[-1]] = values[-1] + for i, value in enumerate(values[1:-1]): + result[ + np.logical_and(thresholds[i] < field, field <= thresholds[i + 1]) + ] = value + return result -# low level functions +def array_boxcox(field, lmbda=1, shift=0): + """ + (Inverse) Box-Cox transformation to denormalize data. + After this transformation, the again Box-Cox transformed field is normal + distributed. + + See: https://en.wikipedia.org/wiki/Power_transform#Box%E2%80%93Cox_transformation -def _zinnharvey(field, conn="high", mean=None, var=None): + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + lmbda : :class:`float`, optional + The lambda parameter of the Box-Cox transformation. + For ``lmbda=0`` one obtains the log-normal transformation. + Default: ``1`` + shift : :class:`float`, optional + The shift parameter from the two-parametric Box-Cox transformation. + The field will be shifted by that value before transformation. + Default: ``0`` + """ + field = np.asarray(field) + result = field + shift + if np.isclose(lmbda, 0): + array_to_lognormal(result) + if np.min(result) < -1 / lmbda: + warn("Box-Cox: Some values will be cut off!") + return (np.maximum(lmbda * result + 1, 0)) ** (1 / lmbda) + + +def array_zinnharvey(field, conn="high", mean=None, var=None): """ Zinn and Harvey transformation to connect low or high values. Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. conn : :class:`str`, optional Desired connectivity. Either "low" or "high". Default: "high" @@ -320,35 +775,32 @@ def _zinnharvey(field, conn="high", mean=None, var=None): Default: :any:`None` var : :class:`float` or :any:`None`, optional Variance of the given field. - If None is given, the mean will be calculated. + If None is given, the variance will be calculated. Default: :any:`None` Returns ------- - :class:`numpy.ndarray` - Transformed field. - """ - if mean is None: - mean = np.mean(field) - if var is None: - var = np.var(field) - field = np.abs((field - mean) / np.sqrt(var)) - field = 2 * erf(field / np.sqrt(2)) - 1 - field = np.sqrt(2) * erfinv(field) + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + result = np.abs((field - mean) / np.sqrt(var)) + result = np.sqrt(2) * erfinv(2 * erf(result / np.sqrt(2)) - 1) if conn == "high": - field = -field - return field * np.sqrt(var) + mean + result = -result + return result * np.sqrt(var) + mean -def _normal_force_moments(field, mean=0, var=1): +def array_force_moments(field, mean=0, var=1): """ Force moments of a normal distributed field. Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. mean : :class:`float`, optional Desired mean of the field. Default: 0 @@ -358,63 +810,61 @@ def _normal_force_moments(field, mean=0, var=1): Returns ------- - :class:`numpy.ndarray` - Transformed field. + :class:`numpy.ndarray` + Transformed field. """ + field = np.asarray(field) var_in = np.var(field) mean_in = np.mean(field) rescale = np.sqrt(var / var_in) return rescale * (field - mean_in) + mean -def _normal_to_lognormal(field): +def array_to_lognormal(field): """ Transform normal distribution to log-normal distribution. Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. Returns ------- - :class:`numpy.ndarray` - Transformed field. + :class:`numpy.ndarray` + Transformed field. """ return np.exp(field) -def _normal_to_uniform(field, mean=None, var=None): +def array_to_uniform(field, mean=None, var=None): """ Transform normal distribution to uniform distribution on [0, 1]. Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. mean : :class:`float` or :any:`None`, optional Mean of the given field. If None is given, the mean will be calculated. Default: :any:`None` var : :class:`float` or :any:`None`, optional Variance of the given field. - If None is given, the mean will be calculated. + If None is given, the variance will be calculated. Default: :any:`None` Returns ------- - :class:`numpy.ndarray` - Transformed field. + :class:`numpy.ndarray` + Transformed field. """ - if mean is None: - mean = np.mean(field) - if var is None: - var = np.var(field) + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) return 0.5 * (1 + erf((field - mean) / np.sqrt(2 * var))) -def _normal_to_arcsin(field, mean=None, var=None, a=0, b=1): +def array_to_arcsin(field, mean=None, var=None, a=None, b=None): """ Transform normal distribution to arcsin distribution. @@ -423,8 +873,7 @@ def _normal_to_arcsin(field, mean=None, var=None, a=0, b=1): Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. mean : :class:`float` or :any:`None`, optional Mean of the given field. If None is given, the mean will be calculated. Default: :any:`None` @@ -433,19 +882,26 @@ def _normal_to_arcsin(field, mean=None, var=None, a=0, b=1): If None is given, the mean will be calculated. Default: :any:`None` a : :class:`float`, optional - Parameter a of the arcsin distribution. Default: 0 + Parameter a of the arcsin distribution (lower bound). + Default: keep mean and variance b : :class:`float`, optional - Parameter b of the arcsin distribution. Default: 1 + Parameter b of the arcsin distribution (upper bound). + Default: keep mean and variance Returns ------- - :class:`numpy.ndarray` - Transformed field. + :class:`numpy.ndarray` + Transformed field. """ - return _uniform_to_arcsin(_normal_to_uniform(field, mean, var), a, b) + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(2.0 * var) if a is None else float(a) + b = mean + np.sqrt(2.0 * var) if b is None else float(b) + return _uniform_to_arcsin(array_to_uniform(field, mean, var), a, b) -def _normal_to_uquad(field, mean=None, var=None, a=0, b=1): +def array_to_uquad(field, mean=None, var=None, a=0, b=1): """ Transform normal distribution to U-quadratic distribution. @@ -454,26 +910,32 @@ def _normal_to_uquad(field, mean=None, var=None, a=0, b=1): Parameters ---------- field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. + Normal distributed values. mean : :class:`float` or :any:`None`, optional Mean of the given field. If None is given, the mean will be calculated. Default: :any:`None` var : :class:`float` or :any:`None`, optional Variance of the given field. - If None is given, the mean will be calculated. + If None is given, the variance will be calculated. Default: :any:`None` a : :class:`float`, optional - Parameter a of the U-quadratic distribution. Default: 0 + Parameter a of the U-quadratic distribution (lower bound). + Default: keep mean and variance b : :class:`float`, optional - Parameter b of the U-quadratic distribution. Default: 1 + Parameter b of the U-quadratic distribution (upper bound). + Default: keep mean and variance Returns ------- - :class:`numpy.ndarray` - Transformed field. + :class:`numpy.ndarray` + Transformed field. """ - return _uniform_to_uquad(_normal_to_uniform(field, mean, var), a, b) + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(5.0 / 3.0 * var) if a is None else float(a) + b = mean + np.sqrt(5.0 / 3.0 * var) if b is None else float(b) + return _uniform_to_uquad(array_to_uniform(field, mean, var), a, b) def _uniform_to_arcsin(field, a=0, b=1): @@ -486,6 +948,7 @@ def _uniform_to_arcsin(field, a=0, b=1): in this case: the arcsin distribution See: https://en.wikipedia.org/wiki/Arcsine_distribution """ + field = np.asarray(field) return (b - a) * np.sin(np.pi * 0.5 * field) ** 2 + a @@ -499,11 +962,12 @@ def _uniform_to_uquad(field, a=0, b=1): in this case: the U-quadratic distribution See: https://en.wikipedia.org/wiki/U-quadratic_distribution """ + field = np.asarray(field) al = 12 / (b - a) ** 3 be = (a + b) / 2 ga = (a - b) ** 3 / 8 y_raw = 3 * field / al + ga - out = np.zeros_like(y_raw) - out[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) - out[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) - return out + be + result = np.zeros_like(y_raw) + result[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) + result[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) + return result + be From a6fec2ccf6f281900ff275bb597094e35bf95bea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 17:50:22 +0200 Subject: [PATCH 02/49] transform: minor fixes --- gstools/transform/field.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/gstools/transform/field.py b/gstools/transform/field.py index efd7387c9..469b75f7f 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -76,7 +76,7 @@ def _get_field_io(field, store): if not isinstance(field, str): raise ValueError( - f"transform: given field name '{field}' is not a string." + f"transform: given field name '{field}' is not a string" ) if isinstance(store, str): if not store.isidentifier(): @@ -193,7 +193,7 @@ def apply(fld, method, field="field", store=True, process=False, **kwargs): return normal_to_uquad(fld, **kwargs) if method.endswith("function"): return apply_function(fld, **kwargs) - raise ValueError(f"transform.apply: unknown method '{method}'.") + raise ValueError(f"transform.apply: unknown method '{method}'") def apply_function( @@ -239,7 +239,7 @@ def apply_function( Transformed field. """ if not callable(function): - raise ValueError("apply_function: function not callable.") + raise ValueError("apply_function: function not callable") output_field, save = _get_field_io(field, store) data = _check_input_field(fld, field) if process: @@ -706,15 +706,14 @@ def array_discrete( else: if len(values) != len(thresholds) + 1: raise ValueError( - "discrete transformation: " - "len(values) != len(thresholds) + 1" + "discrete transformation: len(values) != len(thresholds) + 1" ) values = np.array(values) thresholds = np.array(thresholds) # check thresholds if not np.all(thresholds[:-1] < thresholds[1:]): raise ValueError( - "discrete transformation: thresholds need to be ascending." + "discrete transformation: thresholds need to be ascending" ) # use a separate result so the intermediate results are not affected result = np.empty_like(field) @@ -901,7 +900,7 @@ def array_to_arcsin(field, mean=None, var=None, a=None, b=None): return _uniform_to_arcsin(array_to_uniform(field, mean, var), a, b) -def array_to_uquad(field, mean=None, var=None, a=0, b=1): +def array_to_uquad(field, mean=None, var=None, a=None, b=None): """ Transform normal distribution to U-quadratic distribution. From eb67a8a71d15f8a753833b1ba54efb68ec7b4654 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 17:50:57 +0200 Subject: [PATCH 03/49] Field: add 'transform' method --- gstools/field/base.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/gstools/field/base.py b/gstools/field/base.py index 7c83f370b..3ee23f120 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -20,6 +20,7 @@ fmt_mean_norm_trend, ) from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer +from gstools.transform.field import apply __all__ = ["Field"] @@ -255,6 +256,41 @@ def post_field(self, field, name="field", process=True, save=True): setattr(self, name, field) return field + def transform( + self, method, field="field", store=True, process=False, **kwargs + ): + """ + Apply field transformation. + + Parameters + ---------- + method : :class:`str` + Method to use. See :any:`transform` for available transformations. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + **kwargs + Keyword arguments forwarded to selected method. + + Raises + ------ + ValueError + When method is unknown. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + return apply( + self, method, field=field, store=store, process=process, **kwargs + ) + def to_pyvista( self, field_select="field", fieldname="field" ): # pragma: no cover From f6c669875719ad295bc95fa205b08a6f49d8afde Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 18:24:02 +0200 Subject: [PATCH 04/49] transform: update examples --- examples/07_transformations/00_log_normal.py | 4 +++- examples/07_transformations/01_binary.py | 4 +++- examples/07_transformations/02_discrete.py | 8 ++++--- examples/07_transformations/03_zinn_harvey.py | 4 +++- examples/07_transformations/04_bimodal.py | 4 +++- .../07_transformations/05_combinations.py | 19 ++++++++++++----- examples/07_transformations/README.rst | 21 +++++++++++++++---- 7 files changed, 48 insertions(+), 16 deletions(-) diff --git a/examples/07_transformations/00_log_normal.py b/examples/07_transformations/00_log_normal.py index a44f50fb2..a02673b2a 100755 --- a/examples/07_transformations/00_log_normal.py +++ b/examples/07_transformations/00_log_normal.py @@ -3,6 +3,8 @@ ----------------- Here we transform a field to a log-normal distribution: + +See :any:`transform.normal_to_lognormal` """ import gstools as gs @@ -11,5 +13,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.normal_to_lognormal(srf) +srf.transform("lognormal") srf.plot() diff --git a/examples/07_transformations/01_binary.py b/examples/07_transformations/01_binary.py index 20b419191..a403abb7f 100755 --- a/examples/07_transformations/01_binary.py +++ b/examples/07_transformations/01_binary.py @@ -5,6 +5,8 @@ Here we transform a field to a binary field with only two values. The dividing value is the mean by default and the upper and lower values are derived to preserve the variance. + +See :any:`transform.binary` """ import gstools as gs @@ -13,5 +15,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.binary(srf) +srf.transform("binary") srf.plot() diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py index 325018984..62566ef86 100755 --- a/examples/07_transformations/02_discrete.py +++ b/examples/07_transformations/02_discrete.py @@ -6,6 +6,8 @@ If we do not give thresholds, the pairwise means of the given values are taken as thresholds. If thresholds are given, arbitrary values can be applied to the field. + +See :any:`transform.discrete` """ import numpy as np import gstools as gs @@ -18,14 +20,14 @@ # create 5 equidistanly spaced values, thresholds are the arithmetic means srf.structured([x, y]) discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5) -gs.transform.discrete(srf, discrete_values) +srf.transform("discrete", values=discrete_values) srf.plot() # calculate thresholds for equal shares # but apply different values to the separated classes discrete_values2 = [0, -1, 2, -3, 4] srf.structured([x, y]) -gs.transform.discrete(srf, discrete_values2, thresholds="equal") +srf.transform("discrete", values=discrete_values2, thresholds="equal") srf.plot() # user defined thresholds @@ -33,5 +35,5 @@ # apply different values to the separated classes discrete_values3 = [0, 1, 10] srf.structured([x, y]) -gs.transform.discrete(srf, discrete_values3, thresholds=thresholds) +srf.transform("discrete", values=discrete_values3, thresholds=thresholds) srf.plot() diff --git a/examples/07_transformations/03_zinn_harvey.py b/examples/07_transformations/03_zinn_harvey.py index d3d28a1a7..01e65790c 100755 --- a/examples/07_transformations/03_zinn_harvey.py +++ b/examples/07_transformations/03_zinn_harvey.py @@ -6,6 +6,8 @@ `Zinn & Harvey (2003) `__. With this transformation, one could overcome the restriction that in ordinary Gaussian random fields the mean values are the ones being the most connected. + +See :any:`transform.zinnharvey` """ import gstools as gs @@ -14,5 +16,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.zinnharvey(srf, conn="high") +srf.transform("zinnharvey", conn="high") srf.plot() diff --git a/examples/07_transformations/04_bimodal.py b/examples/07_transformations/04_bimodal.py index e0d4b2ded..1f37d3abb 100755 --- a/examples/07_transformations/04_bimodal.py +++ b/examples/07_transformations/04_bimodal.py @@ -8,6 +8,8 @@ * `uquad `__. Both transformations will preserve the mean and variance of the given field by default. + +See: :any:`transform.normal_to_arcsin` and :any:`transform.normal_to_uquad` """ import gstools as gs @@ -16,5 +18,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) field = srf.structured([x, y]) -gs.transform.normal_to_arcsin(srf) +srf.transform("arcsin") srf.plot() diff --git a/examples/07_transformations/05_combinations.py b/examples/07_transformations/05_combinations.py index 2edb1a180..1e2767ca2 100755 --- a/examples/07_transformations/05_combinations.py +++ b/examples/07_transformations/05_combinations.py @@ -9,7 +9,11 @@ Then we apply the Zinn & Harvey transformation to connect the low values. Afterwards the field is transformed to a binary field and last but not least, we transform it to log-values. + +We can select the desired field by its name and we can define an output name +to store the field. """ +# sphinx_gallery_thumbnail_number = 1 import gstools as gs # structured field with a size of 100x100 and a grid-size of 1x1 @@ -17,13 +21,18 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, mean=-9, seed=20170519) srf.structured([x, y]) -gs.transform.normal_force_moments(srf) -gs.transform.zinnharvey(srf, conn="low") -gs.transform.binary(srf) -gs.transform.normal_to_lognormal(srf) -srf.plot() +srf.transform("force_moments", field="field", store="f_forced") +srf.transform("zinnharvey", field="f_forced", store="f_zinnharvey", conn="low") +srf.transform("binary", field="f_zinnharvey", store="f_binary") +srf.transform("lognormal", field="f_binary", store="f_result") +srf.plot(field="f_result") ############################################################################### # The resulting field could be interpreted as a transmissivity field, where # the values of low permeability are the ones being the most connected # and only two kinds of soil exist. +# +# All stored fields can be accessed and plotted by name: + +print("Max binary value:", srf.f_binary.max()) +srf.plot(field="f_zinnharvey") diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst index 2be183ee0..d93c99307 100644 --- a/examples/07_transformations/README.rst +++ b/examples/07_transformations/README.rst @@ -20,18 +20,31 @@ common transformations: normal_to_uniform normal_to_arcsin normal_to_uquad + apply_function All the transformations take a field class, that holds a generated field, -as input and will manipulate this field inplace. +as input and will manipulate this field inplace or store it with a given name. -Simply import the transform submodule and apply a transformation to the srf class: +Simply apply a transformation to a field class: .. code-block:: python - from gstools import transform as tf + import gstools as gs ... - tf.normal_to_lognormal(srf) + srf = gs.SRF(model) + srf(...) + gs.transform.normal_to_lognormal(srf) + +Or use the provided wrapper: + +.. code-block:: python + + import gstools as gs + ... + srf = gs.SRF(model) + srf(...) + srf.transform("lognormal") Examples -------- From 0b7978b85216bc0783234c5197a6c3c1f739d2c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 18:37:00 +0200 Subject: [PATCH 05/49] transform: update discrete example to store fields --- examples/07_transformations/02_discrete.py | 36 ++++++++++++---------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py index 62566ef86..02f8ef579 100755 --- a/examples/07_transformations/02_discrete.py +++ b/examples/07_transformations/02_discrete.py @@ -12,28 +12,32 @@ import numpy as np import gstools as gs -# structured field with a size of 100x100 and a grid-size of 0.5x0.5 +# Structured field with a size of 100x100 and a grid-size of 0.5x0.5 x = y = np.arange(200) * 0.5 model = gs.Gaussian(dim=2, var=1, len_scale=5) srf = gs.SRF(model, seed=20170519) - -# create 5 equidistanly spaced values, thresholds are the arithmetic means srf.structured([x, y]) -discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5) -srf.transform("discrete", values=discrete_values) -srf.plot() +############################################################################### +# Create 5 equidistanly spaced values, thresholds are the arithmetic means + +values1 = np.linspace(np.min(srf.field), np.max(srf.field), 5) +srf.transform("discrete", store="f1", values=values1) +srf.plot("f1") + +############################################################################### # calculate thresholds for equal shares # but apply different values to the separated classes -discrete_values2 = [0, -1, 2, -3, 4] -srf.structured([x, y]) -srf.transform("discrete", values=discrete_values2, thresholds="equal") -srf.plot() -# user defined thresholds +values2 = [0, -1, 2, -3, 4] +srf.transform("discrete", store="f2", values=values2, thresholds="equal") +srf.plot("f2") + +############################################################################### +# Create user defined thresholds +# and apply different values to the separated classes + +values3 = [0, 1, 10] thresholds = [-1, 1] -# apply different values to the separated classes -discrete_values3 = [0, 1, 10] -srf.structured([x, y]) -srf.transform("discrete", values=discrete_values3, thresholds=thresholds) -srf.plot() +srf.transform("discrete", store="f3", values=values3, thresholds=thresholds) +srf.plot("f3") From abbb86e1df8c4427dcb5d15d09b0664f765d44ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 19:10:42 +0200 Subject: [PATCH 06/49] transform: minor example fixes --- examples/07_transformations/00_log_normal.py | 2 +- examples/07_transformations/02_discrete.py | 2 +- examples/07_transformations/04_bimodal.py | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/07_transformations/00_log_normal.py b/examples/07_transformations/00_log_normal.py index a02673b2a..3d2f45325 100755 --- a/examples/07_transformations/00_log_normal.py +++ b/examples/07_transformations/00_log_normal.py @@ -13,5 +13,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -srf.transform("lognormal") +srf.transform("normal_to_lognormal") # also "lognormal" works srf.plot() diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py index 02f8ef579..a490fb5b5 100755 --- a/examples/07_transformations/02_discrete.py +++ b/examples/07_transformations/02_discrete.py @@ -26,7 +26,7 @@ srf.plot("f1") ############################################################################### -# calculate thresholds for equal shares +# Calculate thresholds for equal shares # but apply different values to the separated classes values2 = [0, -1, 2, -3, 4] diff --git a/examples/07_transformations/04_bimodal.py b/examples/07_transformations/04_bimodal.py index 1f37d3abb..8234a486a 100755 --- a/examples/07_transformations/04_bimodal.py +++ b/examples/07_transformations/04_bimodal.py @@ -1,6 +1,6 @@ """ -bimodal fields ---------------- +Bimodal fields +-------------- We provide two transformations to obtain bimodal distributions: @@ -18,5 +18,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) field = srf.structured([x, y]) -srf.transform("arcsin") +srf.transform("normal_to_arcsin") # also "arcsin" works srf.plot() From 6afabd332c9a809de8b7fc712905bb99bb40fc8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 1 Aug 2021 19:12:58 +0200 Subject: [PATCH 07/49] transform: add clarification to example --- examples/07_transformations/05_combinations.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/07_transformations/05_combinations.py b/examples/07_transformations/05_combinations.py index 1e2767ca2..fb27c5e69 100755 --- a/examples/07_transformations/05_combinations.py +++ b/examples/07_transformations/05_combinations.py @@ -12,6 +12,8 @@ We can select the desired field by its name and we can define an output name to store the field. + +If you don't specify `field` and `store` everything happens inplace. """ # sphinx_gallery_thumbnail_number = 1 import gstools as gs From 128b262263f419f516586318e0d975cf637cb3ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 12:11:33 +0200 Subject: [PATCH 08/49] Field: allow pos=None to reuse present pos tuple --- gstools/field/base.py | 45 +++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 3ee23f120..695f0cdb6 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -89,7 +89,7 @@ def __init__( self.trend = trend def __call__( - self, pos, field=None, mesh_type="unstructured", post_process=True + self, pos=None, field=None, mesh_type="unstructured", process=True ): """Generate the field. @@ -102,7 +102,7 @@ def __call__( the field values. Will be all zeros if :any:`None` is given. mesh_type : :class:`str`, optional 'structured' / 'unstructured'. Default: 'unstructured' - post_process : :class:`bool`, optional + process : :class:`bool`, optional Whether to apply mean, normalizer and trend to the field. Default: `True` @@ -116,13 +116,15 @@ def __call__( field = np.zeros(shape, dtype=np.double) else: field = np.array(field, dtype=np.double).reshape(shape) - return self.post_field(field, process=post_process) + return self.post_field(field, process=process) def structured(self, *args, **kwargs): """Generate a field on a structured mesh. See :any:`__call__` """ + if not (args or "pos" in kwargs) and self.mesh_type == "unstructured": + raise ValueError("Field.structured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="structured") return call(*args, **kwargs) @@ -131,6 +133,8 @@ def unstructured(self, *args, **kwargs): See :any:`__call__` """ + if not (args or "pos" in kwargs) and self.mesh_type != "unstructured": + raise ValueError("Field.unstructured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="unstructured") return call(*args, **kwargs) @@ -174,7 +178,7 @@ def mesh( """ return generate_on_mesh(self, mesh, points, direction, name, **kwargs) - def pre_pos(self, pos, mesh_type="unstructured"): + def pre_pos(self, pos=None, mesh_type="unstructured"): """ Preprocessing positions and mesh_type. @@ -194,17 +198,30 @@ def pre_pos(self, pos, mesh_type="unstructured"): shape : :class:`tuple` Shape of the resulting field. """ - # save mesh-type - self.mesh_type = mesh_type - # save pos tuple - if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, self.dim) - self.pos = pos - pos = generate_grid(pos) + if pos is None: + if self.pos is None: + raise ValueError("Field: no position tuple 'pos' present") + if self.mesh_type is None: + raise ValueError("Field: no 'mesh_type' present") + pos = self.pos + mesh_type = self.mesh_type + if mesh_type != "unstructured": + pos = generate_grid(pos) + shape = tuple(len(p_i) for p_i in pos) + else: + shape = np.shape(pos[0]) else: - pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) - self.pos = pos - shape = np.shape(pos[0]) + # save mesh-type + self.mesh_type = mesh_type + # save pos tuple + if mesh_type != "unstructured": + pos, shape = format_struct_pos_dim(pos, self.dim) + self.pos = pos + pos = generate_grid(pos) + else: + pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) + self.pos = pos + shape = np.shape(pos[0]) # prepend dimension if we have a vector field if self.value_type == "vector": shape = (self.dim,) + shape From 7094cae8f363571265471a7aaf61c08e3c685592 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 17:46:57 +0200 Subject: [PATCH 09/49] Field: enable setting fields by name; allow reusing pos-tuple; enable subscription --- gstools/field/base.py | 199 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 164 insertions(+), 35 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 695f0cdb6..06af503a6 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -11,6 +11,8 @@ """ # pylint: disable=C0103, C0415 from functools import partial +from collections.abc import Iterable +from copy import copy import numpy as np from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, generate_grid @@ -18,6 +20,7 @@ generate_on_mesh, to_vtk_helper, fmt_mean_norm_trend, + _names, ) from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer from gstools.transform.field import apply @@ -28,6 +31,19 @@ """:class:`list` of :class:`str`: valid field value types.""" +def _pos_equal(pos1, pos2): + if pos1 is None or pos2 is None: + return False + if len(pos1) != len(pos2): + return False + for p1, p2 in zip(pos1, pos2): + if len(p1) != len(p2): + return False + if not np.allclose(p1, p2): + return False + return True + + def _set_mean_trend(value, dim): if callable(value) or value is None: return value @@ -71,10 +87,10 @@ def __init__( dim=None, ): # initialize attributes - self.pos = None - self.mesh_type = None - self.field = None - # initialize private attributes + self._mesh_type = "unstructured" # default + self._pos = None + self._field_shape = None + self._field_names = [] self._model = None self._value_type = None self._mean = None @@ -88,8 +104,18 @@ def __init__( self.normalizer = normalizer self.trend = trend + def __getitem__(self, field): + if field in self.field_names: + return getattr(self, field) + raise ValueError(f"Field: requested field '{field}' not present") + def __call__( - self, pos=None, field=None, mesh_type="unstructured", process=True + self, + pos=None, + field=None, + mesh_type="unstructured", + post_process=True, + store=True, ): """Generate the field. @@ -102,27 +128,34 @@ def __call__( the field values. Will be all zeros if :any:`None` is given. mesh_type : :class:`str`, optional 'structured' / 'unstructured'. Default: 'unstructured' - process : :class:`bool`, optional + post_process : :class:`bool`, optional Whether to apply mean, normalizer and trend to the field. Default: `True` + store : :class:`str` or :class:`bool`, optional + Whether to store field (True/False) with default name + or with specified name. + The default is :any:`True` for default name "field". Returns ------- field : :class:`numpy.ndarray` the field values. """ + name, save = self._get_store_config(store) pos, shape = self.pre_pos(pos, mesh_type) if field is None: field = np.zeros(shape, dtype=np.double) else: field = np.array(field, dtype=np.double).reshape(shape) - return self.post_field(field, process=process) + return self.post_field(field, name, post_process, save) def structured(self, *args, **kwargs): """Generate a field on a structured mesh. See :any:`__call__` """ + if self.pos is None: + self.mesh_type = "structured" if not (args or "pos" in kwargs) and self.mesh_type == "unstructured": raise ValueError("Field.structured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="structured") @@ -133,6 +166,8 @@ def unstructured(self, *args, **kwargs): See :any:`__call__` """ + if self.pos is None: + self.mesh_type = "unstructured" if not (args or "pos" in kwargs) and self.mesh_type != "unstructured": raise ValueError("Field.unstructured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="unstructured") @@ -178,7 +213,7 @@ def mesh( """ return generate_on_mesh(self, mesh, points, direction, name, **kwargs) - def pre_pos(self, pos=None, mesh_type="unstructured"): + def pre_pos(self, pos=None, mesh_type="unstructured", info=False): """ Preprocessing positions and mesh_type. @@ -190,47 +225,42 @@ def pre_pos(self, pos=None, mesh_type="unstructured"): mesh_type : :class:`str`, optional 'structured' / 'unstructured' Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information Returns ------- iso_pos : (d, n), :class:`numpy.ndarray` - the isometrized position tuple + Isometrized position tuple. shape : :class:`tuple` Shape of the resulting field. + info : :class:`dict` + Information about settings. """ + info_ret = {"deleted": False} if pos is None: if self.pos is None: raise ValueError("Field: no position tuple 'pos' present") if self.mesh_type is None: raise ValueError("Field: no 'mesh_type' present") - pos = self.pos - mesh_type = self.mesh_type - if mesh_type != "unstructured": - pos = generate_grid(pos) - shape = tuple(len(p_i) for p_i in pos) - else: - shape = np.shape(pos[0]) else: - # save mesh-type - self.mesh_type = mesh_type - # save pos tuple - if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, self.dim) - self.pos = pos - pos = generate_grid(pos) - else: - pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) - self.pos = pos - shape = np.shape(pos[0]) - # prepend dimension if we have a vector field - if self.value_type == "vector": - shape = (self.dim,) + shape - if self.latlon: - raise ValueError("Field: Vector fields not allowed for latlon") - # return isometrized pos tuple and resulting field shape + old_pos = copy(self.pos) + # save pos and mesh-type + self.set_pos(pos, mesh_type) + # remove present fields if new pos is different from current + if not _pos_equal(old_pos, self.pos): + self.delete_fields() + info_ret["deleted"] = True + del old_pos + if self.mesh_type != "unstructured": + pos = generate_grid(self.pos) + else: + pos = self.pos + # return isometrized pos tuple, field shape and possible info + info_ret = (info_ret,) if self.model is None: - return pos, shape - return self.model.isometrize(pos), shape + return (pos, self.field_shape) + info * info_ret + return (self.model.isometrize(pos), self.field_shape) + info * info_ret def post_field(self, field, name="field", process=True, save=True): """ @@ -270,9 +300,25 @@ def post_field(self, field, name="field", process=True, save=True): stacked=False, ) if save: + name = str(name) + if not name.isidentifier() or ( + name not in self.field_names and name in dir(self) + ): + raise ValueError( + f"Field: given field name '{name}' is not valid" + ) + # allow resetting present fields + if name not in self._field_names: + self._field_names.append(name) setattr(self, name, field) return field + def delete_fields(self): + """Delete all present fields""" + for field in self._field_names: + delattr(self, field) + self._field_names = [] + def transform( self, method, field="field", store=True, process=False, **kwargs ): @@ -396,6 +442,89 @@ def plot( return r + def set_pos(self, pos, mesh_type="unstructured"): + """ + Set positions and mesh_type. + + Parameters + ---------- + pos : :any:`iterable` + the position tuple, containing main direction and transversal + directions + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: `"unstructured"` + """ + self.mesh_type = mesh_type + self.pos = pos + + @staticmethod + def _get_store_config(store, default="field", fld_cnt=None): + # single field + if fld_cnt is None: + save = isinstance(store, str) or bool(store) + name = store if isinstance(store, str) else default + return name, save + # multiple fields + default = _names(default, fld_cnt) + save = [True] * fld_cnt + if isinstance(store, str): + store = [store] + if isinstance(store, Iterable): + store = list(store)[:fld_cnt] + store += [True] * (fld_cnt - len(store)) + name = [None] * fld_cnt + for i, val in enumerate(store): + save[i] = isinstance(val, str) or bool(val) + name[i] = val if isinstance(val, str) else default[i] + else: + save = [bool(store)] * fld_cnt + name = copy(default) + return name, save + + @property + def pos(self): + """:class:`tuple`: The position tuple of the field.""" + return self._pos + + @pos.setter + def pos(self, pos): + if self.mesh_type is None: + raise ValueError("Field.pos: can't set 'pos' without 'mesh_type'") + if self.mesh_type == "unstructured": + self._pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) + self._field_shape = np.shape(self._pos[0]) + else: + self._pos, self._field_shape = format_struct_pos_dim(pos, self.dim) + # prepend dimension if we have a vector field + if self.value_type == "vector": + self._field_shape = (self.dim,) + self._field_shape + if self.latlon: + raise ValueError("Field: Vector fields not allowed for latlon") + + @property + def field_names(self): + """:class:`list`: Names of present fields.""" + return self._field_names + + @field_names.deleter + def field_names(self): + self.delete_fields() + + @property + def field_shape(self): + """:class:`tuple`: The shape of the field.""" + return self._field_shape + + @property + def mesh_type(self): + """:class:`str`: The mesh type of the field.""" + return self._mesh_type + + @mesh_type.setter + def mesh_type(self, mesh_type): + self._mesh_type = str(mesh_type) + @property def model(self): """:any:`CovModel`: The covariance model of the field.""" From 6582331d1156a1a2cba6e07b7407d7e085e9719b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 17:48:14 +0200 Subject: [PATCH 10/49] CondSRF: enable storage control; better integration with kriging class (reuse pos, meshtype); only recalc krige if needed --- gstools/field/cond_srf.py | 95 ++++++++++++++++++++++++++++++++++----- 1 file changed, 83 insertions(+), 12 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 609045c83..04479c257 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -48,16 +48,22 @@ def __init__(self, krige, generator="RandMeth", **generator_kwargs): raise ValueError("CondSRF: krige should be an instance of Krige.") self._krige = krige # initialize attributes - self.pos = None - self.mesh_type = None - self.field = None - self.raw_field = None + self._field_names = [] # initialize private attributes self._generator = None # initialize attributes self.set_generator(generator, **generator_kwargs) - def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): + def __call__( + self, + pos=None, + seed=np.nan, + mesh_type="unstructured", + post_process=True, + store=True, + krige_store=True, + **kwargs, + ): """Generate the conditioned spatial random field. The field is saved as `self.field` and is also returned. @@ -71,6 +77,19 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): seed for RNG for reseting. Default: keep seed from generator mesh_type : :class:`str` 'structured' / 'unstructured' + post_process : :class:`bool`, optional + Whether to apply mean, normalizer and trend to the field. + Default: `True` + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store fields (True/False) with default names + or with specified names. + The default is :any:`True` for default names + ["field", "raw_field", "raw_krige"]. + krige_store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store kriging fields (True/False) with default name + or with specified names. + The default is :any:`True` for default names + ["field", "krige_var"]. **kwargs keyword arguments that are forwarded to the kriging routine in use. @@ -79,23 +98,52 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): field : :class:`numpy.ndarray` the conditioned SRF """ + name, save = self._get_store_config( + store, default=["field", "raw_field", "raw_krige"], fld_cnt=3 + ) + krige_name, krige_save = self._get_store_config( + krige_store, default=["field", "krige_var"], fld_cnt=2 + ) kwargs["mesh_type"] = mesh_type kwargs["only_mean"] = False # overwrite if given kwargs["return_var"] = True # overwrite if given kwargs["post_process"] = False # overwrite if given + kwargs["store"] = [False, krige_name[1] if krige_save[1] else False] # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape - iso_pos, shape = self.pre_pos(pos, mesh_type) + iso_pos, shape, info = self.krige.pre_pos(pos, mesh_type, info=True) + # also reset fields in this class (not only krige with call above) + if info["deleted"]: + self.delete_fields() # generate the field - self.raw_field = np.reshape( - self.generator(iso_pos, add_nugget=False), shape - ) - field, krige_var = self.krige(pos, **kwargs) + rawfield = np.reshape(self.generator(iso_pos, add_nugget=False), shape) + # call krige on already set pos (reuse already calculated fields) + if ( + not info["deleted"] + and name[2] in self.field_names + and krige_name[1] in self.krige.field_names + ): + rawkrige, krige_var = self[name[2]], self.krige[krige_name[1]] + else: + rawkrige, krige_var = self.krige(**kwargs) var_scale, nugget = self.get_scaling(krige_var, shape) # need to use a copy to not alter "field" by reference - self.krige.post_field(self.krige.field.copy()) - return self.post_field(field + var_scale * self.raw_field + nugget) + # store krige field + self.krige.post_field( + rawkrige.copy(), krige_name[0], post_process, krige_save[0] + ) + # store raw krige field + self.post_field(rawkrige.copy(), name[2], False, save[2]) + # store raw random field + self.post_field(rawfield, name[1], False, save[1]) + # store cond random field + return self.post_field( + field=rawkrige + var_scale * rawfield + nugget, + name=name[0], + process=post_process, + save=save[0], + ) def get_scaling(self, krige_var, shape): """ @@ -143,6 +191,29 @@ def set_generator(self, generator, **generator_kwargs): else: raise ValueError(f"gstools.CondSRF: Unknown generator {generator}") + @property + def pos(self): + """:class:`tuple`: The position tuple of the field.""" + return self.krige.pos + + @pos.setter + def pos(self, pos): + self.krige.pos = pos + + @property + def field_shape(self): + """:class:`tuple`: The shape of the field.""" + return self.krige.field_shape + + @property + def mesh_type(self): + """:class:`str`: The mesh type of the field.""" + return self.krige.mesh_type + + @mesh_type.setter + def mesh_type(self, mesh_type): + self.krige.mesh_type = mesh_type + @property def krige(self): """:any:`Krige`: The underlying kriging class.""" From 1c0d385aee69e67d61c36537643eab9434872982 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 17:48:40 +0200 Subject: [PATCH 11/49] SRF: enable storage control --- gstools/field/srf.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index de0951241..c572ee8a8 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -100,7 +100,13 @@ def __init__( self.set_generator(generator, **generator_kwargs) def __call__( - self, pos, seed=np.nan, point_volumes=0.0, mesh_type="unstructured" + self, + pos=None, + seed=np.nan, + point_volumes=0.0, + mesh_type="unstructured", + post_process=True, + store=True, ): """Generate the spatial random field. @@ -121,12 +127,20 @@ def __call__( is changed. Default: ``0`` mesh_type : :class:`str` 'structured' / 'unstructured' + post_process : :class:`bool`, optional + Whether to apply mean, normalizer and trend to the field. + Default: `True` + store : :class:`str` or :class:`bool`, optional + Whether to store field (True/False) with default name + or with specified name. + The default is :any:`True` for default name "field". Returns ------- field : :class:`numpy.ndarray` the SRF """ + name, save = self._get_store_config(store) # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape @@ -139,7 +153,7 @@ def __call__( if np.size(scaled_var) > 1: scaled_var = np.reshape(scaled_var, shape) field *= np.sqrt(scaled_var / self.model.sill) - return self.post_field(field) + return self.post_field(field, name, post_process, save) def upscaling_func(self, *args, **kwargs): """Upscaling method applied to the field variance.""" From 585bdbdab0fa1e2f2262432937e6c9ae5991563a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 17:49:52 +0200 Subject: [PATCH 12/49] Krige: enable storage control --- gstools/krige/base.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 80cc21b49..40d387d47 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -144,8 +144,6 @@ def __init__( fit_variogram=False, ): super().__init__(model, mean=mean, normalizer=normalizer, trend=trend) - self.mean_field = None - self.krige_var = None self._unbiased = bool(unbiased) self._exact = bool(exact) self._pseudo_inv = bool(pseudo_inv) @@ -172,13 +170,14 @@ def __init__( def __call__( self, - pos, + pos=None, mesh_type="unstructured", ext_drift=None, chunk_size=None, only_mean=False, return_var=True, post_process=True, + store=True, ): """ Generate the kriging field. @@ -208,6 +207,11 @@ def __call__( post_process : :class:`bool`, optional Whether to apply mean, normalizer and trend to the field. Default: `True` + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store kriging fields (True/False) with default name + or with specified names. + The default is :any:`True` for default names + ["field", "krige_var"] or "mean_field" if `only_mean=True`. Returns ------- @@ -218,6 +222,11 @@ def __call__( (if return_var is True and only_mean is False) """ return_var &= not only_mean # don't return variance when calc. mean + fld_cnt = 2 if return_var else 1 + default = ["mean_field"] if only_mean else ["field", "krige_var"] + name, save = self._get_store_config( + store, default=default, fld_cnt=fld_cnt + ) iso_pos, shape = self.pre_pos(pos, mesh_type) pnt_cnt = len(iso_pos[0]) @@ -248,18 +257,14 @@ def __call__( self._summate(field, krige_var, c_slice, k_vec, return_var) # reshape field if we got a structured mesh field = np.reshape(field, shape) - if only_mean: # care about 'kriging the mean' - return self.post_field(field, "mean_field", process=post_process) # save field to class - field = self.post_field(field, "field", process=post_process) + field = self.post_field(field, name[0], post_process, save[0]) if return_var: # care about the estimated error variance krige_var = np.reshape( np.maximum(self.model.sill - krige_var, 0), shape ) - krige_var = self.post_field(krige_var, "krige_var", process=False) + krige_var = self.post_field(krige_var, name[1], False, save[1]) return field, krige_var - # if we only calculate the field, overwrite the error variance - self.krige_var = None return field def _summate(self, field, krige_var, c_slice, k_vec, return_var): From 07b34d4422806a3255bba7df42195f0ce758ef3c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 17:50:36 +0200 Subject: [PATCH 13/49] Examples: update ensemble example to use new features --- .../06_conditioned_fields/00_condition_ensemble.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/06_conditioned_fields/00_condition_ensemble.py b/examples/06_conditioned_fields/00_condition_ensemble.py index a622fb463..ad11c0c32 100644 --- a/examples/06_conditioned_fields/00_condition_ensemble.py +++ b/examples/06_conditioned_fields/00_condition_ensemble.py @@ -26,15 +26,18 @@ model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5) krige = gs.krige.Ordinary(model, cond_pos, cond_val) cond_srf = gs.CondSRF(krige) +cond_srf.set_pos(gridx) ############################################################################### +# We can specify individual names for each field by the keyword `store`: -fields = [] for i in range(100): - fields.append(cond_srf(gridx, seed=i)) + cond_srf(seed=i, store=f"f{i}") label = "Conditioned ensemble" if i == 0 else None - plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label) -plt.plot(gridx, cond_srf.krige(gridx, only_mean=True), label="estimated mean") + plt.plot(gridx, cond_srf[f"f{i}"], color="k", alpha=0.1, label=label) + +fields = [cond_srf[f"f{i}"] for i in range(100)] +plt.plot(gridx, cond_srf.krige(only_mean=True), label="estimated mean") plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean") plt.plot(gridx, cond_srf.krige.field, linestyle="dashed", label="kriged field") plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") From 007aa2252f286970a1da1ad55b808eb83f12dd62 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 18:14:24 +0200 Subject: [PATCH 14/49] Transform: refactor --- gstools/transform/field.py | 34 +++++----------------------------- 1 file changed, 5 insertions(+), 29 deletions(-) diff --git a/gstools/transform/field.py b/gstools/transform/field.py index 469b75f7f..ad8b32cb4 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -73,30 +73,6 @@ ] -def _get_field_io(field, store): - if not isinstance(field, str): - raise ValueError( - f"transform: given field name '{field}' is not a string" - ) - if isinstance(store, str): - if not store.isidentifier(): - raise ValueError( - f"transform: given store name '{store}' is not valid" - ) - save = True - output_field = store - else: - save = bool(store) - output_field = field - return output_field, save - - -def _check_input_field(fld, field): - if not hasattr(fld, field) or getattr(fld, field) is None: - raise ValueError(f"transform: selected field '{field}' not present") - return getattr(fld, field) - - def _pre_process(fld, data, keep_mean): return remove_trend_norm_mean( pos=fld.pos, @@ -151,7 +127,7 @@ def apply(fld, method, field="field", store=True, process=False, **kwargs): field : :class:`str`, optional Name of field to be transformed. The default is "field". store : :class:`str` or :class:`bool`, optional - Whether to store field inplace (True/False) or under a given name. + Whether to store field inplace (True/False) or with a specified name. The default is True. process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean @@ -239,15 +215,15 @@ def apply_function( Transformed field. """ if not callable(function): - raise ValueError("apply_function: function not callable") - output_field, save = _get_field_io(field, store) - data = _check_input_field(fld, field) + raise ValueError("transform.apply_function: function not a 'callable'") + data = fld[field] + name, save = fld._get_store_config(store, default=field) if process: data = _pre_process(fld, data, keep_mean=keep_mean) data = function(data, **kwargs) if process: data = _post_process(fld, data, keep_mean=keep_mean) - return fld.post_field(data, name=output_field, process=False, save=save) + return fld.post_field(data, name=name, process=False, save=save) def binary( From aa7e6fdb0182c2ddc574567cb9466ee365a2bdd5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 18:22:04 +0200 Subject: [PATCH 15/49] Field: use new features in tools --- gstools/field/plot.py | 10 +++------- gstools/field/tools.py | 9 +-------- 2 files changed, 4 insertions(+), 15 deletions(-) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index ca3a422b4..c0b589b21 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -51,12 +51,10 @@ def plot_field( **kwargs Forwarded to the plotting routine. """ - plt_fld = getattr(fld, field) - assert not (fld.pos is None or plt_fld is None) if fld.dim == 1: - return plot_1d(fld.pos, plt_fld, fig, ax, **kwargs) + return plot_1d(fld.pos, fld[field], fig, ax, **kwargs) return plot_nd( - fld.pos, plt_fld, fld.mesh_type, fig, ax, fld.latlon, **kwargs + fld.pos, fld[field], fld.mesh_type, fig, ax, fld.latlon, **kwargs ) @@ -302,9 +300,7 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover "Only structured vector fields are supported " "for plotting. Please create one on a structured grid." ) - plt_fld = getattr(fld, field) - assert not (fld.pos is None or plt_fld is None) - + plt_fld = fld[field] norm = np.sqrt(plt_fld[0, :].T ** 2 + plt_fld[1, :].T ** 2) fig, ax = get_fig_ax(fig, ax) diff --git a/gstools/field/tools.py b/gstools/field/tools.py index fba20f94f..f46a0d305 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -70,11 +70,8 @@ def to_vtk_helper( fieldname : :class:`str`, optional Name of the field in the VTK file. Default: "field" """ + field = f_cls[field_select] if field_select in f_cls.field_names else None if f_cls.value_type == "vector": - if hasattr(f_cls, field_select): - field = getattr(f_cls, field_select) - else: - field = None if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): suf = ["_X", "_Y", "_Z"] fields = {} @@ -85,10 +82,6 @@ def to_vtk_helper( return vtk_export(filename, f_cls.pos, fields, f_cls.mesh_type) raise ValueError(f"Field.to_vtk: '{field_select}' not available.") if f_cls.value_type == "scalar": - if hasattr(f_cls, field_select): - field = getattr(f_cls, field_select) - else: - field = None if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): if filename is None: return to_vtk(f_cls.pos, {fieldname: field}, f_cls.mesh_type) From 08f4b2e4affc06b181a8b80c0e478bae05355f0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 18:29:28 +0200 Subject: [PATCH 16/49] Field: make get_store_config public --- gstools/field/base.py | 25 +++++++++++++++++++++++-- gstools/field/cond_srf.py | 4 ++-- gstools/field/srf.py | 2 +- gstools/krige/base.py | 2 +- gstools/transform/field.py | 2 +- 5 files changed, 28 insertions(+), 7 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 06af503a6..73c0a9a5a 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -141,7 +141,7 @@ def __call__( field : :class:`numpy.ndarray` the field values. """ - name, save = self._get_store_config(store) + name, save = self.get_store_config(store) pos, shape = self.pre_pos(pos, mesh_type) if field is None: field = np.zeros(shape, dtype=np.double) @@ -459,7 +459,28 @@ def set_pos(self, pos, mesh_type="unstructured"): self.pos = pos @staticmethod - def _get_store_config(store, default="field", fld_cnt=None): + def get_store_config(store, default="field", fld_cnt=None): + """ + Get sotrage configuration from given selection. + + Parameters + ---------- + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store fields (True/False) with default names + or with specified names. + The default is :any:`True` for default names. + default : :class:`str` or :class:`list`, optional + Default field names. The default is "field". + fld_cnt : :any:`None` or :class:`int`, optional + Number of fields when using lists. The default is None. + + Returns + ------- + name : :class:`str` or :class:`list` + Name(s) of field. + save : :class:`bool` or :class:`list` + Whether to save field(s). + """ # single field if fld_cnt is None: save = isinstance(store, str) or bool(store) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 04479c257..ab7ca1287 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -98,10 +98,10 @@ def __call__( field : :class:`numpy.ndarray` the conditioned SRF """ - name, save = self._get_store_config( + name, save = self.get_store_config( store, default=["field", "raw_field", "raw_krige"], fld_cnt=3 ) - krige_name, krige_save = self._get_store_config( + krige_name, krige_save = self.get_store_config( krige_store, default=["field", "krige_var"], fld_cnt=2 ) kwargs["mesh_type"] = mesh_type diff --git a/gstools/field/srf.py b/gstools/field/srf.py index c572ee8a8..d5e273f16 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -140,7 +140,7 @@ def __call__( field : :class:`numpy.ndarray` the SRF """ - name, save = self._get_store_config(store) + name, save = self.get_store_config(store) # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 40d387d47..1470954d0 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -224,7 +224,7 @@ def __call__( return_var &= not only_mean # don't return variance when calc. mean fld_cnt = 2 if return_var else 1 default = ["mean_field"] if only_mean else ["field", "krige_var"] - name, save = self._get_store_config( + name, save = self.get_store_config( store, default=default, fld_cnt=fld_cnt ) iso_pos, shape = self.pre_pos(pos, mesh_type) diff --git a/gstools/transform/field.py b/gstools/transform/field.py index ad8b32cb4..1fade9d95 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -217,7 +217,7 @@ def apply_function( if not callable(function): raise ValueError("transform.apply_function: function not a 'callable'") data = fld[field] - name, save = fld._get_store_config(store, default=field) + name, save = fld.get_store_config(store, default=field) if process: data = _pre_process(fld, data, keep_mean=keep_mean) data = function(data, **kwargs) From d7a2490d36b727207e722ea4395c5ecea300ce75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 19:27:19 +0200 Subject: [PATCH 17/49] Tests: reduce tolerance for an occasionally failing test --- tests/test_covmodel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 15123af20..f62fccaf3 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -260,7 +260,7 @@ def test_fitting(self): model.fit_variogram( self.gamma_x, self.gamma_y, weights=lambda x: 1 / (1 + x) ) - self.assertAlmostEqual(model.len_scale, len_save) + self.assertAlmostEqual(model.len_scale, len_save, places=6) # check ValueErrors with self.assertRaises(ValueError): model.fit_variogram(self.gamma_x, self.gamma_y, sill=2, var=3) From 6efc09ea61097129a3778aa2efddd14e323a7d64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 19:33:55 +0200 Subject: [PATCH 18/49] Examples: use new storage control for ensemble example --- examples/01_random_field/01_srf_ensemble.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/examples/01_random_field/01_srf_ensemble.py b/examples/01_random_field/01_srf_ensemble.py index 6174d4ce2..6f5398cc6 100644 --- a/examples/01_random_field/01_srf_ensemble.py +++ b/examples/01_random_field/01_srf_ensemble.py @@ -4,6 +4,8 @@ Creating an ensemble of random fields would also be a great idea. Let's reuse most of the previous code. + +We will set the position tuple `pos` before generation to reuse it afterwards. """ import numpy as np @@ -14,26 +16,26 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model) +srf.set_pos([x, y], "structured") ############################################################################### # This time, we did not provide a seed to :any:`SRF`, as the seeds will used # during the actual computation of the fields. We will create four ensemble -# members, for better visualisation and save them in a list and in a first +# members, for better visualisation, save them in to srf class and in a first # step, we will be using the loop counter as the seeds. ens_no = 4 -field = [] for i in range(ens_no): - field.append(srf.structured([x, y], seed=i)) + srf(seed=i, store=f"field{i}") ############################################################################### -# Now let's have a look at the results: +# Now let's have a look at the results. We can access the fields by name: fig, ax = pt.subplots(2, 2, sharex=True, sharey=True) ax = ax.flatten() for i in range(ens_no): - ax[i].imshow(field[i].T, origin="lower") + ax[i].imshow(srf[f"field{i}"].T, origin="lower") pt.show() ############################################################################### @@ -49,4 +51,4 @@ seed = MasterRNG(20170519) for i in range(ens_no): - field.append(srf.structured([x, y], seed=seed())) + srf(seed=seed(), store=f"better_field{i}") From 6db666ab4451041ab4b678deb16594430ba58def Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 19:45:03 +0200 Subject: [PATCH 19/49] Examples: use storage control in DWD example --- examples/08_geo_coordinates/01_dwd_krige.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index bac4b31cb..2ddc55ebd 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -130,8 +130,9 @@ def north_south_drift(lat, lon): g_lat = np.arange(47, 56.1, 0.1) g_lon = np.arange(5, 16.1, 0.1) -field, k_var = uk((g_lat, g_lon), mesh_type="structured") -mean = uk((g_lat, g_lon), mesh_type="structured", only_mean=True) +uk.set_pos((g_lat, g_lon), mesh_type="structured") +uk(return_var=False, store="temp_field") +uk(only_mean=True, store="mean_field") ############################################################################### # And that's it. Now let's have a look at the generated field and the input @@ -140,8 +141,8 @@ def north_south_drift(lat, lon): levels = np.linspace(5, 23, 64) fig, ax = plt.subplots(1, 3, figsize=[10, 5], sharey=True) sca = ax[0].scatter(lon, lat, c=temp, vmin=5, vmax=23, cmap="coolwarm") -co1 = ax[1].contourf(g_lon, g_lat, field, levels, cmap="coolwarm") -co2 = ax[2].contourf(g_lon, g_lat, mean, levels, cmap="coolwarm") +co1 = ax[1].contourf(g_lon, g_lat, uk["temp_field"], levels, cmap="coolwarm") +co2 = ax[2].contourf(g_lon, g_lat, uk["mean_field"], levels, cmap="coolwarm") [ax[i].plot(border[:, 0], border[:, 1], color="k") for i in range(3)] [ax[i].set_xlim([5, 16]) for i in range(3)] @@ -160,8 +161,8 @@ def north_south_drift(lat, lon): # a look at a cross-section at a longitude of 10 degree: fig, ax = plt.subplots() -ax.plot(g_lat, field[:, 50], label="Interpolated temperature") -ax.plot(g_lat, mean[:, 50], label="North-South mean drift") +ax.plot(g_lat, uk["temp_field"][:, 50], label="Interpolated temperature") +ax.plot(g_lat, uk["mean_field"][:, 50], label="North-South mean drift") ax.set_xlabel("Lat in deg") ax.set_ylabel("T in [°C]") ax.set_title("North-South cross-section at 10°") From d6aba851c6e0e206cf36774d1ccff418ce25a9c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 2 Aug 2021 19:53:15 +0200 Subject: [PATCH 20/49] Examples: use storage control in 2D cond ensemble --- .../01_2D_condition_ensemble.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py index 35d00bfd9..1ecd5f616 100644 --- a/examples/06_conditioned_fields/01_2D_condition_ensemble.py +++ b/examples/06_conditioned_fields/01_2D_condition_ensemble.py @@ -20,14 +20,14 @@ model = gs.Gaussian(dim=2, var=0.5, len_scale=5, anis=0.5, angles=-0.5) krige = gs.Krige(model, cond_pos=cond_pos, cond_val=cond_val) cond_srf = gs.CondSRF(krige) +cond_srf.set_pos([x, y], "structured") ############################################################################### # We create a list containing the generated conditioned fields. ens_no = 4 -field = [] for i in range(ens_no): - field.append(cond_srf.structured([x, y], seed=i)) + cond_srf(seed=i, store=f"f{i}") ############################################################################### # Now let's have a look at the pairwise differences between the generated @@ -35,19 +35,22 @@ fig, ax = plt.subplots(ens_no + 1, ens_no + 1, figsize=(8, 8)) # plotting kwargs for scatter and image -sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=np.max(field)) -im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=np.max(field)) +vmax = np.max([cond_srf[f"f{i}"] for i in range(ens_no)]) +sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax) +im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax) for i in range(ens_no): # conditioned fields and conditions - ax[i + 1, 0].imshow(field[i].T, **im_kwargs) + ax[i + 1, 0].imshow(cond_srf[f"f{i}"].T, **im_kwargs) ax[i + 1, 0].scatter(*cond_pos, **sc_kwargs) ax[i + 1, 0].set_ylabel(f"Field {i+1}", fontsize=10) - ax[0, i + 1].imshow(field[i].T, **im_kwargs) + ax[0, i + 1].imshow(cond_srf[f"f{i}"].T, **im_kwargs) ax[0, i + 1].scatter(*cond_pos, **sc_kwargs) ax[0, i + 1].set_title(f"Field {i+1}", fontsize=10) # absolute differences for j in range(ens_no): - ax[i + 1, j + 1].imshow(np.abs(field[i] - field[j]).T, **im_kwargs) + ax[i + 1, j + 1].imshow( + np.abs(cond_srf[f"f{i}"] - cond_srf[f"f{j}"]).T, **im_kwargs + ) # beautify plots ax[0, 0].axis("off") From bf94de84e523c39b41b82c5e8b009ba8d4c0a144 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 01:02:49 +0200 Subject: [PATCH 21/49] Field: allow indexing and slicing --- gstools/field/base.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 73c0a9a5a..fc02ecb45 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -104,10 +104,16 @@ def __init__( self.normalizer = normalizer self.trend = trend - def __getitem__(self, field): - if field in self.field_names: - return getattr(self, field) - raise ValueError(f"Field: requested field '{field}' not present") + def __getitem__(self, subscript): + if subscript in self.field_names: + return getattr(self, subscript) + if isinstance(subscript, int): + return self[self.field_names[subscript]] + if isinstance(subscript, slice): + return [self[f] for f in self.field_names[subscript]] + if isinstance(subscript, Iterable): + return [self[f] for f in subscript] + raise ValueError(f"Field: requested field '{subscript}' not present") def __call__( self, @@ -523,6 +529,11 @@ def pos(self, pos): if self.latlon: raise ValueError("Field: Vector fields not allowed for latlon") + @property + def all_fields(self): + """:class:`list`: All fields as stacked list.""" + return self[self.field_names] + @property def field_names(self): """:class:`list`: Names of present fields.""" From a02ddef5f6de3e2dabfeae1232c2e1532ff26a15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 01:03:57 +0200 Subject: [PATCH 22/49] Examples: use indexing in ensemble examples --- examples/01_random_field/01_srf_ensemble.py | 7 ++-- .../01_2D_condition_ensemble.py | 34 ++++++++++++------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/examples/01_random_field/01_srf_ensemble.py b/examples/01_random_field/01_srf_ensemble.py index 6f5398cc6..2563d9601 100644 --- a/examples/01_random_field/01_srf_ensemble.py +++ b/examples/01_random_field/01_srf_ensemble.py @@ -24,18 +24,18 @@ # members, for better visualisation, save them in to srf class and in a first # step, we will be using the loop counter as the seeds. - ens_no = 4 for i in range(ens_no): srf(seed=i, store=f"field{i}") ############################################################################### -# Now let's have a look at the results. We can access the fields by name: +# Now let's have a look at the results. We can access the fields by name or +# index: fig, ax = pt.subplots(2, 2, sharex=True, sharey=True) ax = ax.flatten() for i in range(ens_no): - ax[i].imshow(srf[f"field{i}"].T, origin="lower") + ax[i].imshow(srf[i].T, origin="lower") pt.show() ############################################################################### @@ -46,7 +46,6 @@ # provides a seed generator :any:`MasterRNG`. The loop, in which the fields are # generated would then look like - from gstools.random import MasterRNG seed = MasterRNG(20170519) diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py index 1ecd5f616..97c906f2d 100644 --- a/examples/06_conditioned_fields/01_2D_condition_ensemble.py +++ b/examples/06_conditioned_fields/01_2D_condition_ensemble.py @@ -24,10 +24,14 @@ ############################################################################### # We create a list containing the generated conditioned fields. +# By specifying ``store=[f"fld{i}", False, False]``, only the conditioned field +# is stored with the specified name. The raw random field and the raw kriging +# field is not stored. This way, we can access each conditioned field by index +# ``cond_srf[i]``: ens_no = 4 for i in range(ens_no): - cond_srf(seed=i, store=f"f{i}") + cond_srf(seed=i, store=[f"fld{i}", False, False]) ############################################################################### # Now let's have a look at the pairwise differences between the generated @@ -35,22 +39,20 @@ fig, ax = plt.subplots(ens_no + 1, ens_no + 1, figsize=(8, 8)) # plotting kwargs for scatter and image -vmax = np.max([cond_srf[f"f{i}"] for i in range(ens_no)]) -sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax) -im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax) +vmax = np.max(cond_srf.all_fields) +sc_kw = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax) +im_kw = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax) for i in range(ens_no): # conditioned fields and conditions - ax[i + 1, 0].imshow(cond_srf[f"f{i}"].T, **im_kwargs) - ax[i + 1, 0].scatter(*cond_pos, **sc_kwargs) - ax[i + 1, 0].set_ylabel(f"Field {i+1}", fontsize=10) - ax[0, i + 1].imshow(cond_srf[f"f{i}"].T, **im_kwargs) - ax[0, i + 1].scatter(*cond_pos, **sc_kwargs) - ax[0, i + 1].set_title(f"Field {i+1}", fontsize=10) + ax[i + 1, 0].imshow(cond_srf[i].T, **im_kw) + ax[i + 1, 0].scatter(*cond_pos, **sc_kw) + ax[i + 1, 0].set_ylabel(f"Field {i}", fontsize=10) + ax[0, i + 1].imshow(cond_srf[i].T, **im_kw) + ax[0, i + 1].scatter(*cond_pos, **sc_kw) + ax[0, i + 1].set_title(f"Field {i}", fontsize=10) # absolute differences for j in range(ens_no): - ax[i + 1, j + 1].imshow( - np.abs(cond_srf[f"f{i}"] - cond_srf[f"f{j}"]).T, **im_kwargs - ) + ax[i + 1, j + 1].imshow(np.abs(cond_srf[i] - cond_srf[j]).T, **im_kw) # beautify plots ax[0, 0].axis("off") @@ -59,3 +61,9 @@ a.set_xticks([]), a.set_yticks([]) fig.subplots_adjust(wspace=0, hspace=0) fig.show() + +############################################################################### +# To check if the generated fields are correct, we can have a look at their +# names: + +print(cond_srf.field_names) From 84bc379b163834bd25b9615bc2fc20b491f2c41e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 01:15:07 +0200 Subject: [PATCH 23/49] Krige: cleanup --- gstools/krige/base.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 1470954d0..bd6f00ed7 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -224,9 +224,8 @@ def __call__( return_var &= not only_mean # don't return variance when calc. mean fld_cnt = 2 if return_var else 1 default = ["mean_field"] if only_mean else ["field", "krige_var"] - name, save = self.get_store_config( - store, default=default, fld_cnt=fld_cnt - ) + name, save = self.get_store_config(store, default, fld_cnt) + iso_pos, shape = self.pre_pos(pos, mesh_type) pnt_cnt = len(iso_pos[0]) From 41f1fa6a9a3396a34bfd1519b0e4813c034d1569 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 12:51:42 +0200 Subject: [PATCH 24/49] Field: bigfix in getitem; add delitem; apply shape in post_field --- gstools/field/base.py | 53 +++++++++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index fc02ecb45..bf4439804 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -104,16 +104,38 @@ def __init__( self.normalizer = normalizer self.trend = trend - def __getitem__(self, subscript): - if subscript in self.field_names: - return getattr(self, subscript) - if isinstance(subscript, int): - return self[self.field_names[subscript]] - if isinstance(subscript, slice): - return [self[f] for f in self.field_names[subscript]] - if isinstance(subscript, Iterable): - return [self[f] for f in subscript] - raise ValueError(f"Field: requested field '{subscript}' not present") + def __getitem__(self, key): + if key in self.field_names: + return getattr(self, key) + if isinstance(key, int): + return self[self.field_names[key]] + if isinstance(key, slice): + return [self[f] for f in self.field_names[key]] + if isinstance(key, Iterable) and not isinstance(key, str): + return [self[f] for f in key] + raise KeyError(f"{self.name}: requested field '{key}' not present") + + def __delitem__(self, key): + names = [] + if key in self.field_names: + names = [key] + elif isinstance(key, int): + names = [self.field_names[key]] + elif isinstance(key, slice): + names = self.field_names[key] + elif isinstance(key, Iterable) and not isinstance(key, str): + for k in key: + k = self.field_names[k] if isinstance(key, int) else k + names.append(k) + else: + raise KeyError(f"{self.name}: requested field '{key}' not present") + for name in names: + if name not in self.field_names: + raise KeyError( + f"{self.name}: requested field '{name}' not present" + ) + delattr(self, name) + del self._field_names[self._field_names.index(name)] def __call__( self, @@ -291,6 +313,9 @@ def post_field(self, field, name="field", process=True, save=True): field : :class:`numpy.ndarray` Processed field values. """ + if self.field_shape is None: + raise ValueError("post_field: no 'field_shape' present.") + field = np.asarray(field, dtype=np.double).reshape(self.field_shape) if process: if self.pos is None: raise ValueError("post_field: no 'pos' tuple set for field.") @@ -319,11 +344,9 @@ def post_field(self, field, name="field", process=True, save=True): setattr(self, name, field) return field - def delete_fields(self): - """Delete all present fields""" - for field in self._field_names: - delattr(self, field) - self._field_names = [] + def delete_fields(self, select=None): + """Delete selected fields.""" + del self[self.field_names if select is None else select] def transform( self, method, field="field", store=True, process=False, **kwargs From 4807ef653c407bd0461266babf2271eae931a33f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 12:52:31 +0200 Subject: [PATCH 25/49] CondSRF: optimize storage if reusing fields --- gstools/field/cond_srf.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index ab7ca1287..140abe3ca 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -124,17 +124,20 @@ def __call__( and name[2] in self.field_names and krige_name[1] in self.krige.field_names ): + reuse = True rawkrige, krige_var = self[name[2]], self.krige[krige_name[1]] else: + reuse = False rawkrige, krige_var = self.krige(**kwargs) var_scale, nugget = self.get_scaling(krige_var, shape) - # need to use a copy to not alter "field" by reference - # store krige field - self.krige.post_field( - rawkrige.copy(), krige_name[0], post_process, krige_save[0] - ) + # store krige field (need a copy to not alter field by reference) + if not reuse or krige_name[0] not in self.krige.field_names: + self.krige.post_field( + rawkrige.copy(), krige_name[0], post_process, krige_save[0] + ) # store raw krige field - self.post_field(rawkrige.copy(), name[2], False, save[2]) + if not reuse: + self.post_field(rawkrige, name[2], False, save[2]) # store raw random field self.post_field(rawfield, name[1], False, save[1]) # store cond random field From 03e28b278d8b2b615177dade0e978dcbcaf9bc5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 13:24:26 +0200 Subject: [PATCH 26/49] Refactor: prevent unnecessary array copies with np.asarray --- gstools/covmodel/base.py | 10 +++++----- gstools/covmodel/fit.py | 8 ++++---- gstools/covmodel/models.py | 28 ++++++++++++++-------------- gstools/covmodel/tools.py | 8 ++++---- gstools/field/base.py | 6 +++--- gstools/field/generator.py | 4 ++-- gstools/krige/base.py | 8 +++++--- gstools/krige/tools.py | 4 ++-- gstools/normalizer/base.py | 4 ++-- gstools/normalizer/tools.py | 4 ++-- gstools/tools/geometric.py | 32 ++++++++++++++++---------------- gstools/tools/misc.py | 6 +++--- gstools/tools/special.py | 8 ++++---- gstools/transform/field.py | 6 +++--- gstools/variogram/binning.py | 4 ++-- gstools/variogram/variogram.py | 8 ++++---- 16 files changed, 75 insertions(+), 73 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index f8c86589c..76de7d197 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -265,7 +265,7 @@ def cor_spatial(self, pos): def vario_nugget(self, r): """Isotropic variogram of the model respecting the nugget at r=0.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.variogram(r[r_gz]) @@ -274,7 +274,7 @@ def vario_nugget(self, r): def cov_nugget(self, r): """Isotropic covariance of the model respecting the nugget at r=0.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.covariance(r[r_gz]) @@ -501,7 +501,7 @@ def spectral_density(self, k): k : :class:`float` Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert` """ - k = np.array(np.abs(k), dtype=np.double) + k = np.asarray(np.abs(k), dtype=np.double) return self._sft.transform(self.correlation, k, ret_err=False) def spectral_rad_pdf(self, r): @@ -525,14 +525,14 @@ def _has_ppf(self): def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" - pos = np.array(pos, dtype=np.double).reshape((self.field_dim, -1)) + pos = np.asarray(pos, dtype=np.double).reshape((self.field_dim, -1)) if self.latlon: return latlon2pos(pos) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) def anisometrize(self, pos): """Bring a position tuple into the anisotropic coordinate-system.""" - pos = np.array(pos, dtype=np.double).reshape((self.dim, -1)) + pos = np.asarray(pos, dtype=np.double).reshape((self.dim, -1)) if self.latlon: return pos2latlon(pos) return np.dot( diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index c67c9839b..8629cdffa 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -321,8 +321,8 @@ def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0): def _check_vario(model, x_data, y_data): # prepare variogram data - x_data = np.array(x_data).reshape(-1) - y_data = np.array(y_data).reshape(-1) + x_data = np.asarray(x_data).reshape(-1) + y_data = np.asarray(y_data).reshape(-1) # if multiple variograms are given, they will be interpreted # as directional variograms along the main rotated axes of the model is_dir_vario = False @@ -355,9 +355,9 @@ def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario): elif is_dir_vario: if weights.size * model.dim == x_data.size: weights = np.tile(weights, model.dim) - weights = 1.0 / np.array(weights).reshape(-1) + weights = 1.0 / np.asarray(weights).reshape(-1) else: - weights = 1.0 / np.array(weights).reshape(-1) + weights = 1.0 / np.asarray(weights).reshape(-1) curve_fit_kwargs["sigma"] = weights curve_fit_kwargs["absolute_sigma"] = True diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index c05a92dcb..9c770ea59 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -73,14 +73,14 @@ def default_rescale(self): return np.sqrt(np.pi) / 2.0 def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) return (self.len_rescaled / 2.0 / np.sqrt(np.pi)) ** self.dim * np.exp( -((k * self.len_rescaled / 2.0) ** 2) ) def spectral_rad_cdf(self, r): """Gaussian radial spectral cdf.""" - r = np.array(r, dtype=np.double) + r = np.asarray(r, dtype=np.double) if self.dim == 1: return sps.erf(r * self.len_rescaled / 2.0) if self.dim == 2: @@ -100,7 +100,7 @@ def spectral_rad_ppf(self, u): ----- Not defined for 3D. """ - u = np.array(u, dtype=np.double) + u = np.asarray(u, dtype=np.double) if self.dim == 1: return 2.0 / self.len_rescaled * sps.erfinv(u) if self.dim == 2: @@ -143,7 +143,7 @@ def cor(self, h): return np.exp(-h) def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) return ( self.len_rescaled ** self.dim * sps.gamma((self.dim + 1) / 2.0) @@ -153,7 +153,7 @@ def spectral_density(self, k): # noqa: D102 def spectral_rad_cdf(self, r): """Exponential radial spectral cdf.""" - r = np.array(r, dtype=np.double) + r = np.asarray(r, dtype=np.double) if self.dim == 1: return np.arctan(r * self.len_rescaled) * 2.0 / np.pi if self.dim == 2: @@ -178,7 +178,7 @@ def spectral_rad_ppf(self, u): ----- Not defined for 3D. """ - u = np.array(u, dtype=np.double) + u = np.asarray(u, dtype=np.double) if self.dim == 1: return np.tan(np.pi / 2 * u) / self.len_rescaled if self.dim == 2: @@ -341,7 +341,7 @@ def default_opt_arg_bounds(self): def cor(self, h): """Matérn normalized correlation function.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) # for nu > 20 we just use the gaussian model if self.nu > 20.0: return np.exp(-((h / 2.0) ** 2)) @@ -360,7 +360,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) # for nu > 20 we just use an approximation of the gaussian model if self.nu > 20.0: return ( @@ -570,7 +570,7 @@ class Circular(CovModel): def cor(self, h): """Circular normalized correlation function.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) res = np.zeros_like(h) # arccos is instable around h=1 h_l1 = h < 1.0 @@ -661,7 +661,7 @@ class HyperSpherical(CovModel): def cor(self, h): """Hyper-Spherical normalized correlation function.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) res = np.zeros_like(h) h_l1 = h < 1 nu = (self.dim - 1.0) / 2.0 @@ -670,7 +670,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) res = np.empty_like(k) kl = k * self.len_rescaled kl_gz = np.logical_not(np.isclose(k, 0)) @@ -751,7 +751,7 @@ def default_opt_arg_bounds(self): def cor(self, h): """Super-Spherical normalized correlation function.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) res = np.zeros_like(h) h_l1 = h < 1 fac = 1.0 / sps.hyp2f1(0.5, -self.nu, 1.5, 1.0) @@ -846,7 +846,7 @@ def check_opt_arg(self): def cor(self, h): """J-Bessel correlation.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) h_gz = np.logical_not(np.isclose(h, 0)) hh = h[h_gz] res = np.ones_like(h) @@ -855,7 +855,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) k_ll = k < 1.0 / self.len_rescaled kk = k[k_ll] res = np.zeros_like(k) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index b3abf1ab2..56e65ec0e 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -71,12 +71,12 @@ def correlation(self, r): def correlation_from_cor(self, r): """Correlation function of the model.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) return self.cor(r / self.len_rescaled) def cor_from_correlation(self, h): """Correlation taking a non-dimensional range.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) return self.correlation(h * self.len_rescaled) abstract = True @@ -270,7 +270,7 @@ def check_arg_in_bounds(model, arg, val=None): raise ValueError("check bounds: unknown argument: {}".format(arg)) bnd = list(model.arg_bounds[arg]) val = getattr(model, arg) if val is None else val - val = np.array(val) + val = np.asarray(val) error_case = 0 if len(bnd) == 2: bnd.append("cc") # use closed intervals by default @@ -332,7 +332,7 @@ def spectral_rad_pdf(model, r): PDF values. """ - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) if model.dim > 1: r_gz = np.logical_not(np.isclose(r, 0)) # to prevent numerical errors, we just calculate where r>0 diff --git a/gstools/field/base.py b/gstools/field/base.py index bf4439804..590c359ed 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -47,7 +47,7 @@ def _pos_equal(pos1, pos2): def _set_mean_trend(value, dim): if callable(value) or value is None: return value - value = np.array(value, dtype=np.double).ravel() + value = np.asarray(value, dtype=np.double).ravel() if value.size > 1 and value.size != dim: # vector mean raise ValueError(f"Mean/Trend: Wrong size ({value})") return value if value.size > 1 else value.item() @@ -174,7 +174,7 @@ def __call__( if field is None: field = np.zeros(shape, dtype=np.double) else: - field = np.array(field, dtype=np.double).reshape(shape) + field = np.asarray(field, dtype=np.double).reshape(shape) return self.post_field(field, name, post_process, save) def structured(self, *args, **kwargs): @@ -542,7 +542,7 @@ def pos(self, pos): if self.mesh_type is None: raise ValueError("Field.pos: can't set 'pos' without 'mesh_type'") if self.mesh_type == "unstructured": - self._pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) + self._pos = np.asarray(pos, dtype=np.double).reshape(self.dim, -1) self._field_shape = np.shape(self._pos[0]) else: self._pos, self._field_shape = format_struct_pos_dim(pos, self.dim) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 9b3eeca67..9c2761bf6 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -125,7 +125,7 @@ def __call__(self, pos, add_nugget=True): :class:`numpy.ndarray` the random modes """ - pos = np.array(pos, dtype=np.double) + pos = np.asarray(pos, dtype=np.double) summed_modes = summate(self._cov_sample, self._z_1, self._z_2, pos) nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0 return np.sqrt(self.model.var / self._mode_no) * summed_modes + nugget @@ -415,7 +415,7 @@ def __call__(self, pos): :class:`numpy.ndarray` the random modes """ - pos = np.array(pos, dtype=np.double) + pos = np.asarray(pos, dtype=np.double) summed_modes = summate_incompr( self._cov_sample, self._z_1, self._z_2, pos ) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index bd6f00ed7..b823a5f81 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -368,7 +368,9 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): the drift values at the given positions """ if ext_drift is not None: - ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) + ext_drift = np.array( + ext_drift, dtype=np.double, ndmin=2, copy=False + ) if ext_drift.size == 0: # treat empty array as no ext_drift return np.array([]) if set_cond: @@ -381,7 +383,7 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): raise ValueError("Krige: wrong number of external drifts.") if np.prod(ext_shape) != np.prod(shape): raise ValueError("Krige: wrong number of ext. drift values.") - return np.array(ext_drift, dtype=np.double).reshape(shape) + return np.asarray(ext_drift, dtype=np.double).reshape(shape) if not set_cond and self._cond_ext_drift.size > 0: raise ValueError("Krige: wrong number of ext. drift values.") return np.array([]) @@ -607,7 +609,7 @@ def cond_err(self, value): "krige.cond_err: measurement errors can't be given, " "when interpolator should be exact." ) - value = np.array(value, dtype=np.double).reshape(-1) + value = np.asarray(value, dtype=np.double).reshape(-1) if value.size == 1: self._cond_err = value.item() else: diff --git a/gstools/krige/tools.py b/gstools/krige/tools.py index 5d220f034..d84c44855 100644 --- a/gstools/krige/tools.py +++ b/gstools/krige/tools.py @@ -44,8 +44,8 @@ def set_condition(cond_pos, cond_val, dim): the error checked cond_val """ # convert the input for right shapes and dimension checks - cond_val = np.array(cond_val, dtype=np.double).reshape(-1) - cond_pos = np.array(cond_pos, dtype=np.double).reshape(dim, -1) + cond_val = np.asarray(cond_val, dtype=np.double).reshape(-1) + cond_pos = np.asarray(cond_pos, dtype=np.double).reshape(dim, -1) if len(cond_pos[0]) != len(cond_val): raise ValueError( "Please check your 'cond_pos' and 'cond_val' parameters. " diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index da6e0a9cb..a9fcf9b0c 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -67,10 +67,10 @@ def _kernel_loglikelihood(self, data): return res + np.sum(np.log(np.maximum(1e-16, self._derivative(data)))) def _check_input(self, data, data_range=None, return_output_template=True): - is_data = np.array(np.logical_not(np.isnan(data))) + is_data = np.logical_not(np.isnan(data)) if return_output_template: out = np.full_like(data, np.nan, dtype=np.double) - data = np.array(data, dtype=np.double)[is_data] + data = np.asarray(data, dtype=np.double)[is_data] if data_range is not None and np.min(np.abs(data_range)) < np.inf: dat_in = np.logical_and(data > data_range[0], data < data_range[1]) if not np.all(dat_in): diff --git a/gstools/normalizer/tools.py b/gstools/normalizer/tools.py index f71b758f0..0ec10f773 100644 --- a/gstools/normalizer/tools.py +++ b/gstools/normalizer/tools.py @@ -90,7 +90,7 @@ def apply_mean_norm_trend( pos, shape, dim = format_unstruct_pos_shape( pos, field.shape, check_stacked_shape=stacked ) - field = np.array(field, dtype=np.double).reshape(shape) + field = np.asarray(field, dtype=np.double).reshape(shape) else: dim = len(pos) if not stacked: @@ -169,7 +169,7 @@ def remove_trend_norm_mean( pos, shape, dim = format_unstruct_pos_shape( pos, field.shape, check_stacked_shape=stacked ) - field = np.array(field, dtype=np.double).reshape(shape) + field = np.asarray(field, dtype=np.double).reshape(shape) else: dim = len(pos) if not stacked: diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 204df5aea..843114479 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -79,7 +79,7 @@ def set_angles(dim, angles): ----- If too few angles are given, they are filled up with `0`. """ - out_angles = np.array(angles, dtype=np.double) + out_angles = np.asarray(angles, dtype=np.double) out_angles = np.atleast_1d(out_angles)[: no_of_angles(dim)] # fill up the rotation angle array with zeros out_angles = np.pad( @@ -110,7 +110,7 @@ def set_anis(dim, anis): ----- If too few anisotropy ratios are given, they are filled up with `1`. """ - out_anis = np.array(anis, dtype=np.double) + out_anis = np.asarray(anis, dtype=np.double) out_anis = np.atleast_1d(out_anis)[: dim - 1] if len(out_anis) < dim - 1: # fill up the anisotropies with ones, such that len()==dim-1 @@ -351,9 +351,9 @@ def generate_grid(pos): :class:`numpy.ndarray` Unstructured position tuple. """ - return np.array(np.meshgrid(*pos, indexing="ij"), dtype=np.double).reshape( - (len(pos), -1) - ) + return np.asarray( + np.meshgrid(*pos, indexing="ij"), dtype=np.double + ).reshape((len(pos), -1)) def generate_st_grid(pos, time, mesh_type="unstructured"): @@ -379,14 +379,14 @@ def generate_st_grid(pos, time, mesh_type="unstructured"): ----- Time dimension will be the last one. """ - time = np.array(time, dtype=np.double).reshape(-1) + time = np.asarray(time, dtype=np.double).reshape(-1) if mesh_type != "unstructured": pos = generate_grid(pos) else: - pos = np.array(pos, dtype=np.double, ndmin=2) + pos = np.array(pos, dtype=np.double, ndmin=2, copy=False) out = [np.repeat(p.reshape(-1), np.size(time)) for p in pos] out.append(np.tile(time, np.size(pos[0]))) - return np.array(out, dtype=np.double) + return np.asarray(out, dtype=np.double) # conversion ################################################################## @@ -416,12 +416,12 @@ def format_struct_pos_dim(pos, dim): Shape of the resulting field. """ if dim == 1: - pos = (np.array(pos, dtype=np.double).reshape(-1),) + pos = (np.asarray(pos, dtype=np.double).reshape(-1),) elif len(pos) != dim: raise ValueError("Formatting: position tuple doesn't match dimension.") else: - pos = tuple(np.array(p_i, dtype=np.double).reshape(-1) for p_i in pos) - shape = tuple(len(p_i) for p_i in pos) + pos = tuple(np.asarray(p, dtype=np.double).reshape(-1) for p in pos) + shape = tuple(len(p) for p in pos) return pos, shape @@ -479,7 +479,7 @@ def format_struct_pos_shape(pos, shape, check_stacked_shape=False): else: wrong_shape = True else: - struct_size = np.prod([p_i.size for p_i in check_pos]) + struct_size = np.prod([p.size for p in check_pos]) # case: 1D unstacked if check_pos.size == shape_size: dim = 1 @@ -552,7 +552,7 @@ def format_unstruct_pos_shape(pos, shape, check_stacked_shape=False): # now we try to be smart pre_len = len(np.atleast_1d(pos)) # care about 1D: pos can be given as 1D array here -> convert to 2D array - pos = np.array(pos, dtype=np.double, ndmin=2) + pos = np.array(pos, dtype=np.double, ndmin=2, copy=False) post_len = len(pos) # first array dimension should be spatial dimension (1D is special case) dim = post_len if pre_len == post_len else 1 @@ -606,7 +606,7 @@ def ang2dir(angles, dtype=np.double, dim=None): the array of direction vectors """ pre_dim = np.asanyarray(angles).ndim - angles = np.array(angles, ndmin=2, dtype=dtype) + angles = np.array(angles, ndmin=2, dtype=dtype, copy=False) if len(angles.shape) > 2: raise ValueError(f"Can't interpret angles array {angles}") dim = angles.shape[1] + 1 if dim is None else dim @@ -644,7 +644,7 @@ def latlon2pos(latlon, radius=1.0, dtype=np.double): :class:`numpy.ndarray` the 3D position array """ - latlon = np.array(latlon, dtype=dtype).reshape((2, -1)) + latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1)) lat, lon = np.deg2rad(latlon) return np.array( ( @@ -675,7 +675,7 @@ def pos2latlon(pos, radius=1.0, dtype=np.double): :class:`numpy.ndarray` the 3D position array """ - pos = np.array(pos, dtype=dtype).reshape((3, -1)) + pos = np.asarray(pos, dtype=dtype).reshape((3, -1)) # prevent numerical errors in arcsin lat = np.arcsin(np.maximum(np.minimum(pos[2] / radius, 1.0), -1.0)) lon = np.arctan2(pos[1], pos[0]) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 53e3c9cd6..c3e994684 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -104,7 +104,7 @@ def eval_func( # care about scalar inputs func_val = 0 if func_val is None else func_val if broadcast and not callable(func_val) and np.size(func_val) == 1: - return np.array(func_val, dtype=np.double).item() + return np.asarray(func_val, dtype=np.double).item() if not callable(func_val): func_val = _func_from_single_val(func_val, dim, value_type=value_type) # care about mesh and function call @@ -112,7 +112,7 @@ def eval_func( pos, shape = format_struct_pos_dim(pos, dim) pos = generate_grid(pos) else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) + pos = np.asarray(pos, dtype=np.double).reshape(dim, -1) shape = np.shape(pos[0]) # prepend dimension if we have a vector field if value_type == "vector": @@ -125,7 +125,7 @@ def _func_from_single_val(value, dim=None, value_type="scalar"): v_d = dim if value_type == "vector" else 1 # value dim if v_d is None: # pragma: no cover raise ValueError("_func_from_single_val: dim needed for vector value.") - value = np.array(value, dtype=np.double).ravel()[:v_d] + value = np.asarray(value, dtype=np.double).ravel()[:v_d] # fill up vector valued output to dimension with last value value = np.pad( value, (0, v_d - len(value)), "constant", constant_values=value[-1] diff --git a/gstools/tools/special.py b/gstools/tools/special.py index 63940c08f..b6990ec89 100644 --- a/gstools/tools/special.py +++ b/gstools/tools/special.py @@ -86,7 +86,7 @@ def exp_int(s, x): return sps.exp1(x) if np.isclose(s, np.around(s)) and s > -0.5: return sps.expn(int(np.around(s)), x) - x = np.array(x, dtype=np.double) + x = np.asarray(x, dtype=np.double) x_neg = x < 0 x = np.abs(x) x_compare = x ** min((10, max(((1 - s), 1)))) @@ -146,7 +146,7 @@ def tplstable_cor(r, len_scale, hurst, alpha): alpha : :class:`float`, optional Shape parameter of the stable model. """ - r = np.array(np.abs(r / len_scale), dtype=np.double) + r = np.asarray(np.abs(r / len_scale), dtype=np.double) r[np.isclose(r, 0)] = 0 # hack to prevent numerical errors res = np.ones_like(r) res[r > 0] = (2 * hurst / alpha) * exp_int( @@ -179,7 +179,7 @@ def tpl_exp_spec_dens(k, dim, len_scale, hurst, len_low=0.0): spectal density of the TPLExponential model """ if np.isclose(len_low, 0.0): - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) z = (k * len_scale) ** 2 a = hurst + dim / 2.0 b = hurst + 0.5 @@ -218,7 +218,7 @@ def tpl_gau_spec_dens(k, dim, len_scale, hurst, len_low=0.0): spectal density of the TPLExponential model """ if np.isclose(len_low, 0.0): - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) z = np.array((k * len_scale / 2.0) ** 2) res = np.empty_like(z) z_gz = z > 0.1 # greater zero diff --git a/gstools/transform/field.py b/gstools/transform/field.py index 1fade9d95..90ee6c6ec 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -673,7 +673,7 @@ def array_discrete( elif thresholds == "equal": mean = np.mean(field) if mean is None else float(mean) var = np.var(field) if var is None else float(var) - values = np.array(values) + values = np.asarray(values) n = len(values) p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] rescale = np.sqrt(var * 2) @@ -684,8 +684,8 @@ def array_discrete( raise ValueError( "discrete transformation: len(values) != len(thresholds) + 1" ) - values = np.array(values) - thresholds = np.array(thresholds) + values = np.asarray(values) + thresholds = np.asarray(thresholds) # check thresholds if not np.all(thresholds[:-1] < thresholds[1:]): raise ValueError( diff --git a/gstools/variogram/binning.py b/gstools/variogram/binning.py index c7be0eca0..c2ea7a6c2 100644 --- a/gstools/variogram/binning.py +++ b/gstools/variogram/binning.py @@ -80,13 +80,13 @@ def standard_bins( if mesh_type != "unstructured": pos = generate_grid(format_struct_pos_dim(pos, dim)[0]) else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) + pos = np.asarray(pos, dtype=np.double).reshape(dim, -1) pos = latlon2pos(pos) if latlon else pos pnt_cnt = len(pos[0]) box = [] for axis in pos: box.append([np.min(axis), np.max(axis)]) - box = np.array(box) + box = np.asarray(box) diam = np.linalg.norm(box[:, 0] - box[:, 1]) # convert diameter to great-circle distance if using latlon diam = chordal_to_great_circle(diam) if latlon else diam diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index c6d98e9fc..b25dacb47 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -239,7 +239,7 @@ def vario_estimate( John Wiley & Sons. (2007) """ if bin_edges is not None: - bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) + bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 # allow multiple fields at same positions (ndmin=2: first axis -> field ID) # need to convert to ma.array, since list of ma.array is not recognised @@ -292,7 +292,7 @@ def vario_estimate( # set directions dir_no = 0 if direction is not None and dim > 1: - direction = np.array(direction, ndmin=2, dtype=np.double) + direction = np.array(direction, ndmin=2, dtype=np.double, copy=False) if len(direction.shape) > 2: raise ValueError(f"Can't interpret directions: {direction}") if direction.shape[1] != dim: @@ -441,9 +441,9 @@ def vario_estimate_axis( field = np.ma.array(field, ndmin=1, dtype=np.double) if missing: field.mask = np.logical_or(field.mask, missing_mask) - mask = np.array(np.ma.getmaskarray(field), dtype=np.int32) + mask = np.asarray(np.ma.getmaskarray(field), dtype=np.int32) else: - field = np.array(field, ndmin=1, dtype=np.double) + field = np.array(field, ndmin=1, dtype=np.double, copy=False) missing_mask = None # free space axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction) From aa579b9ed5d5bf158a206c8675903b87189e0e0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 14:56:08 +0200 Subject: [PATCH 27/49] Changlog: update --- CHANGELOG.md | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e6d92ace..b7a2239b8 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,31 @@ All notable changes to **GSTools** will be documented in this file. + +## [1.3.3] - Pure Pink - ? + +### Enhancements +See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197) +- `gstools.transform`: + - add keywords `field`, `store` and `process` keyword to all transformations to control storage and respect `normalizer` + - added `apply_function` transformation + - added `apply` as wrapper for all transformations + - added `transform` method to all `Field` (sub)classes as interface to `transform.apply` + - added checks for normal fields to work smoothly with recently added `normalizer` submodule +- `Field`: + - allow naming fields when generating and control storage with `store` keyword + - all subclasses now have the `post_process` keyword (apply mean, normalizer, trend) + - added subscription to access fields by name (`Field["field"]`) + - added `set_pos` method to set position tuple + - allow reusing present `pos` tuple + - added `pos`, `mesh_type`, `field_names`, `field_shape`, `all_fields` properties +- `CondSRF`: + - memory optimization by forwarding `pos` from underlying `krige` instance + - only recalculate kriging field if `pos` tuple changed (optimized ensemble generation) +- performance improvement by using `np.asarray` instead of `np.array` where possible +- updated examples to use new features + + ## [1.3.2] - Pure Pink - 2021-07 ### Bugfixes @@ -280,7 +305,7 @@ All notable changes to **GSTools** will be documented in this file. First release of GSTools. -[Unreleased]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD +[1.3.3]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD [1.3.2]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.1...v1.3.2 [1.3.1]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.0...v1.3.1 [1.3.0]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.1...v1.3.0 From dca926bca288e1145eb797443f81ec9bb2d9f876 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 15:07:16 +0200 Subject: [PATCH 28/49] CovModel.fit: minor refactor --- gstools/covmodel/fit.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index 8629cdffa..76dfb1375 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -352,11 +352,9 @@ def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario): weights = 1.0 / weights(x_data) elif isinstance(weights, str) and weights == "inv": weights = 1.0 + x_data - elif is_dir_vario: - if weights.size * model.dim == x_data.size: - weights = np.tile(weights, model.dim) - weights = 1.0 / np.asarray(weights).reshape(-1) else: + if is_dir_vario and weights.size * model.dim == x_data.size: + weights = np.tile(weights, model.dim) weights = 1.0 / np.asarray(weights).reshape(-1) curve_fit_kwargs["sigma"] = weights curve_fit_kwargs["absolute_sigma"] = True From fc7e52941e4009336e6ca15b67a9ad4ce3278e95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 16:30:05 +0200 Subject: [PATCH 29/49] CondSRF: fix subtle bug when setting pos tuple --- gstools/field/base.py | 31 +++++++++++++++++++++---------- gstools/field/cond_srf.py | 31 +++++++++++++++++++++++++++---- 2 files changed, 48 insertions(+), 14 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 590c359ed..332589dfb 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -262,7 +262,7 @@ def pre_pos(self, pos=None, mesh_type="unstructured", info=False): Isometrized position tuple. shape : :class:`tuple` Shape of the resulting field. - info : :class:`dict` + info : :class:`dict`, optional Information about settings. """ info_ret = {"deleted": False} @@ -272,14 +272,7 @@ def pre_pos(self, pos=None, mesh_type="unstructured", info=False): if self.mesh_type is None: raise ValueError("Field: no 'mesh_type' present") else: - old_pos = copy(self.pos) - # save pos and mesh-type - self.set_pos(pos, mesh_type) - # remove present fields if new pos is different from current - if not _pos_equal(old_pos, self.pos): - self.delete_fields() - info_ret["deleted"] = True - del old_pos + info_ret = self.set_pos(pos, mesh_type, info=True) if self.mesh_type != "unstructured": pos = generate_grid(self.pos) else: @@ -471,7 +464,7 @@ def plot( return r - def set_pos(self, pos, mesh_type="unstructured"): + def set_pos(self, pos, mesh_type="unstructured", info=False): """ Set positions and mesh_type. @@ -483,9 +476,27 @@ def set_pos(self, pos, mesh_type="unstructured"): mesh_type : :class:`str`, optional 'structured' / 'unstructured' Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information + + Returns + ------- + info : :class:`dict` + Information about settings. """ + info_ret = {"deleted": False} + old_type = copy(self.mesh_type) + old_pos = copy(self.pos) + # save pos and mesh-type self.mesh_type = mesh_type self.pos = pos + # remove present fields if new pos is different from current + if old_type != self.mesh_type or not _pos_equal(old_pos, self.pos): + self.delete_fields() + info_ret["deleted"] = True + del old_pos + if info: + return info_ret @staticmethod def get_store_config(store, default="field", fld_cnt=None): diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 140abe3ca..f06729f73 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -112,10 +112,7 @@ def __call__( # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape - iso_pos, shape, info = self.krige.pre_pos(pos, mesh_type, info=True) - # also reset fields in this class (not only krige with call above) - if info["deleted"]: - self.delete_fields() + iso_pos, shape, info = self.pre_pos(pos, mesh_type, info=True) # generate the field rawfield = np.reshape(self.generator(iso_pos, add_nugget=False), shape) # call krige on already set pos (reuse already calculated fields) @@ -194,6 +191,32 @@ def set_generator(self, generator, **generator_kwargs): else: raise ValueError(f"gstools.CondSRF: Unknown generator {generator}") + def set_pos(self, pos, mesh_type="unstructured", info=False): + """ + Set positions and mesh_type. + + Parameters + ---------- + pos : :any:`iterable` + the position tuple, containing main direction and transversal + directions + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information + + Returns + ------- + info : :class:`dict` + Information about settings. + """ + info_ret = super().set_pos(pos, mesh_type, info=True) + if info_ret["deleted"]: + self.krige.delete_fields() + if info: + return info_ret + @property def pos(self): """:class:`tuple`: The position tuple of the field.""" From c76fd6d6b0295c96837b878aba809d7abcf3010f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 17:13:44 +0200 Subject: [PATCH 30/49] Field: add default_names attr; add __len__ and __contains__ magic methods --- gstools/field/base.py | 17 +++++++++++++++-- gstools/field/cond_srf.py | 11 ++++++----- gstools/krige/base.py | 5 ++++- 3 files changed, 25 insertions(+), 8 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 332589dfb..11a3350df 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -77,6 +77,9 @@ class Field: Dimension of the field if no model is given. """ + default_names = ["field"] + """:class:`list`: Default field names.""" + def __init__( self, model=None, @@ -104,6 +107,12 @@ def __init__( self.normalizer = normalizer self.trend = trend + def __len__(self): + return len(self.field_names) + + def __contains__(self, item): + return item in self.field_names + def __getitem__(self, key): if key in self.field_names: return getattr(self, key) @@ -498,8 +507,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): if info: return info_ret - @staticmethod - def get_store_config(store, default="field", fld_cnt=None): + def get_store_config(self, store, default=None, fld_cnt=None): """ Get sotrage configuration from given selection. @@ -521,6 +529,11 @@ def get_store_config(store, default="field", fld_cnt=None): save : :class:`bool` or :class:`list` Whether to save field(s). """ + if default is None: + if fld_cnt is None: + default = self.default_names[0] + else: + default = self.default_names # single field if fld_cnt is None: save = isinstance(store, str) or bool(store) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index f06729f73..4e24947b6 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -43,6 +43,9 @@ class CondSRF(Field): Have a look at the provided generators for further information. """ + default_names = ["field", "raw_field", "raw_krige"] + """:class:`list`: Default field names.""" + def __init__(self, krige, generator="RandMeth", **generator_kwargs): if not isinstance(krige, Krige): raise ValueError("CondSRF: krige should be an instance of Krige.") @@ -98,11 +101,9 @@ def __call__( field : :class:`numpy.ndarray` the conditioned SRF """ - name, save = self.get_store_config( - store, default=["field", "raw_field", "raw_krige"], fld_cnt=3 - ) - krige_name, krige_save = self.get_store_config( - krige_store, default=["field", "krige_var"], fld_cnt=2 + name, save = self.get_store_config(store=store, fld_cnt=3) + krige_name, krige_save = self.krige.get_store_config( + store=krige_store, fld_cnt=2 ) kwargs["mesh_type"] = mesh_type kwargs["only_mean"] = False # overwrite if given diff --git a/gstools/krige/base.py b/gstools/krige/base.py index b823a5f81..339c69eaa 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -125,6 +125,9 @@ class Krige(Field): Springer, Berlin, Heidelberg (2003) """ + default_names = ["field", "krige_var", "mean_field"] + """:class:`list`: Default field names.""" + def __init__( self, model, @@ -223,7 +226,7 @@ def __call__( """ return_var &= not only_mean # don't return variance when calc. mean fld_cnt = 2 if return_var else 1 - default = ["mean_field"] if only_mean else ["field", "krige_var"] + default = self.default_names[2] if only_mean else None name, save = self.get_store_config(store, default, fld_cnt) iso_pos, shape = self.pre_pos(pos, mesh_type) From 3eb0efcaf7aa87290cfeaabb5d656abf33eca46c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 17:19:55 +0200 Subject: [PATCH 31/49] Field: better name for default attribute --- gstools/field/base.py | 9 ++++----- gstools/field/cond_srf.py | 5 ++--- gstools/krige/base.py | 4 ++-- 3 files changed, 8 insertions(+), 10 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 11a3350df..cf8375aa1 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -77,7 +77,7 @@ class Field: Dimension of the field if no model is given. """ - default_names = ["field"] + default_field_names = ["field"] """:class:`list`: Default field names.""" def __init__( @@ -504,8 +504,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): self.delete_fields() info_ret["deleted"] = True del old_pos - if info: - return info_ret + return info_ret if info else None def get_store_config(self, store, default=None, fld_cnt=None): """ @@ -531,9 +530,9 @@ def get_store_config(self, store, default=None, fld_cnt=None): """ if default is None: if fld_cnt is None: - default = self.default_names[0] + default = self.default_field_names[0] else: - default = self.default_names + default = self.default_field_names # single field if fld_cnt is None: save = isinstance(store, str) or bool(store) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 4e24947b6..4c97e2c77 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -43,7 +43,7 @@ class CondSRF(Field): Have a look at the provided generators for further information. """ - default_names = ["field", "raw_field", "raw_krige"] + default_field_names = ["field", "raw_field", "raw_krige"] """:class:`list`: Default field names.""" def __init__(self, krige, generator="RandMeth", **generator_kwargs): @@ -215,8 +215,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): info_ret = super().set_pos(pos, mesh_type, info=True) if info_ret["deleted"]: self.krige.delete_fields() - if info: - return info_ret + return info_ret if info else None @property def pos(self): diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 339c69eaa..2d74972f9 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -125,7 +125,7 @@ class Krige(Field): Springer, Berlin, Heidelberg (2003) """ - default_names = ["field", "krige_var", "mean_field"] + default_field_names = ["field", "krige_var", "mean_field"] """:class:`list`: Default field names.""" def __init__( @@ -226,7 +226,7 @@ def __call__( """ return_var &= not only_mean # don't return variance when calc. mean fld_cnt = 2 if return_var else 1 - default = self.default_names[2] if only_mean else None + default = self.default_field_names[2] if only_mean else None name, save = self.get_store_config(store, default, fld_cnt) iso_pos, shape = self.pre_pos(pos, mesh_type) From e74221b36a72bdafbd6956abab835df0cb021f98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 17:50:08 +0200 Subject: [PATCH 32/49] Test: check reuse of field in CondSRF --- tests/test_condition.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/test_condition.py b/tests/test_condition.py index 587469d13..edb5b9669 100644 --- a/tests/test_condition.py +++ b/tests/test_condition.py @@ -1,8 +1,6 @@ # -*- coding: utf-8 -*- -""" -This is the unittest of CovModel class. -""" - +"""This is the unittest of CondSRF class.""" +from copy import copy import numpy as np import unittest import gstools as gs @@ -62,6 +60,10 @@ def test_simple(self): crf = gs.CondSRF(krige, seed=19970221) field_1 = crf.unstructured(self.pos[:dim]) field_2 = crf.structured(self.pos[:dim]) + # check reuse + raw_kr2 = copy(crf["raw_krige"]) + crf(seed=19970222) + self.assertTrue(np.allclose(raw_kr2, crf["raw_krige"])) for i, val in enumerate(self.cond_val): self.assertAlmostEqual(val, field_1[i], places=2) self.assertAlmostEqual(val, field_2[dim * (i,)], places=2) From 9d85f209e80a94e555c79fe2ba2ece0f428f8d6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 18:12:31 +0200 Subject: [PATCH 33/49] Field: remove mesh_type checks (never None) --- gstools/field/base.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index cf8375aa1..c4a94fcb4 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -278,8 +278,6 @@ def pre_pos(self, pos=None, mesh_type="unstructured", info=False): if pos is None: if self.pos is None: raise ValueError("Field: no position tuple 'pos' present") - if self.mesh_type is None: - raise ValueError("Field: no 'mesh_type' present") else: info_ret = self.set_pos(pos, mesh_type, info=True) if self.mesh_type != "unstructured": @@ -562,8 +560,6 @@ def pos(self): @pos.setter def pos(self, pos): - if self.mesh_type is None: - raise ValueError("Field.pos: can't set 'pos' without 'mesh_type'") if self.mesh_type == "unstructured": self._pos = np.asarray(pos, dtype=np.double).reshape(self.dim, -1) self._field_shape = np.shape(self._pos[0]) From 002ad6d8e0f2c75afaab3e37a48e6c00f9f9d08c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 18:12:48 +0200 Subject: [PATCH 34/49] test: more Field tests --- tests/test_field.py | 61 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) diff --git a/tests/test_field.py b/tests/test_field.py index 53e7c3367..9628e8f38 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -46,6 +46,67 @@ def test_raise(self): with self.assertRaises(ValueError): gs.field.Field(dim=3, mean=[1, 2]) + def test_pos_compare(self): + fld = gs.field.Field(dim=1) + fld.set_pos([1, 2]) + fld._dim = 2 + info = fld.set_pos([[1], [2]], info=True) + self.assertTrue(info["deleted"]) + info = fld.set_pos([[2], [3]], info=True) + self.assertTrue(info["deleted"]) + + def test_magic(self): + fld = gs.field.Field(dim=1) + f1 = np.array([0, 0], dtype=np.double) + f2 = np.array([2, 3], dtype=np.double) + fld([1, 2], store="f1") # default field with zeros + fld([1, 2], f2, store="f2") + fields1 = fld[:] + fields2 = fld[[0, 1]] + fields3 = fld[["f1", "f2"]] + fields4 = fld.all_fields + self.assertTrue(np.allclose([f1, f2], fields1)) + self.assertTrue(np.allclose([f1, f2], fields2)) + self.assertTrue(np.allclose([f1, f2], fields3)) + self.assertTrue(np.allclose([f1, f2], fields4)) + self.assertEqual(len(fld), 2) + self.assertTrue("f1" in fld) + self.assertTrue("f2" in fld) + self.assertFalse("f3" in fld) + # subscription + with self.assertRaises(KeyError): + fld["f3"] + with self.assertRaises(KeyError): + del fld["f3"] + with self.assertRaises(KeyError): + del fld[["f3"]] + del fld["f1"] + self.assertFalse("f1" in fld) + fld([1, 2], f1, store="f1") + del fld[-1] + self.assertFalse("f1" in fld) + fld([1, 2], f1, store="f1") + del fld[:] + self.assertEqual(len(fld), 0) + fld([1, 2], f1, store="f1") + del fld.field_names + self.assertEqual(len(fld), 0) + + def test_reuse(self): + fld = gs.field.Field(dim=1) + with self.assertRaises(ValueError): + fld() + with self.assertRaises(ValueError): + fld.post_field([1, 2]) + with self.assertRaises(ValueError): + fld.post_field([1, 2], process=False, name=0) + fld.set_pos([1, 2]) + with self.assertRaises(ValueError): + fld.structured() + fld.set_pos([1, 2], "structured") + with self.assertRaises(ValueError): + fld.unstructured() + if __name__ == "__main__": unittest.main() From 56560e3f063b8e930725f0eda354cf9a58252df1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 19:21:29 +0200 Subject: [PATCH 35/49] transform: add keep_mean to all methods with True by default --- gstools/transform/field.py | 138 ++++++++++++++++++++++++++++++------- 1 file changed, 113 insertions(+), 25 deletions(-) diff --git a/gstools/transform/field.py b/gstools/transform/field.py index 90ee6c6ec..0535adff4 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -178,7 +178,7 @@ def apply_function( field="field", store=True, process=False, - keep_mean=False, + keep_mean=True, **kwargs, ): """ @@ -200,7 +200,7 @@ def apply_function( of given Field instance. The default is False. keep_mean : :class:`bool`, optional Whether to keep the mean of the field if process=True. - The default is False. + The default is True. **kwargs Keyword arguments forwarded to given function. @@ -234,6 +234,7 @@ def binary( field="field", store=True, process=False, + keep_mean=True, ): """ Binary transformation. @@ -261,6 +262,9 @@ def binary( process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -270,7 +274,7 @@ def binary( if not process and divide is None: _check_for_default_normal(fld) if divide is None: - mean = 0.0 if process else fld.mean + mean = 0.0 if process and not keep_mean else fld.mean divide = mean if divide is None else divide upper = mean + np.sqrt(fld.model.sill) if upper is None else upper lower = mean - np.sqrt(fld.model.sill) if lower is None else lower @@ -284,7 +288,7 @@ def binary( field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) @@ -296,6 +300,7 @@ def discrete( field="field", store=True, process=False, + keep_mean=True, ): """ Discrete transformation. @@ -324,6 +329,9 @@ def discrete( process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -335,7 +343,7 @@ def discrete( kw = dict( values=values, thresholds=thresholds, - mean=0.0 if process else fld.mean, + mean=0.0 if process and not keep_mean else fld.mean, var=fld.model.sill, ) return apply_function( @@ -344,12 +352,20 @@ def discrete( field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) -def boxcox(fld, lmbda=1, shift=0, field="field", store=True, process=False): +def boxcox( + fld, + lmbda=1, + shift=0, + field="field", + store=True, + process=False, + keep_mean=True, +): """ (Inverse) Box-Cox transformation to denormalize data. @@ -378,6 +394,9 @@ def boxcox(fld, lmbda=1, shift=0, field="field", store=True, process=False): process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -391,12 +410,19 @@ def boxcox(fld, lmbda=1, shift=0, field="field", store=True, process=False): field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) -def zinnharvey(fld, conn="high", field="field", store=True, process=False): +def zinnharvey( + fld, + conn="high", + field="field", + store=True, + process=False, + keep_mean=True, +): """ Zinn and Harvey transformation to connect low or high values. @@ -417,6 +443,9 @@ def zinnharvey(fld, conn="high", field="field", store=True, process=False): process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -425,19 +454,29 @@ def zinnharvey(fld, conn="high", field="field", store=True, process=False): """ if not process: _check_for_default_normal(fld) - kw = dict(conn=conn, mean=0.0 if process else fld.mean, var=fld.model.sill) + kw = dict( + conn=conn, + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + ) return apply_function( fld=fld, function=array_zinnharvey, field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) -def normal_force_moments(fld, field="field", store=True, process=False): +def normal_force_moments( + fld, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Force moments of a normal distributed field. @@ -455,6 +494,9 @@ def normal_force_moments(fld, field="field", store=True, process=False): process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -463,19 +505,23 @@ def normal_force_moments(fld, field="field", store=True, process=False): """ if not process: _check_for_default_normal(fld) - kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, var=fld.model.sill + ) return apply_function( fld=fld, function=array_force_moments, field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) -def normal_to_lognormal(fld, field="field", store=True, process=False): +def normal_to_lognormal( + fld, field="field", store=True, process=False, keep_mean=True +): """ Transform normal distribution to log-normal distribution. @@ -493,6 +539,9 @@ def normal_to_lognormal(fld, field="field", store=True, process=False): process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -505,11 +554,17 @@ def normal_to_lognormal(fld, field="field", store=True, process=False): field=field, store=store, process=process, - keep_mean=True, # apply to normal field including mean + keep_mean=keep_mean, ) -def normal_to_uniform(fld, field="field", store=True, process=False): +def normal_to_uniform( + fld, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Transform normal distribution to uniform distribution on [0, 1]. @@ -519,23 +574,34 @@ def normal_to_uniform(fld, field="field", store=True, process=False): ---------- fld : :any:`Field` Field class containing a generated field. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. """ if not process: _check_for_default_normal(fld) - kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, var=fld.model.sill + ) return apply_function( fld=fld, function=array_to_uniform, field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) def normal_to_arcsin( - fld, a=None, b=None, field="field", store=True, process=False + fld, + a=None, + b=None, + field="field", + store=True, + process=False, + keep_mean=True, ): """ Transform normal distribution to the bimodal arcsin distribution. @@ -562,6 +628,9 @@ def normal_to_arcsin( process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -570,20 +639,31 @@ def normal_to_arcsin( """ if not process: _check_for_default_normal(fld) - kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill, a=a, b=b) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + a=a, + b=b, + ) return apply_function( fld=fld, function=array_to_arcsin, field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) def normal_to_uquad( - fld, a=None, b=None, field="field", store=True, process=False + fld, + a=None, + b=None, + field="field", + store=True, + process=False, + keep_mean=True, ): """ Transform normal distribution to U-quadratic distribution. @@ -610,6 +690,9 @@ def normal_to_uquad( process : :class:`bool`, optional Whether to process in/out fields with trend, normalizer and mean of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- @@ -618,14 +701,19 @@ def normal_to_uquad( """ if not process: _check_for_default_normal(fld) - kw = dict(mean=0.0 if process else fld.mean, var=fld.model.sill, a=a, b=b) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + a=a, + b=b, + ) return apply_function( fld=fld, function=array_to_uquad, field=field, store=store, process=process, - keep_mean=False, + keep_mean=keep_mean, **kw, ) From c049c1e15954c97d4eeb27a361ccbde6c76d9166 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 19:21:41 +0200 Subject: [PATCH 36/49] tests: add transform tests --- tests/test_srf.py | 40 --------- tests/test_transform.py | 184 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 184 insertions(+), 40 deletions(-) create mode 100644 tests/test_transform.py diff --git a/tests/test_srf.py b/tests/test_srf.py index 7cca60bf4..0a83067de 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -256,46 +256,6 @@ def test_mesh_pyvista(self): _ = srf.mesh(pv_mesh, seed=self.seed, points="points") self.assertTrue(np.allclose(field, pv_mesh["field"])) - def test_transform(self): - self.cov_model.dim = 2 - srf = gs.SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_force_moments(srf) # force ergodicity of the given field - self.assertAlmostEqual(srf.field.mean(), srf.mean) - self.assertAlmostEqual(srf.field.var(), srf.model.var) - tf.zinnharvey(srf) # make high values mostly connected - tf.normal_force_moments(srf) # force ergodicity of the given field - tf.normal_to_lognormal(srf) # log-normal - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_arcsin(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_uquad(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_uniform(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.binary(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.boxcox(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = np.linspace(np.min(srf.field), np.max(srf.field), 3) - tf.discrete(srf, values) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 1] - thresholds = [-0.9, 0.1] - tf.discrete(srf, values, thresholds) - np.testing.assert_array_equal(np.unique(srf.field), [-1, 0, 1]) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 1] - tf.discrete(srf, values, thresholds="arithmetic") - np.testing.assert_array_equal(np.unique(srf.field), [-1.0, 0.0, 1.0]) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 0.5, 1] - tf.discrete(srf, values, thresholds="equal") - np.testing.assert_array_equal(np.unique(srf.field), values) - def test_incomprrandmeth(self): self.cov_model = gs.Gaussian(dim=2, var=0.5, len_scale=1.0) srf = gs.SRF( diff --git a/tests/test_transform.py b/tests/test_transform.py new file mode 100644 index 000000000..089e94aeb --- /dev/null +++ b/tests/test_transform.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +This is the unittest of SRF class. +""" + +import unittest +import numpy as np +import gstools as gs +from gstools import transform as tf +import meshio + +HAS_PYVISTA = False +try: + import pyvista as pv + + HAS_PYVISTA = True +except ImportError: + pass + + +class TestSRF(unittest.TestCase): + def setUp(self): + self.cov_model = gs.Gaussian(dim=2, var=1.5, len_scale=4.0) + self.mean = 0.3 + self.mode_no = 100 + + self.seed = 825718662 + self.x_grid = np.linspace(0.0, 12.0, 48) + self.y_grid = np.linspace(0.0, 10.0, 46) + + self.methods = [ + "binary", + "boxcox", + "zinnharvey", + "normal_force_moments", + "normal_to_lognormal", + "normal_to_uniform", + "normal_to_arcsin", + "normal_to_uquad", + ] + + def test_transform_normal(self): + srf = gs.SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no) + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + for method in self.methods: + srf.transform(method, store=method) + std = np.sqrt(srf.model.var) + self.assertTrue(set(self.methods) == set(srf.field_names[1:])) + # force moments + self.assertAlmostEqual(srf["normal_force_moments"].mean(), srf.mean) + self.assertAlmostEqual(srf["normal_force_moments"].var(), std ** 2) + # binary + np.testing.assert_allclose( + np.unique(srf.binary), srf.mean + np.array([-std, std]) + ) + # boxcox + np.testing.assert_allclose( + srf.field, gs.normalizer.BoxCox().normalize(srf.boxcox) + ) + # lognormal + np.testing.assert_allclose(srf.field, np.log(srf.normal_to_lognormal)) + # unifrom + self.assertTrue(np.all(srf.normal_to_uniform < 1)) + self.assertTrue(np.all(srf.normal_to_uniform > 0)) + # how to test arcsin and uquad?! + + # discrete + values = [-1, 0, 1] + thresholds = [-0.9, 0.1] + srf.transform( + "discrete", values=values, thresholds=thresholds, store="f1" + ) + np.testing.assert_array_equal(np.unique(srf.f1), [-1, 0, 1]) + + values = [-1, 0, 1] + srf.transform( + "discrete", values=values, thresholds="arithmetic", store="f2" + ) + np.testing.assert_array_equal(np.unique(srf.f2), [-1.0, 0.0, 1.0]) + + values = [-1, 0, 0.5, 1] + srf.transform( + "discrete", values=values, thresholds="equal", store="f3" + ) + np.testing.assert_array_equal(np.unique(srf.f3), values) + + def test_transform_denormal(self): + srf_fail = gs.SRF( + model=self.cov_model, + mean=self.mean, + mode_no=self.mode_no, + trend=lambda x, y: x, + ) + srf_fail((self.x_grid, self.y_grid), mesh_type="structured") + with self.assertRaises(ValueError): + srf_fail.transform("zinnharvey") + + srf_fail = gs.SRF( + model=self.cov_model, + mean=lambda x, y: x, + mode_no=self.mode_no, + ) + srf_fail((self.x_grid, self.y_grid), mesh_type="structured") + with self.assertRaises(ValueError): + srf_fail.transform("zinnharvey") + + srf = gs.SRF( + model=self.cov_model, + mean=self.mean, + mode_no=self.mode_no, + normalizer=gs.normalizer.LogNormal, + ) + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + + for method in self.methods: + if method in ("normal_to_lognormal", "boxcox"): + continue + with self.assertRaises(ValueError): + srf.transform(method, store=method) + + for method in self.methods: + srf.transform(method, store=method, process=True) + std = np.sqrt(srf.model.var) + self.assertTrue(set(self.methods) == set(srf.field_names[1:])) + # force moments + self.assertAlmostEqual( + np.log(srf["normal_force_moments"]).mean(), srf.mean + ) + self.assertAlmostEqual( + np.log(srf["normal_force_moments"]).var(), std ** 2 + ) + # binary + np.testing.assert_allclose( + np.unique(np.log(srf.binary)), srf.mean + np.array([-std, std]) + ) + # boxcox + np.testing.assert_allclose( + np.log(srf.field), + gs.normalizer.BoxCox().normalize(np.log(srf.boxcox)), + ) + # lognormal + np.testing.assert_allclose(srf.field, np.log(srf.normal_to_lognormal)) + # unifrom + self.assertTrue(np.all(np.log(srf.normal_to_uniform) < 1)) + self.assertTrue(np.all(np.log(srf.normal_to_uniform) > 0)) + + # discrete + values = [-1, 0, 1] + thresholds = [-0.9, 0.1] + srf.transform( + "discrete", + values=values, + thresholds=thresholds, + store="f1", + process=True, + ) + np.testing.assert_array_equal(np.unique(np.log(srf.f1)), [-1, 0, 1]) + + values = [-1, 0, 1] + srf.transform( + "discrete", + values=values, + thresholds="arithmetic", + store="f2", + process=True, + ) + np.testing.assert_array_equal( + np.unique(np.log(srf.f2)), [-1.0, 0.0, 1.0] + ) + + values = [-1, 0, 0.5, 1] + srf.transform( + "discrete", + values=values, + thresholds="equal", + store="f3", + process=True, + ) + np.testing.assert_array_equal(np.unique(np.log(srf.f3)), values) + + +if __name__ == "__main__": + unittest.main() From 3923da04e8806e7b0aded58b7c67cb96f24cfa4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 3 Aug 2021 19:28:38 +0200 Subject: [PATCH 37/49] transform: own module for array transformations --- gstools/transform/__init__.py | 2 + gstools/transform/array.py | 351 ++++++++++++++++++++++++++++++++++ gstools/transform/field.py | 351 +--------------------------------- tests/test_transform.py | 16 +- 4 files changed, 364 insertions(+), 356 deletions(-) create mode 100644 gstools/transform/array.py diff --git a/gstools/transform/__init__.py b/gstools/transform/__init__.py index 3a2a50061..0000a5c73 100644 --- a/gstools/transform/__init__.py +++ b/gstools/transform/__init__.py @@ -59,6 +59,8 @@ normal_to_uniform, normal_to_arcsin, normal_to_uquad, +) +from gstools.transform.array import ( array_discrete, array_boxcox, array_zinnharvey, diff --git a/gstools/transform/array.py b/gstools/transform/array.py new file mode 100644 index 000000000..683cac7ed --- /dev/null +++ b/gstools/transform/array.py @@ -0,0 +1,351 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing array transformations. + +.. currentmodule:: gstools.transform.array + +The following functions are provided + +Transformations +^^^^^^^^^^^^^^^ + +.. autosummary:: + array_discrete + array_boxcox + array_zinnharvey + array_force_moments + array_to_lognormal + array_to_uniform + array_to_arcsin + array_to_uquad +""" +# pylint: disable=C0103, C0123, R0911 +from warnings import warn +import numpy as np +from scipy.special import erf, erfinv + +__all__ = [ + "array_discrete", + "array_boxcox", + "array_zinnharvey", + "array_force_moments", + "array_to_lognormal", + "array_to_uniform", + "array_to_arcsin", + "array_to_uquad", +] + + +def array_discrete( + field, values, thresholds="arithmetic", mean=None, var=None +): + """ + Discrete transformation. + + After this transformation, the field has only `len(values)` discrete + values. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + values : :any:`numpy.ndarray` + The discrete values the field will take + thresholds : :class:`str` or :any:`numpy.ndarray`, optional + the thresholds, where the value classes are separated + possible values are: + * "arithmetic": the mean of the 2 neighbouring values + * "equal": devide the field into equal parts + * an array of explicitly given thresholds + Default: "arithmetic" + mean : :class:`float`or :any:`None` + Mean of the field for "equal" thresholds. Default: np.mean(field) + var : :class:`float`or :any:`None` + Variance of the field for "equal" thresholds. Default: np.var(field) + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + if thresholds == "arithmetic": + # just in case, sort the values + values = np.sort(values) + thresholds = (values[1:] + values[:-1]) / 2 + elif thresholds == "equal": + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + values = np.asarray(values) + n = len(values) + p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] + rescale = np.sqrt(var * 2) + # use quantile of the normal distribution to get equal ratios + thresholds = mean + rescale * erfinv(2 * p - 1) + else: + if len(values) != len(thresholds) + 1: + raise ValueError( + "discrete transformation: len(values) != len(thresholds) + 1" + ) + values = np.asarray(values) + thresholds = np.asarray(thresholds) + # check thresholds + if not np.all(thresholds[:-1] < thresholds[1:]): + raise ValueError( + "discrete transformation: thresholds need to be ascending" + ) + # use a separate result so the intermediate results are not affected + result = np.empty_like(field) + # handle edge cases + result[field <= thresholds[0]] = values[0] + result[field > thresholds[-1]] = values[-1] + for i, value in enumerate(values[1:-1]): + result[ + np.logical_and(thresholds[i] < field, field <= thresholds[i + 1]) + ] = value + return result + + +def array_boxcox(field, lmbda=1, shift=0): + """ + (Inverse) Box-Cox transformation to denormalize data. + + After this transformation, the again Box-Cox transformed field is normal + distributed. + + See: https://en.wikipedia.org/wiki/Power_transform#Box%E2%80%93Cox_transformation + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + lmbda : :class:`float`, optional + The lambda parameter of the Box-Cox transformation. + For ``lmbda=0`` one obtains the log-normal transformation. + Default: ``1`` + shift : :class:`float`, optional + The shift parameter from the two-parametric Box-Cox transformation. + The field will be shifted by that value before transformation. + Default: ``0`` + """ + field = np.asarray(field) + result = field + shift + if np.isclose(lmbda, 0): + array_to_lognormal(result) + if np.min(result) < -1 / lmbda: + warn("Box-Cox: Some values will be cut off!") + return (np.maximum(lmbda * result + 1, 0)) ** (1 / lmbda) + + +def array_zinnharvey(field, conn="high", mean=None, var=None): + """ + Zinn and Harvey transformation to connect low or high values. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + conn : :class:`str`, optional + Desired connectivity. Either "low" or "high". + Default: "high" + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + result = np.abs((field - mean) / np.sqrt(var)) + result = np.sqrt(2) * erfinv(2 * erf(result / np.sqrt(2)) - 1) + if conn == "high": + result = -result + return result * np.sqrt(var) + mean + + +def array_force_moments(field, mean=0, var=1): + """ + Force moments of a normal distributed field. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float`, optional + Desired mean of the field. + Default: 0 + var : :class:`float` or :any:`None`, optional + Desired variance of the field. + Default: 1 + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + var_in = np.var(field) + mean_in = np.mean(field) + rescale = np.sqrt(var / var_in) + return rescale * (field - mean_in) + mean + + +def array_to_lognormal(field): + """ + Transform normal distribution to log-normal distribution. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + return np.exp(field) + + +def array_to_uniform(field, mean=None, var=None): + """ + Transform normal distribution to uniform distribution on [0, 1]. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + return 0.5 * (1 + erf((field - mean) / np.sqrt(2 * var))) + + +def array_to_arcsin(field, mean=None, var=None, a=None, b=None): + """ + Transform normal distribution to arcsin distribution. + + See: https://en.wikipedia.org/wiki/Arcsine_distribution + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the mean will be calculated. + Default: :any:`None` + a : :class:`float`, optional + Parameter a of the arcsin distribution (lower bound). + Default: keep mean and variance + b : :class:`float`, optional + Parameter b of the arcsin distribution (upper bound). + Default: keep mean and variance + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(2.0 * var) if a is None else float(a) + b = mean + np.sqrt(2.0 * var) if b is None else float(b) + return _uniform_to_arcsin(array_to_uniform(field, mean, var), a, b) + + +def array_to_uquad(field, mean=None, var=None, a=None, b=None): + """ + Transform normal distribution to U-quadratic distribution. + + See: https://en.wikipedia.org/wiki/U-quadratic_distribution + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + a : :class:`float`, optional + Parameter a of the U-quadratic distribution (lower bound). + Default: keep mean and variance + b : :class:`float`, optional + Parameter b of the U-quadratic distribution (upper bound). + Default: keep mean and variance + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(5.0 / 3.0 * var) if a is None else float(a) + b = mean + np.sqrt(5.0 / 3.0 * var) if b is None else float(b) + return _uniform_to_uquad(array_to_uniform(field, mean, var), a, b) + + +def _uniform_to_arcsin(field, a=0, b=1): + """ + PPF of your desired distribution. + + The PPF is the inverse of the CDF and is used to sample a distribution + from uniform distributed values on [0, 1] + + in this case: the arcsin distribution + See: https://en.wikipedia.org/wiki/Arcsine_distribution + """ + field = np.asarray(field) + return (b - a) * np.sin(np.pi * 0.5 * field) ** 2 + a + + +def _uniform_to_uquad(field, a=0, b=1): + """ + PPF of your desired distribution. + + The PPF is the inverse of the CDF and is used to sample a distribution + from uniform distributed values on [0, 1] + + in this case: the U-quadratic distribution + See: https://en.wikipedia.org/wiki/U-quadratic_distribution + """ + field = np.asarray(field) + al = 12 / (b - a) ** 3 + be = (a + b) / 2 + ga = (a - b) ** 3 / 8 + y_raw = 3 * field / al + ga + result = np.zeros_like(y_raw) + result[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) + result[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) + return result + be diff --git a/gstools/transform/field.py b/gstools/transform/field.py index 0535adff4..d855630e6 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -26,29 +26,24 @@ normal_to_uniform normal_to_arcsin normal_to_uquad - -Low-Level Routines -^^^^^^^^^^^^^^^^^^ - -.. autosummary:: - array_discrete - array_boxcox - array_zinnharvey - array_force_moments - array_to_lognormal - array_to_uniform - array_to_arcsin - array_to_uquad """ # pylint: disable=C0103, C0123, R0911 -from warnings import warn import numpy as np -from scipy.special import erf, erfinv from gstools.normalizer import ( Normalizer, remove_trend_norm_mean, apply_mean_norm_trend, ) +from gstools.transform.array import ( + array_discrete, + array_boxcox, + array_zinnharvey, + array_force_moments, + array_to_lognormal, + array_to_uniform, + array_to_arcsin, + array_to_uquad, +) __all__ = [ "apply", @@ -62,14 +57,6 @@ "normal_to_uniform", "normal_to_arcsin", "normal_to_uquad", - "array_discrete", - "array_boxcox", - "array_zinnharvey", - "array_force_moments", - "array_to_lognormal", - "array_to_uniform", - "array_to_arcsin", - "array_to_uquad", ] @@ -716,321 +703,3 @@ def normal_to_uquad( keep_mean=keep_mean, **kw, ) - - -# low level functions - - -def array_discrete( - field, values, thresholds="arithmetic", mean=None, var=None -): - """ - Discrete transformation. - - After this transformation, the field has only `len(values)` discrete - values. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - values : :any:`numpy.ndarray` - The discrete values the field will take - thresholds : :class:`str` or :any:`numpy.ndarray`, optional - the thresholds, where the value classes are separated - possible values are: - * "arithmetic": the mean of the 2 neighbouring values - * "equal": devide the field into equal parts - * an array of explicitly given thresholds - Default: "arithmetic" - mean : :class:`float`or :any:`None` - Mean of the field for "equal" thresholds. Default: np.mean(field) - var : :class:`float`or :any:`None` - Variance of the field for "equal" thresholds. Default: np.var(field) - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - if thresholds == "arithmetic": - # just in case, sort the values - values = np.sort(values) - thresholds = (values[1:] + values[:-1]) / 2 - elif thresholds == "equal": - mean = np.mean(field) if mean is None else float(mean) - var = np.var(field) if var is None else float(var) - values = np.asarray(values) - n = len(values) - p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] - rescale = np.sqrt(var * 2) - # use quantile of the normal distribution to get equal ratios - thresholds = mean + rescale * erfinv(2 * p - 1) - else: - if len(values) != len(thresholds) + 1: - raise ValueError( - "discrete transformation: len(values) != len(thresholds) + 1" - ) - values = np.asarray(values) - thresholds = np.asarray(thresholds) - # check thresholds - if not np.all(thresholds[:-1] < thresholds[1:]): - raise ValueError( - "discrete transformation: thresholds need to be ascending" - ) - # use a separate result so the intermediate results are not affected - result = np.empty_like(field) - # handle edge cases - result[field <= thresholds[0]] = values[0] - result[field > thresholds[-1]] = values[-1] - for i, value in enumerate(values[1:-1]): - result[ - np.logical_and(thresholds[i] < field, field <= thresholds[i + 1]) - ] = value - return result - - -def array_boxcox(field, lmbda=1, shift=0): - """ - (Inverse) Box-Cox transformation to denormalize data. - - After this transformation, the again Box-Cox transformed field is normal - distributed. - - See: https://en.wikipedia.org/wiki/Power_transform#Box%E2%80%93Cox_transformation - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - lmbda : :class:`float`, optional - The lambda parameter of the Box-Cox transformation. - For ``lmbda=0`` one obtains the log-normal transformation. - Default: ``1`` - shift : :class:`float`, optional - The shift parameter from the two-parametric Box-Cox transformation. - The field will be shifted by that value before transformation. - Default: ``0`` - """ - field = np.asarray(field) - result = field + shift - if np.isclose(lmbda, 0): - array_to_lognormal(result) - if np.min(result) < -1 / lmbda: - warn("Box-Cox: Some values will be cut off!") - return (np.maximum(lmbda * result + 1, 0)) ** (1 / lmbda) - - -def array_zinnharvey(field, conn="high", mean=None, var=None): - """ - Zinn and Harvey transformation to connect low or high values. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - conn : :class:`str`, optional - Desired connectivity. Either "low" or "high". - Default: "high" - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the variance will be calculated. - Default: :any:`None` - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - mean = np.mean(field) if mean is None else float(mean) - var = np.var(field) if var is None else float(var) - result = np.abs((field - mean) / np.sqrt(var)) - result = np.sqrt(2) * erfinv(2 * erf(result / np.sqrt(2)) - 1) - if conn == "high": - result = -result - return result * np.sqrt(var) + mean - - -def array_force_moments(field, mean=0, var=1): - """ - Force moments of a normal distributed field. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - mean : :class:`float`, optional - Desired mean of the field. - Default: 0 - var : :class:`float` or :any:`None`, optional - Desired variance of the field. - Default: 1 - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - var_in = np.var(field) - mean_in = np.mean(field) - rescale = np.sqrt(var / var_in) - return rescale * (field - mean_in) + mean - - -def array_to_lognormal(field): - """ - Transform normal distribution to log-normal distribution. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - return np.exp(field) - - -def array_to_uniform(field, mean=None, var=None): - """ - Transform normal distribution to uniform distribution on [0, 1]. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the variance will be calculated. - Default: :any:`None` - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - mean = np.mean(field) if mean is None else float(mean) - var = np.var(field) if var is None else float(var) - return 0.5 * (1 + erf((field - mean) / np.sqrt(2 * var))) - - -def array_to_arcsin(field, mean=None, var=None, a=None, b=None): - """ - Transform normal distribution to arcsin distribution. - - See: https://en.wikipedia.org/wiki/Arcsine_distribution - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the mean will be calculated. - Default: :any:`None` - a : :class:`float`, optional - Parameter a of the arcsin distribution (lower bound). - Default: keep mean and variance - b : :class:`float`, optional - Parameter b of the arcsin distribution (upper bound). - Default: keep mean and variance - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - mean = np.mean(field) if mean is None else float(mean) - var = np.var(field) if var is None else float(var) - a = mean - np.sqrt(2.0 * var) if a is None else float(a) - b = mean + np.sqrt(2.0 * var) if b is None else float(b) - return _uniform_to_arcsin(array_to_uniform(field, mean, var), a, b) - - -def array_to_uquad(field, mean=None, var=None, a=None, b=None): - """ - Transform normal distribution to U-quadratic distribution. - - See: https://en.wikipedia.org/wiki/U-quadratic_distribution - - Parameters - ---------- - field : :class:`numpy.ndarray` - Normal distributed values. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the variance will be calculated. - Default: :any:`None` - a : :class:`float`, optional - Parameter a of the U-quadratic distribution (lower bound). - Default: keep mean and variance - b : :class:`float`, optional - Parameter b of the U-quadratic distribution (upper bound). - Default: keep mean and variance - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - field = np.asarray(field) - mean = np.mean(field) if mean is None else float(mean) - var = np.var(field) if var is None else float(var) - a = mean - np.sqrt(5.0 / 3.0 * var) if a is None else float(a) - b = mean + np.sqrt(5.0 / 3.0 * var) if b is None else float(b) - return _uniform_to_uquad(array_to_uniform(field, mean, var), a, b) - - -def _uniform_to_arcsin(field, a=0, b=1): - """ - PPF of your desired distribution. - - The PPF is the inverse of the CDF and is used to sample a distribution - from uniform distributed values on [0, 1] - - in this case: the arcsin distribution - See: https://en.wikipedia.org/wiki/Arcsine_distribution - """ - field = np.asarray(field) - return (b - a) * np.sin(np.pi * 0.5 * field) ** 2 + a - - -def _uniform_to_uquad(field, a=0, b=1): - """ - PPF of your desired distribution. - - The PPF is the inverse of the CDF and is used to sample a distribution - from uniform distributed values on [0, 1] - - in this case: the U-quadratic distribution - See: https://en.wikipedia.org/wiki/U-quadratic_distribution - """ - field = np.asarray(field) - al = 12 / (b - a) ** 3 - be = (a + b) / 2 - ga = (a - b) ** 3 / 8 - y_raw = 3 * field / al + ga - result = np.zeros_like(y_raw) - result[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) - result[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) - return result + be diff --git a/tests/test_transform.py b/tests/test_transform.py index 089e94aeb..02b2fc47a 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -1,22 +1,8 @@ -#!/usr/bin/env python # -*- coding: utf-8 -*- -""" -This is the unittest of SRF class. -""" - +"""This is the unittest of the transform submodule.""" import unittest import numpy as np import gstools as gs -from gstools import transform as tf -import meshio - -HAS_PYVISTA = False -try: - import pyvista as pv - - HAS_PYVISTA = True -except ImportError: - pass class TestSRF(unittest.TestCase): From 4b44bbd95f93afe2fad822f358cd2c9702c9994c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 10:04:18 +0200 Subject: [PATCH 38/49] transform: bugfix in boxcox for lmbda=0 --- gstools/transform/array.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/transform/array.py b/gstools/transform/array.py index 683cac7ed..6ec9fd132 100644 --- a/gstools/transform/array.py +++ b/gstools/transform/array.py @@ -131,7 +131,7 @@ def array_boxcox(field, lmbda=1, shift=0): field = np.asarray(field) result = field + shift if np.isclose(lmbda, 0): - array_to_lognormal(result) + return array_to_lognormal(result) if np.min(result) < -1 / lmbda: warn("Box-Cox: Some values will be cut off!") return (np.maximum(lmbda * result + 1, 0)) ** (1 / lmbda) From a363919afa2294f17d6aacd6f0935981b072b605 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 10:04:48 +0200 Subject: [PATCH 39/49] Field: remove redundant pos check in post_field --- gstools/field/base.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index c4a94fcb4..d06ceb150 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -317,8 +317,6 @@ def post_field(self, field, name="field", process=True, save=True): raise ValueError("post_field: no 'field_shape' present.") field = np.asarray(field, dtype=np.double).reshape(self.field_shape) if process: - if self.pos is None: - raise ValueError("post_field: no 'pos' tuple set for field.") field = apply_mean_norm_trend( pos=self.pos, field=field, From f0f3f59577b0a45f81123b77240ea9a5b6a000bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 10:05:19 +0200 Subject: [PATCH 40/49] tests: add tests for some missing lines --- tests/test_field.py | 10 +++++++++- tests/test_transform.py | 19 +++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/tests/test_field.py b/tests/test_field.py index 9628e8f38..a5209dcb9 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -91,16 +91,24 @@ def test_magic(self): fld([1, 2], f1, store="f1") del fld.field_names self.assertEqual(len(fld), 0) + # store config (missing check) + name, save = fld.get_store_config(store="fld", fld_cnt=1) + self.assertEqual(name, ["fld"]) + self.assertTrue(save[0]) def test_reuse(self): fld = gs.field.Field(dim=1) + # no pos tuple with self.assertRaises(ValueError): fld() + # no field shape with self.assertRaises(ValueError): fld.post_field([1, 2]) + # bad name + fld.set_pos([1, 2]) with self.assertRaises(ValueError): fld.post_field([1, 2], process=False, name=0) - fld.set_pos([1, 2]) + # incompatible reuse with self.assertRaises(ValueError): fld.structured() fld.set_pos([1, 2], "structured") diff --git a/tests/test_transform.py b/tests/test_transform.py index 02b2fc47a..641f7d013 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -44,8 +44,12 @@ def test_transform_normal(self): np.testing.assert_allclose( srf.field, gs.normalizer.BoxCox().normalize(srf.boxcox) ) + with self.assertWarns(Warning): + srf.transform("boxcox", store="boxcox_warn", lmbda=2) # lognormal np.testing.assert_allclose(srf.field, np.log(srf.normal_to_lognormal)) + srf.transform("boxcox", store="boxcox2", lmbda=0) + np.testing.assert_allclose(srf.boxcox2, srf.normal_to_lognormal) # unifrom self.assertTrue(np.all(srf.normal_to_uniform < 1)) self.assertTrue(np.all(srf.normal_to_uniform > 0)) @@ -70,6 +74,21 @@ def test_transform_normal(self): "discrete", values=values, thresholds="equal", store="f3" ) np.testing.assert_array_equal(np.unique(srf.f3), values) + # checks + with self.assertRaises(ValueError): + srf.transform("discrete", values=values, thresholds=[1]) + with self.assertRaises(ValueError): + srf.transform("discrete", values=values, thresholds=[1, 3, 2]) + + # function + srf.transform("function", function=lambda x: 2 * x, store="f4") + np.testing.assert_array_equal(2 * srf.field, srf.f4) + with self.assertRaises(ValueError): + srf.transform("function", function=None) + + # unknown method + with self.assertRaises(ValueError): + srf.transform("foobar") def test_transform_denormal(self): srf_fail = gs.SRF( From ef049d8ee2c853fababbfeddf8653c50832abf51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 10:18:21 +0200 Subject: [PATCH 41/49] Field: some repr fixes in subclasses --- gstools/field/cond_srf.py | 4 ++-- gstools/krige/base.py | 5 ----- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 4c97e2c77..f869ff165 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -297,6 +297,6 @@ def value_type(self, value_type): def __repr__(self): """Return String representation.""" - return "CondSRF(krige={0}, generator={1})".format( - self.krige, self.generator.name + return "{0}(krige={1}, generator={2})".format( + self.name, self.krige, self.generator.name ) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 2d74972f9..8006fb72a 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -698,11 +698,6 @@ def ext_drift_no(self): """:class:`int`: Number of external drift values per point.""" return self.cond_ext_drift.shape[0] - @property - def name(self): - """:class:`str`: The name of the kriging class.""" - return self.__class__.__name__ - def __repr__(self): """Return String representation.""" return "{0}(model={1}, cond_no={2}{3})".format( From f3d36ca7e888ee46bdbf06fd7eba33c36aba909c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 10:21:47 +0200 Subject: [PATCH 42/49] transform: correct doc link --- gstools/field/base.py | 3 ++- gstools/transform/field.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index d06ceb150..9323c0adc 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -355,7 +355,8 @@ def transform( Parameters ---------- method : :class:`str` - Method to use. See :any:`transform` for available transformations. + Method to use. + See :any:`gstools.transform` for available transformations. field : :class:`str`, optional Name of field to be transformed. The default is "field". store : :class:`str` or :class:`bool`, optional diff --git a/gstools/transform/field.py b/gstools/transform/field.py index d855630e6..6b9076543 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -110,7 +110,8 @@ def apply(fld, method, field="field", store=True, process=False, **kwargs): fld : :any:`Field` Field class containing a generated field. method : :class:`str` - Method to use. See :any:`transform` for available transformations. + Method to use. + See :any:`gstools.transform` for available transformations. field : :class:`str`, optional Name of field to be transformed. The default is "field". store : :class:`str` or :class:`bool`, optional From 2808a1bcb3228c02316ef0f53ed45f0c07e646d8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 23:03:07 +0200 Subject: [PATCH 43/49] Special: bugfix for inc_gamma for integer s<0; added inc_gamma_low --- CHANGELOG.md | 4 ++++ gstools/tools/special.py | 25 +++++++++++++++++++++++-- 2 files changed, 27 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b7a2239b8..be0e69b14 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,10 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197) - only recalculate kriging field if `pos` tuple changed (optimized ensemble generation) - performance improvement by using `np.asarray` instead of `np.array` where possible - updated examples to use new features +- added incomplete lower gamma function `inc_gamma_low` (for TPLGaussian spectral density) + +### Bugfixes +- `inc_gamma` was defined wrong for integer `s < 0` ## [1.3.2] - Pure Pink - 2021-07 diff --git a/gstools/tools/special.py b/gstools/tools/special.py index b6990ec89..f71ae12da 100644 --- a/gstools/tools/special.py +++ b/gstools/tools/special.py @@ -8,6 +8,7 @@ .. autosummary:: inc_gamma + inc_gamma_low exp_int inc_beta tplstable_cor @@ -21,6 +22,7 @@ __all__ = [ "confidence_scaling", "inc_gamma", + "inc_gamma_low", "exp_int", "inc_beta", "tplstable_cor", @@ -64,12 +66,31 @@ def inc_gamma(s, x): if np.isclose(s, 0): return sps.exp1(x) if np.isclose(s, np.around(s)) and s < -0.5: - return x ** (s - 1) * sps.expn(int(1 - np.around(s)), x) + return x ** s * sps.expn(int(1 - np.around(s)), x) if s < 0: return (inc_gamma(s + 1, x) - x ** s * np.exp(-x)) / s return sps.gamma(s) * sps.gammaincc(s, x) +def inc_gamma_low(s, x): + r"""Calculate the lower incomplete gamma function. + + Given by: :math:`\gamma(s,x) = \int_0^x t^{s-1}\,e^{-t}\,{\rm d}t` + + Parameters + ---------- + s : :class:`float` + exponent in the integral + x : :class:`numpy.ndarray` + input values + """ + if np.isclose(s, np.around(s)) and s < 0.5: + return np.full_like(x, np.inf, dtype=np.double) + if s < 0: + return (inc_gamma_low(s + 1, x) + x ** s * np.exp(-x)) / s + return sps.gamma(s) * sps.gammainc(s, x) + + def exp_int(s, x): r"""Calculate the exponential integral :math:`E_s(x)`. @@ -225,7 +246,7 @@ def tpl_gau_spec_dens(k, dim, len_scale, hurst, len_low=0.0): z_nz = np.logical_not(z_gz) # near zero a = hurst + dim / 2.0 fac = (len_scale / 2.0) ** dim * hurst / np.pi ** (dim / 2.0) - res[z_gz] = fac * (sps.gamma(a) - inc_gamma(a, z[z_gz])) / z[z_gz] ** a + res[z_gz] = fac * inc_gamma_low(a, z[z_gz]) / z[z_gz] ** a # first order approximation for z near zero res[z_nz] = fac * (1.0 / a - z[z_nz] / (a + 1.0)) return res From bb57cd187a19a94ac743adbdef14b845494ae372 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 4 Aug 2021 23:24:05 +0200 Subject: [PATCH 44/49] Special: add inc_gamma_low to tools; update doc gen --- docs/source/random.rst | 2 -- docs/source/tools.rst | 2 -- gstools/random/__init__.py | 6 ++++++ gstools/tools/__init__.py | 9 +++++++++ 4 files changed, 15 insertions(+), 4 deletions(-) diff --git a/docs/source/random.rst b/docs/source/random.rst index 4bce61b0e..c43314996 100644 --- a/docs/source/random.rst +++ b/docs/source/random.rst @@ -2,8 +2,6 @@ gstools.random ============== .. automodule:: gstools.random - :members: - :undoc-members: .. raw:: latex diff --git a/docs/source/tools.rst b/docs/source/tools.rst index f6c86fe4a..e3e91643a 100644 --- a/docs/source/tools.rst +++ b/docs/source/tools.rst @@ -2,8 +2,6 @@ gstools.tools ============= .. automodule:: gstools.tools - :members: - :undoc-members: .. raw:: latex diff --git a/gstools/random/__init__.py b/gstools/random/__init__.py index ec78b43c0..6c02da31f 100644 --- a/gstools/random/__init__.py +++ b/gstools/random/__init__.py @@ -8,18 +8,24 @@ ^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + RNG Seed Generator ^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + MasterRNG Distribution factory ^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + dist_gen ---- diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index e03529c3b..15e3c2c3f 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -8,6 +8,8 @@ ^^^^^^ .. autosummary:: + :toctree: generated + vtk_export vtk_export_structured vtk_export_unstructured @@ -19,8 +21,11 @@ ^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + confidence_scaling inc_gamma + inc_gamma_low exp_int inc_beta tplstable_cor @@ -31,6 +36,8 @@ ^^^^^^^^^ .. autosummary:: + :toctree: generated + rotated_main_axes set_angles set_anis @@ -68,6 +75,7 @@ from gstools.tools.special import ( confidence_scaling, inc_gamma, + inc_gamma_low, exp_int, inc_beta, tplstable_cor, @@ -107,6 +115,7 @@ "to_vtk_unstructured", "confidence_scaling", "inc_gamma", + "inc_gamma_low", "exp_int", "inc_beta", "tplstable_cor", From 12e88c2eeb8836a8b7db9cbf3c8ddd5b615efc6a Mon Sep 17 00:00:00 2001 From: LSchueler Date: Fri, 6 Aug 2021 17:27:18 +0200 Subject: [PATCH 45/49] Fix type --- gstools/field/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 9323c0adc..2af1f6a1e 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -505,7 +505,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): def get_store_config(self, store, default=None, fld_cnt=None): """ - Get sotrage configuration from given selection. + Get storage configuration from given selection. Parameters ---------- From ce7cbd70c94ec1d0a0c4b161d5264990381b8abf Mon Sep 17 00:00:00 2001 From: LSchueler Date: Fri, 6 Aug 2021 17:27:49 +0200 Subject: [PATCH 46/49] Correct docs --- gstools/field/base.py | 2 +- gstools/field/cond_srf.py | 2 +- gstools/field/srf.py | 2 +- gstools/krige/base.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 2af1f6a1e..bbee34828 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -158,7 +158,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions field : :class:`numpy.ndarray` or :any:`None`, optional diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index f869ff165..ffba2c32c 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -73,7 +73,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions seed : :class:`int`, optional diff --git a/gstools/field/srf.py b/gstools/field/srf.py index d5e273f16..86d4ef6db 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -114,7 +114,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions seed : :class:`int`, optional diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 8006fb72a..0c87d0f1a 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -190,7 +190,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions (x, [y, z]) mesh_type : :class:`str`, optional From 691275b440c281b8aa49b6f65b1a6c4b722f55e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 7 Aug 2021 17:56:10 +0200 Subject: [PATCH 47/49] Examples: use a seed-generator for ensemble examples --- examples/06_conditioned_fields/00_condition_ensemble.py | 4 +++- examples/06_conditioned_fields/01_2D_condition_ensemble.py | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/examples/06_conditioned_fields/00_condition_ensemble.py b/examples/06_conditioned_fields/00_condition_ensemble.py index ad11c0c32..551e867b2 100644 --- a/examples/06_conditioned_fields/00_condition_ensemble.py +++ b/examples/06_conditioned_fields/00_condition_ensemble.py @@ -29,10 +29,12 @@ cond_srf.set_pos(gridx) ############################################################################### +# To generate the ensemble we will use a seed-generator. # We can specify individual names for each field by the keyword `store`: +seed = gs.random.MasterRNG(20170519) for i in range(100): - cond_srf(seed=i, store=f"f{i}") + cond_srf(seed=seed(), store=f"f{i}") label = "Conditioned ensemble" if i == 0 else None plt.plot(gridx, cond_srf[f"f{i}"], color="k", alpha=0.1, label=label) diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py index 97c906f2d..b77ee0f9a 100644 --- a/examples/06_conditioned_fields/01_2D_condition_ensemble.py +++ b/examples/06_conditioned_fields/01_2D_condition_ensemble.py @@ -23,15 +23,16 @@ cond_srf.set_pos([x, y], "structured") ############################################################################### -# We create a list containing the generated conditioned fields. +# To generate the ensemble we will use a seed-generator. # By specifying ``store=[f"fld{i}", False, False]``, only the conditioned field # is stored with the specified name. The raw random field and the raw kriging # field is not stored. This way, we can access each conditioned field by index # ``cond_srf[i]``: +seed = gs.random.MasterRNG(20170519) ens_no = 4 for i in range(ens_no): - cond_srf(seed=i, store=[f"fld{i}", False, False]) + cond_srf(seed=seed(), store=[f"fld{i}", False, False]) ############################################################################### # Now let's have a look at the pairwise differences between the generated From b3ffa2d98b9c49fba07c0b3359fba586ea10662d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 7 Aug 2021 18:04:37 +0200 Subject: [PATCH 48/49] Field: add warning section in set_pos and pre_pos --- gstools/field/base.py | 10 ++++++++++ gstools/field/cond_srf.py | 5 +++++ 2 files changed, 15 insertions(+) diff --git a/gstools/field/base.py b/gstools/field/base.py index bbee34828..de8659c42 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -273,6 +273,11 @@ def pre_pos(self, pos=None, mesh_type="unstructured", info=False): Shape of the resulting field. info : :class:`dict`, optional Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. """ info_ret = {"deleted": False} if pos is None: @@ -489,6 +494,11 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): ------- info : :class:`dict` Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. """ info_ret = {"deleted": False} old_type = copy(self.mesh_type) diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index ffba2c32c..fd8db8c07 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -211,6 +211,11 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): ------- info : :class:`dict` Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. """ info_ret = super().set_pos(pos, mesh_type, info=True) if info_ret["deleted"]: From c555ba78b904661341b9d1db7376a2b3bc0172e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 7 Aug 2021 18:17:41 +0200 Subject: [PATCH 49/49] Field.set_pos: add optional tag for return --- gstools/field/base.py | 2 +- gstools/field/cond_srf.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index de8659c42..10c22126b 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -492,7 +492,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): Returns ------- - info : :class:`dict` + info : :class:`dict`, optional Information about settings. Warnings diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index fd8db8c07..5d424c6ec 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -209,7 +209,7 @@ def set_pos(self, pos, mesh_type="unstructured", info=False): Returns ------- - info : :class:`dict` + info : :class:`dict`, optional Information about settings. Warnings