diff --git a/pyproject.toml b/pyproject.toml index 090c7c63..88185bf3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -155,7 +155,7 @@ target-version = [ max-locals = 50 max-branches = 30 max-statements = 85 - max-attributes = 25 + max-attributes = 30 max-public-methods = 80 [tool.cibuildwheel] diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 11e63a2b..7488dc9c 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -54,6 +54,7 @@ .. autosummary:: CovModel + SumModel Covariance Models ^^^^^^^^^^^^^^^^^ @@ -62,6 +63,7 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~ .. autosummary:: + Nugget Gaussian Exponential Matern @@ -153,9 +155,11 @@ JBessel, Linear, Matern, + Nugget, Rational, Spherical, Stable, + SumModel, SuperSpherical, TPLExponential, TPLGaussian, @@ -198,6 +202,8 @@ __all__ += ["transform", "normalizer", "config"] __all__ += [ "CovModel", + "SumModel", + "Nugget", "Gaussian", "Exponential", "Matern", diff --git a/src/gstools/covmodel/__init__.py b/src/gstools/covmodel/__init__.py index 28ab81f2..76920c61 100644 --- a/src/gstools/covmodel/__init__.py +++ b/src/gstools/covmodel/__init__.py @@ -19,6 +19,7 @@ :toctree: CovModel + SumModel Covariance Models ^^^^^^^^^^^^^^^^^ @@ -27,6 +28,7 @@ .. autosummary:: :toctree: + Nugget Gaussian Exponential Matern @@ -53,7 +55,7 @@ TPLSimple """ -from gstools.covmodel.base import CovModel +from gstools.covmodel.base import CovModel, SumModel from gstools.covmodel.models import ( Circular, Cubic, @@ -64,6 +66,7 @@ JBessel, Linear, Matern, + Nugget, Rational, Spherical, Stable, @@ -78,6 +81,8 @@ __all__ = [ "CovModel", + "SumModel", + "Nugget", "Gaussian", "Exponential", "Matern", diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 23e19881..38a47491 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -9,7 +9,7 @@ CovModel """ -# pylint: disable=C0103, R0201, E1101, C0302, W0613 +# pylint: disable=C0103, R0201, E1101, C0302, W0613, W0231 import copy import numpy as np @@ -18,10 +18,19 @@ from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram +from gstools.covmodel.sum_tools import ( + default_mod_kwargs, + sum_check, + sum_default_arg_bounds, + sum_default_opt_arg_bounds, + sum_model_repr, + sum_set_norm_len_ratios, + sum_set_norm_var_ratios, +) from gstools.covmodel.tools import ( _init_subclass, check_arg_bounds, - check_bounds, + check_arg_in_bounds, compare, default_arg_from_bounds, model_repr, @@ -127,12 +136,6 @@ class CovModel: If given, the model dimension will be determined from this spatial dimension and the possible temporal dimension if temporal is ture. Default: None - var_raw : :class:`float` or :any:`None`, optional - raw variance of the model which will be multiplied with - :any:`CovModel.var_factor` to result in the actual variance. - If given, ``var`` will be ignored. - (This is just for models that override :any:`CovModel.var_factor`) - Default: :any:`None` hankel_kw: :class:`dict` or :any:`None`, optional Modify the init-arguments of :any:`hankel.SymmetricFourierTransform` @@ -144,6 +147,8 @@ class CovModel: If present, they are described in the section `Other Parameters`. """ + _add_doc = True + def __init__( self, dim=3, @@ -159,7 +164,6 @@ def __init__( geo_scale=RADIAN_SCALE, temporal=False, spatial_dim=None, - var_raw=None, hankel_kw=None, **opt_arg, ): @@ -169,6 +173,8 @@ def __init__( if not hasattr(self, "variogram"): raise TypeError("Don't instantiate 'CovModel' directly!") + # indicator for initialization status (True after __init__) + self._init = False # prepare dim setting self._dim = None self._hankel_kw = None @@ -205,9 +211,12 @@ def __init__( # set parameters self.rescale = rescale + self._var = float(var) self._nugget = float(nugget) # set anisotropy and len_scale, disable anisotropy for latlon models + if integral_scale is not None: + len_scale = integral_scale self._len_scale, self._anis = set_len_anis( self.dim, len_scale, anis, self.latlon ) @@ -215,26 +224,17 @@ def __init__( self.dim, angles, self.latlon, self.temporal ) - # set var at last, because of the var_factor (to be right initialized) - if var_raw is None: - self._var = None - self.var = var - else: - self._var = float(var_raw) - self._integral_scale = None - self.integral_scale = integral_scale - # set var again, if int_scale affects var_factor - if var_raw is None: - self._var = None - self.var = var - else: - self._var = float(var_raw) + if integral_scale is not None: + self.integral_scale = integral_scale + # final check for parameter bounds self.check_arg_bounds() # additional checks for the optional arguments (provided by user) self.check_opt_arg() # precision for printing self._prec = 3 + # initialized + self._init = True # one of these functions needs to be overridden def __init_subclass__(cls): @@ -244,7 +244,8 @@ def __init_subclass__(cls): # modify the docstrings: class docstring gets attributes added if cls.__doc__ is None: cls.__doc__ = "User defined GSTools Covariance-Model." - cls.__doc__ += CovModel.__doc__[45:] + if cls._add_doc: + cls.__doc__ += CovModel.__doc__[45:] # overridden functions get standard doc if no new doc was created ign = ["__", "variogram", "covariance", "cor"] for att, attr_cls in cls.__dict__.items(): @@ -474,10 +475,6 @@ def fix_dim(self): """Set a fix dimension for the model.""" return None - def var_factor(self): - """Factor for the variance.""" - return 1.0 - def default_rescale(self): """Provide default rescaling factor.""" return 1.0 @@ -486,8 +483,7 @@ def default_rescale(self): def calc_integral_scale(self): """Calculate the integral scale of the isotrope model.""" - self._integral_scale = integral(self.correlation, 0, np.inf)[0] - return self._integral_scale + return integral(self.correlation, 0, np.inf)[0] def percentile_scale(self, per=0.9): """Calculate the percentile scale of the isotrope model. @@ -765,8 +761,8 @@ def default_arg_bounds(self): Given as a dictionary. """ res = { - "var": (0.0, np.inf, "oo"), - "len_scale": (0.0, np.inf, "oo"), + "var": (0.0, np.inf, "co"), + "len_scale": (0.0, np.inf, "co"), "nugget": (0.0, np.inf, "co"), "anis": (0.0, np.inf, "oo"), } @@ -809,11 +805,7 @@ def var_bounds(self): @var_bounds.setter def var_bounds(self, bounds): - if not check_bounds(bounds): - raise ValueError( - f"Given bounds for 'var' are not valid, got: {bounds}" - ) - self._var_bounds = bounds + self.set_arg_bounds(var=bounds) @property def len_scale_bounds(self): @@ -829,11 +821,7 @@ def len_scale_bounds(self): @len_scale_bounds.setter def len_scale_bounds(self, bounds): - if not check_bounds(bounds): - raise ValueError( - f"Given bounds for 'len_scale' are not valid, got: {bounds}" - ) - self._len_scale_bounds = bounds + self.set_arg_bounds(len_scale=bounds) @property def nugget_bounds(self): @@ -849,11 +837,7 @@ def nugget_bounds(self): @nugget_bounds.setter def nugget_bounds(self, bounds): - if not check_bounds(bounds): - raise ValueError( - f"Given bounds for 'nugget' are not valid, got: {bounds}" - ) - self._nugget_bounds = bounds + self.set_arg_bounds(nugget=bounds) @property def anis_bounds(self): @@ -869,11 +853,7 @@ def anis_bounds(self): @anis_bounds.setter def anis_bounds(self, bounds): - if not check_bounds(bounds): - raise ValueError( - f"Given bounds for 'anis' are not valid, got: {bounds}" - ) - self._anis_bounds = bounds + self.set_arg_bounds(anis=bounds) @property def opt_arg_bounds(self): @@ -949,24 +929,11 @@ def dim(self, dim): @property def var(self): """:class:`float`: The variance of the model.""" - return self._var * self.var_factor() + return self._var @var.setter def var(self, var): - self._var = float(var) / self.var_factor() - self.check_arg_bounds() - - @property - def var_raw(self): - """:class:`float`: The raw variance of the model without factor. - - (See. CovModel.var_factor) - """ - return self._var - - @var_raw.setter - def var_raw(self, var_raw): - self._var = float(var_raw) + self._var = float(var) self.check_arg_bounds() @property @@ -989,10 +956,7 @@ def len_scale(self, len_scale): self._len_scale, anis = set_len_anis( self.dim, len_scale, self.anis, self.latlon ) - if self.latlon: - self._anis = np.array((self.dim - 1) * [1], dtype=np.double) - else: - self._anis = anis + self._anis = anis self.check_arg_bounds() @property @@ -1008,7 +972,7 @@ def rescale(self, rescale): @property def len_rescaled(self): """:class:`float`: The rescaled main length scale of the model.""" - return self._len_scale / self._rescale + return self.len_scale / self.rescale @property def anis(self): @@ -1043,24 +1007,21 @@ def integral_scale(self): ValueError If integral scale is not setable. """ - self._integral_scale = self.calc_integral_scale() - return self._integral_scale + return self.calc_integral_scale() @integral_scale.setter def integral_scale(self, integral_scale): - if integral_scale is not None: - # format int-scale right - self.len_scale = integral_scale - integral_scale = self.len_scale - # reset len_scale - self.len_scale = 1.0 - int_tmp = self.calc_integral_scale() - self.len_scale = integral_scale / int_tmp - if not np.isclose(self.integral_scale, integral_scale, rtol=1e-3): - raise ValueError( - f"{self.name}: Integral scale could not be set correctly! " - "Please just provide a 'len_scale'!" - ) + int_scale, anis = set_len_anis( + self.dim, integral_scale, self.anis, self.latlon + ) + old_scale = self.integral_scale + self.anis = anis + self.len_scale = self.len_scale * int_scale / old_scale + if not np.isclose(self.integral_scale, int_scale, rtol=1e-3): + raise ValueError( + f"{self.name}: Integral scale could not be set correctly! " + "Please just provide a 'len_scale'!" + ) @property def hankel_kw(self): @@ -1115,7 +1076,7 @@ def sill(self): @property def arg(self): """:class:`list` of :class:`str`: Names of all arguments.""" - return ["var", "len_scale", "nugget", "anis", "angles"] + self._opt_arg + return ["var", "len_scale", "nugget", "anis", "angles"] + self.opt_arg @property def arg_list(self): @@ -1128,7 +1089,7 @@ def arg_list(self): @property def iso_arg(self): """:class:`list` of :class:`str`: Names of isotropic arguments.""" - return ["var", "len_scale", "nugget"] + self._opt_arg + return ["var", "len_scale", "nugget"] + self.opt_arg @property def iso_arg_list(self): @@ -1156,7 +1117,7 @@ def len_scale_vec(self): """ res = np.zeros(self.dim, dtype=np.double) res[0] = self.len_scale - for i in range(1, self._dim): + for i in range(1, self.dim): res[i] = self.len_scale * self.anis[i - 1] return res @@ -1202,9 +1163,548 @@ def __setattr__(self, name, value): """Set an attribute.""" super().__setattr__(name, value) # if an optional variogram argument was given, check bounds - if hasattr(self, "_opt_arg") and name in self._opt_arg: + if getattr(self, "_init", False) and name in self.opt_arg: self.check_arg_bounds() + def __add__(self, other): + """Add two covariance models and return a SumModel.""" + return SumModel(self, other) + + def __radd__(self, other): + """Add two covariance models and return a SumModel.""" + return SumModel(self, other) + def __repr__(self): """Return String representation.""" return model_repr(self) + + +class SumModel(CovModel): + r"""Class for sums of covariance models. + + This class represents sums of covariance models. The nugget of + each contained model will be set to zero and the sum model will + have its own nugget. + The variance of the sum model is the sum of the sub model variances + and the length scale of the sum model is the variance-weighted sum + of the length scales of the sub models. This is motivated by the fact, + that the integral scale of the sum model is equal to the variance-weighted + sum of the integral scales of the sub models. + An empty sum represents a pure Nugget model. + Resetting the total variance or the total length scale will evenly + scale the variances or length scales of the sub models. + Sum models will have a constant rescale factor of one. + + Parameters + ---------- + *models + tuple of :any:`CovModel` instances of subclasses to sum. + All models will get a nugget of zero and the nugget will + be set in the SumModel directly. + Models need to have matching temporal setting, + latlon setting, anis, angles and geo_scale. + The model dimension will be specified by the first given model. + dim : :class:`int`, optional + dimension of the model. + Includes the temporal dimension if temporal is true. + To specify only the spatial dimension in that case, use `spatial_dim`. + Default: ``3`` or dimension of the first given model (if instance). + vars : :class:`list` of :class:`float`, optional + variances of the models. Will reset present variances. + len_scales : :class:`list` of :class:`float`, optional + length scale of the models. Will reset present length scales. + nugget : :class:`float`, optional + nugget of the sum-model. All summed models will a nugget of zero. + Default: ``0.0`` + anis : :class:`float` or :class:`list`, optional + anisotropy ratios in the transversal directions [e_y, e_z]. + + * e_y = l_y / l_x + * e_z = l_z / l_x + + If only one value is given in 3D, e_y will be set to 1. + This value will be ignored, if multiple len_scales are given. + Default: ``1.0`` or anis of the first given model (if instance). + angles : :class:`float` or :class:`list`, optional + angles of rotation (given in rad): + + * in 2D: given as rotation around z-axis + * in 3D: given by yaw, pitch, and roll (known as Tait–Bryan angles) + + Default: ``0.0`` or angles of the first given model (if instance). + integral_scale : :class:`float` or :class:`list` or :any:`None`, optional + If given, ``len_scale`` will be ignored and recalculated, + so that the integral scale of the model matches the given one. + Default: :any:`None` + latlon : :class:`bool`, optional + Whether the model is describing 2D fields on earths surface described + by latitude and longitude. When using this, the model will internally + use the associated 'Yadrenko' model to represent a valid model. + This means, the spatial distance :math:`r` will be replaced by + :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle + distance, which is equal to the spatial distance of two points in 3D. + As a consequence, `dim` will be set to `3` and anisotropy will be + disabled. `geo_scale` can be set to e.g. earth's radius, + to have a meaningful `len_scale` parameter. + Default: False or latlon config of the first given model (if instance). + geo_scale : :class:`float`, optional + Geographic unit scaling in case of latlon coordinates to get a + meaningful length scale unit. + By default, len_scale is assumed to be in radians with latlon=True. + Can be set to :any:`KM_SCALE` to have len_scale in km or + :any:`DEGREE_SCALE` to have len_scale in degrees. + Default: :any:`RADIAN_SCALE` or geo_scale of the first given model (if instance). + temporal : :class:`bool`, optional + Create a metric spatio-temporal covariance model. + Setting this to true will increase `dim` and `field_dim` by 1. + `spatial_dim` will be `field_dim - 1`. + The time-dimension is appended, meaning the pos tuple is (x,y,z,...,t). + Default: False or temporal config of the first given model (if instance). + spatial_dim : :class:`int`, optional + spatial dimension of the model. + If given, the model dimension will be determined from this spatial dimension + and the possible temporal dimension if temporal is ture. + Default: None + var : :class:`float`, optional + Total variance of the sum-model. Will evenly scale present variances. + len_scale : :class:`float` or :class:`list`, optional + Total length scale of the sum-model. Will evenly scale present length scales. + **opt_arg + Optional arguments of the sub-models will have and added index of the sub-model. + Also covers ``var_`` and ``length_scale_`` but they should preferably be + set by ``vars`` and ``length_scales``. + """ + + _add_doc = False + + def __init__(self, *models, **kwargs): + self._init = False + self._models = [] + add_nug = 0.0 + to_init = None + imsg = ( + "SumModel: either all models are CovModel instances or subclasses." + ) + for mod in models: + if isinstance(mod, type) and issubclass(mod, SumModel): + if to_init is not None and not to_init: + raise ValueError(imsg) + to_init = True + continue # treat un-init sum-model as nugget model with 0 nugget + if isinstance(mod, SumModel): + if to_init is not None and to_init: + raise ValueError(imsg) + to_init = False + self._models += mod.models + add_nug += mod.nugget + elif isinstance(mod, CovModel): + if to_init is not None and to_init: + raise ValueError(imsg) + to_init = False + self._models.append(mod) + elif isinstance(mod, type) and issubclass(mod, CovModel): + if to_init is not None and not to_init: + raise ValueError(imsg) + to_init = True + self._models.append(mod(**default_mod_kwargs(kwargs))) + else: + msg = "SumModel: models need to be instances or subclasses of CovModel." + raise ValueError(msg) + # handle parameters when only Nugget models given + if models and not self.models: + for mod in models: + if not isinstance(mod, type): + kwargs.setdefault("dim", mod.dim) + kwargs.setdefault("latlon", mod.latlon) + kwargs.setdefault("temporal", mod.temporal) + kwargs.setdefault("geo_scale", mod.geo_scale) + kwargs.setdefault("anis", mod.anis) + kwargs.setdefault("angles", mod.angles) + break + # pop entries that can't be re-set + self._latlon = bool( + kwargs.pop( + "latlon", self.models[0].latlon if self.models else False + ) + ) + self._temporal = bool( + kwargs.pop( + "temporal", self.models[0].temporal if self.models else False + ) + ) + self._geo_scale = float( + kwargs.pop( + "geo_scale", + self.models[0].geo_scale if self.models else RADIAN_SCALE, + ) + ) + var_set = kwargs.pop("var", None) + len_set = kwargs.pop("len_scale", None) + # convert nugget + self._nugget = float( + kwargs.pop( + "nugget", + sum((mod.nugget for mod in self.models), 0) + add_nug, + ) + ) + for mod in self.models: + mod._nugget = 0.0 + # prepare dim setting + if "spatial_dim" in kwargs: + spatial_dim = kwargs.pop("spatial_dim") + if spatial_dim is not None: + kwargs["dim"] = spatial_dim + int(self.temporal) + self._dim = None + self._hankel_kw = {} + self._sft = None + self.dim = kwargs.get("dim", self.models[0].dim if self.models else 3) + # prepare parameters (they are checked in dim setting) + anis = kwargs.get("anis", self.models[0].anis if self.models else 1) + angles = kwargs.get( + "angles", self.models[0].angles if self.models else 0 + ) + _, self._anis = set_len_anis(self.dim, 1.0, anis, self.latlon) + self._angles = set_model_angles( + self.dim, angles, self.latlon, self.temporal + ) + # prepare parameter boundaries + self._var_bounds = None + self._len_scale_bounds = None + self._nugget_bounds = None + self._anis_bounds = None + self._opt_arg_bounds = {} + bounds = self.default_arg_bounds() + bounds.update(self.default_opt_arg_bounds()) + self.set_arg_bounds(check_args=False, **bounds) + # finalize init + self._prec = 3 + self._init = True + # set remaining args + for arg, val in kwargs.items(): + setattr(self, arg, val) + # reset total variance and length scale last + if var_set is not None: + self.var = var_set + if len_set is not None: + self.len_scale = len_set + # check consistency of sub models + self.check() + + def __iter__(self): + return iter(self.models) + + def __len__(self): + return len(self.models) + + def __contains__(self, item): + return item in self.models + + def __getitem__(self, key): + return self.models[key] + + def check(self): + """Check consistency of contained models.""" + sum_check(self) + + def default_arg_bounds(self): + """Default boundaries for arguments as dict.""" + return sum_default_arg_bounds(self) + + def default_opt_arg_bounds(self): + """Defaults boundaries for optional arguments as dict.""" + return sum_default_opt_arg_bounds(self) + + def set_norm_var_ratios(self, ratios, skip=None, var=None): + """ + Set variances of contained models by normalized ratios in [0, 1]. + + Ratios are given as normalized ratios in [0, 1] as relative ratio of + variance to remaining difference to total variance of the Sum-Model. + + Parameters + ---------- + ratios : iterable + Ratios to set. Should have a length of len(models) - len(exclude) - 1 + skip : iterable, optional + Model indices to skip. Should have compatible lenth, by default None + var : float, optional + Desired variance, by default current variance + + Raises + ------ + ValueError + If number of ratios is not matching. + """ + sum_set_norm_var_ratios(self, ratios, skip, var) + + def set_norm_len_ratios(self, ratios, skip=None, len_scale=None): + """ + Set length scales of contained models by normalized ratios in [0, 1]. + + Ratios are given as normalized ratios in [0, 1] as relative ratio of + len_scale * var / total_var to remaining difference to + total len_scale of the Sum-Model. + + Parameters + ---------- + ratios : iterable + Ratios to set. Should have a length of len(models) - len(exclude) - 1 + skip : iterable, optional + Model indices to skip. Should have compatible lenth, by default None + len_scale : float, optional + Desired len_scale, by default current len_scale + + Raises + ------ + ValueError + If number of ratios is not matching. + """ + sum_set_norm_len_ratios(self, ratios, skip, len_scale) + + @property + def hankel_kw(self): + """:class:`dict`: :any:`hankel.SymmetricFourierTransform` kwargs.""" + return self._hankel_kw + + @property + def models(self): + """:class:`tuple`: The summed models.""" + return self._models + + @property + def rescale(self): + """:class:`float`: SumModel has a constant rescale factor of one.""" + return 1.0 + + @property + def geo_scale(self): + """:class:`float`: Geographic scaling for geographical coords.""" + return self._geo_scale + + @geo_scale.setter + def geo_scale(self, geo_scale): + self._geo_scale = abs(float(geo_scale)) + for mod in self.models: + mod.geo_scale = geo_scale + + @property + def dim(self): + """:class:`int`: The dimension of the model.""" + return self._dim + + @dim.setter + def dim(self, dim): + set_dim(self, dim) + for mod in self.models: + mod.dim = dim + + @property + def var(self): + """:class:`float`: The variance of the model.""" + return sum((var for var in self.vars), 0.0) + + @var.setter + def var(self, var): + for mod, rat in zip(self.models, self.ratios): + mod.var = rat * var + if not self.models and not np.isclose(var, 0): + msg = f"{self.name}: variance needs to be 0." + raise ValueError(msg) + check_arg_in_bounds(self, "var", error=True) + check_arg_in_bounds(self, "len_scale", error=True) + + @property + def vars(self): + """:class:`list`: The variances of the models.""" + return [mod.var for mod in self.models] + + @vars.setter + def vars(self, variances): + if len(variances) != len(self): + msg = "SumModel: number of given variances not matching" + raise ValueError(msg) + for mod, var in zip(self.models, variances): + mod.var = var + check_arg_in_bounds(self, "var", error=True) + check_arg_in_bounds(self, "len_scale", error=True) + + @property + def len_scale(self): + """:class:`float`: The main length scale of the model.""" + return sum( + ( + mod.len_scale * rat + for mod, rat in zip(self.models, self.ratios) + ), + 0.0, + ) + + @len_scale.setter + def len_scale(self, len_scale): + len_scale, anis = set_len_anis( + self.dim, len_scale, self.anis, self.latlon + ) + old_scale = self.len_scale + self.anis = anis + for mod in self.models: + mod.len_scale = mod.len_scale * len_scale / old_scale + if not self.models and not np.isclose(len_scale, 0): + msg = f"{self.name}: length scale needs to be 0." + raise ValueError(msg) + check_arg_in_bounds(self, "len_scale", error=True) + + @property + def len_scales(self): + """:class:`list`: The main length scales of the models.""" + return [mod.len_scale for mod in self.models] + + @len_scales.setter + def len_scales(self, len_scales): + if len(len_scales) != len(self): + msg = "SumModel: number of given length scales not matching" + raise ValueError(msg) + for mod, len_scale in zip(self.models, len_scales): + mod.len_scale = len_scale + check_arg_in_bounds(self, "len_scale", error=True) + + @property + def anis(self): + """:class:`numpy.ndarray`: The anisotropy factors of the model.""" + return self._anis + + @anis.setter + def anis(self, anis): + _, self._anis = set_len_anis(self.dim, 1.0, anis, self.latlon) + for mod in self.models: + mod.anis = anis + check_arg_in_bounds(self, "anis", error=True) + + @property + def angles(self): + """:class:`numpy.ndarray`: Rotation angles (in rad) of the model.""" + return self._angles + + @angles.setter + def angles(self, angles): + self._angles = set_model_angles( + self.dim, angles, self.latlon, self.temporal + ) + for mod in self.models: + mod.angles = angles + + @property + def ratios(self): + """:class:`numpy.ndarray`: Variance ratios of the sub-models.""" + var = self.var + if np.isclose(var, 0) and len(self) > 0: + return np.full(len(self), 1 / len(self)) + return np.array([mod.var / var for mod in self.models]) + + @ratios.setter + def ratios(self, ratios): + if len(ratios) != len(self): + msg = "SumModel.ratios: wrong number of given ratios." + raise ValueError(msg) + if ratios and not np.isclose(np.sum(ratios), 1): + msg = "SumModel.ratios: given ratios do not sum to 1." + raise ValueError(msg) + var = self.var + for mod, rat in zip(self.models, ratios): + mod.var = var * rat + check_arg_in_bounds(self, "var", error=True) + check_arg_in_bounds(self, "len_scale", error=True) + + def calc_integral_scale(self): + return sum( + ( + mod.integral_scale * rat + for mod, rat in zip(self.models, self.ratios) + ), + 0.0, + ) + + def spectral_density(self, k): + return sum( + ( + mod.spectral_density(k) * rat + for mod, rat in zip(self.models, self.ratios) + ), + np.zeros_like(k, dtype=np.double), + ) + + def correlation(self, r): + """SumModel correlation function.""" + return sum( + ( + mod.correlation(r) * rat + for mod, rat in zip(self.models, self.ratios) + ), + np.zeros_like(r, dtype=np.double), + ) + + @property + def opt_arg(self): + """:class:`list` of :class:`str`: Names of the optional arguments.""" + return sum( + [ + [f"{opt}_{i}" for opt in mod.opt_arg] + for i, mod in enumerate(self.models) + ], + [], + ) + + @property + def sub_arg(self): + """:class:`list` of :class:`str`: Names of the sub-arguments for var and len_scale.""" + return sum( + [ + [f"{o}_{i}" for o in ["var", "len_scale"]] + for i, mod in enumerate(self.models) + ], + [], + ) + + def __setattr__(self, name, value): + """Set an attribute.""" + sub_arg = False + if getattr(self, "_init", False): + for i, mod in enumerate(self.models): + if name == f"var_{i}": + mod.var = value + sub_arg = True + if name == f"len_scale_{i}": + mod.len_scale = value + sub_arg = True + if name == f"integral_scale_{i}": + mod.integral_scale = value + sub_arg = True + for opt in mod.opt_arg: + if name == f"{opt}_{i}": + setattr(mod, opt, value) + sub_arg = True + if sub_arg: + break + if sub_arg: + self.check_arg_bounds() + else: + super().__setattr__(name, value) + + def __getattr__(self, name): + """Get an attribute.""" + # __getattr__ is only called when an attribute is not found in the usual places + if name != "_init" and getattr(self, "_init", False): + for i, mod in enumerate(self.models): + if name == f"var_{i}": + return mod.var + if name == f"len_scale_{i}": + return mod.len_scale + if name == f"integral_scale_{i}": + return mod.integral_scale + for opt in mod.opt_arg: + if name == f"{opt}_{i}": + return getattr(mod, opt) + raise AttributeError(f"'{self.name}' object has no attribute '{name}'") + + def __repr__(self): + """Return String representation.""" + return sum_model_repr(self) diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 8b19f497..5ba398ef 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -19,9 +19,6 @@ __all__ = ["fit_variogram"] -DEFAULT_PARA = ["var", "len_scale", "nugget"] - - def fit_variogram( model, x_data, @@ -169,8 +166,10 @@ def fit_variogram( para, sill, constrain_sill, anis = _pre_para( model, para_select, sill, anis ) + if not any(para.values()): + raise ValueError("fit: no parameters selected for fitting.") # check curve_fit kwargs - curve_fit_kwargs = {} if curve_fit_kwargs is None else curve_fit_kwargs + curve_fit_kwargs = curve_fit_kwargs or {} # check method if method not in ["trf", "dogbox"]: raise ValueError("fit: method needs to be either 'trf' or 'dogbox'") @@ -213,21 +212,13 @@ def fit_variogram( def _pre_para(model, para_select, sill, anis): """Preprocess selected parameters.""" - var_last = False - var_tmp = 0.0 # init value + # if values given, set them in the model, afterwards all entries are bool for par in para_select: if par not in model.arg_bounds: raise ValueError(f"fit: unknown parameter in selection: {par}") if not isinstance(para_select[par], bool): - if par == "var": - var_last = True - var_tmp = float(para_select[par]) - else: - setattr(model, par, float(para_select[par])) + setattr(model, par, float(para_select[par])) para_select[par] = False - # set variance last due to possible recalculations - if var_last: - model.var = var_tmp # remove those that were set to True para_select = {k: v for k, v in para_select.items() if not v} # handling the sill @@ -268,8 +259,7 @@ def _pre_para(model, para_select, sill, anis): else: constrain_sill = False # select all parameters to be fitted - para = {par: True for par in DEFAULT_PARA} - para.update({opt: True for opt in model.opt_arg}) + para = {par: True for par in model.iso_arg} # now deselect unwanted parameters para.update(para_select) # check if anisotropy should be fitted or set @@ -368,7 +358,7 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): low_bounds = [] top_bounds = [] init_guess_list = [] - for par in DEFAULT_PARA: + for par in model.iso_arg: if para[par]: low_bounds.append(model.arg_bounds[par][0]) if par == "var" and constrain_sill: # var <= sill in this case @@ -381,16 +371,6 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): default=init_guess[par], ) ) - for opt in model.opt_arg: - if para[opt]: - low_bounds.append(model.arg_bounds[opt][0]) - top_bounds.append(model.arg_bounds[opt][1]) - init_guess_list.append( - _init_guess( - bounds=[low_bounds[-1], top_bounds[-1]], - default=init_guess[opt], - ) - ) if anis: for i in range(model.dim - 1): low_bounds.append(model.anis_bounds[0]) @@ -413,40 +393,23 @@ def _init_guess(bounds, default): def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): """Create the curve for scipys curve_fit.""" - var_save = model.var # we need arg1, otherwise curve_fit throws an error (bug?!) def curve(x, arg1, *args): """Adapted Variogram function.""" args = (arg1,) + args para_skip = 0 - opt_skip = 0 - if para["var"]: - var_tmp = args[para_skip] - if constrain_sill: - nugget_tmp = sill - var_tmp - # punishment, if resulting nugget out of range for fixed sill - if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0: - return np.full_like(x, np.inf) - # nugget estimation deselected in this case - model.nugget = nugget_tmp - para_skip += 1 - if para["len_scale"]: - model.len_scale = args[para_skip] - para_skip += 1 - if para["nugget"]: - model.nugget = args[para_skip] - para_skip += 1 - for opt in model.opt_arg: - if para[opt]: - setattr(model, opt, args[para_skip + opt_skip]) - opt_skip += 1 - # set var at last because of var_factor (other parameter needed) - if para["var"]: - model.var = var_tmp - # needs to be reset for TPL models when len_scale was changed - else: - model.var = var_save + for par in model.iso_arg: + if para[par]: + setattr(model, par, args[para_skip]) + para_skip += 1 + if constrain_sill: + nugget_tmp = sill - model.var + # punishment, if resulting nugget out of range for fixed sill + if check_arg_in_bounds(model, "nugget", nugget_tmp) > 0: + return np.full_like(x, np.inf) + # nugget estimation deselected in this case + model.nugget = nugget_tmp if is_dir_vario: if anis: model.anis = args[1 - model.dim :] @@ -464,32 +427,17 @@ def _post_fitting(model, para, popt, anis, is_dir_vario): """Postprocess fitting results and application to model.""" fit_para = {} para_skip = 0 - opt_skip = 0 - var_tmp = 0.0 # init value - for par in DEFAULT_PARA: + for par in model.iso_arg: if para[par]: - if par == "var": # set variance last - var_tmp = popt[para_skip] - else: - setattr(model, par, popt[para_skip]) + setattr(model, par, popt[para_skip]) fit_para[par] = popt[para_skip] para_skip += 1 else: fit_para[par] = getattr(model, par) - for opt in model.opt_arg: - if para[opt]: - setattr(model, opt, popt[para_skip + opt_skip]) - fit_para[opt] = popt[para_skip + opt_skip] - opt_skip += 1 - else: - fit_para[opt] = getattr(model, opt) if is_dir_vario: if anis: model.anis = popt[1 - model.dim :] fit_para["anis"] = model.anis - # set var at last because of var_factor (other parameter needed) - if para["var"]: - model.var = var_tmp return fit_para diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index b1a9d68e..7c1cd836 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -6,6 +6,7 @@ The following classes are provided .. autosummary:: + Nugget Gaussian Exponential Matern @@ -21,17 +22,18 @@ JBessel """ -# pylint: disable=C0103, E1101, R0201 +# pylint: disable=C0103, C0302, E1101, R0201 import warnings import numpy as np from scipy import special as sps -from gstools.covmodel.base import CovModel +from gstools.covmodel.base import CovModel, SumModel from gstools.covmodel.tools import AttributeWarning from gstools.tools.special import exp_int, inc_gamma_low __all__ = [ + "Nugget", "Gaussian", "Exponential", "Matern", @@ -48,6 +50,69 @@ ] +class Nugget(SumModel): + r"""Pure nugget model. + + Parameters + ---------- + dim : :class:`int`, optional + dimension of the model. + Includes the temporal dimension if temporal is true. + To specify only the spatial dimension in that case, use `spatial_dim`. + Default: ``3`` + nugget : :class:`float`, optional + nugget of the model. Default: ``0.0`` + anis : :class:`float` or :class:`list`, optional + anisotropy ratios in the transversal directions [e_y, e_z]. + + * e_y = l_y / l_x + * e_z = l_z / l_x + + If only one value is given in 3D, e_y will be set to 1. + This value will be ignored, if multiple len_scales are given. + Default: ``1.0`` + angles : :class:`float` or :class:`list`, optional + angles of rotation (given in rad): + + * in 2D: given as rotation around z-axis + * in 3D: given by yaw, pitch, and roll (known as Tait–Bryan angles) + + Default: ``0.0`` + latlon : :class:`bool`, optional + Whether the model is describing 2D fields on earths surface described + by latitude and longitude. When using this, the model will internally + use the associated 'Yadrenko' model to represent a valid model. + This means, the spatial distance :math:`r` will be replaced by + :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle + distance, which is equal to the spatial distance of two points in 3D. + As a consequence, `dim` will be set to `3` and anisotropy will be + disabled. `geo_scale` can be set to e.g. earth's radius, + to have a meaningful `len_scale` parameter. + Default: False + geo_scale : :class:`float`, optional + Geographic unit scaling in case of latlon coordinates to get a + meaningful length scale unit. + By default, len_scale is assumed to be in radians with latlon=True. + Can be set to :any:`KM_SCALE` to have len_scale in km or + :any:`DEGREE_SCALE` to have len_scale in degrees. + Default: :any:`RADIAN_SCALE` + temporal : :class:`bool`, optional + Create a metric spatio-temporal covariance model. + Setting this to true will increase `dim` and `field_dim` by 1. + `spatial_dim` will be `field_dim - 1`. + The time-dimension is appended, meaning the pos tuple is (x,y,z,...,t). + Default: False + spatial_dim : :class:`int`, optional + spatial dimension of the model. + If given, the model dimension will be determined from this spatial dimension + and the possible temporal dimension if temporal is ture. + Default: None + """ + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + class Gaussian(CovModel): r"""The Gaussian covariance model. diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index 32148c14..f7fa58db 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -68,7 +68,10 @@ def plot_vario_spatial( ): # pragma: no cover """Plot spatial variogram of a given CovModel.""" if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(-x_max, x_max) + x_min pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) @@ -83,7 +86,10 @@ def plot_cov_spatial( ): # pragma: no cover """Plot spatial covariance of a given CovModel.""" if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(-x_max, x_max) + x_min pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) @@ -98,7 +104,10 @@ def plot_cor_spatial( ): # pragma: no cover """Plot spatial correlation of a given CovModel.""" if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(-x_max, x_max) + x_min pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) @@ -114,7 +123,10 @@ def plot_variogram( """Plot variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} variogram") ax.plot(x_s, model.variogram(x_s), **kwargs) @@ -129,7 +141,10 @@ def plot_covariance( """Plot covariance of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} covariance") ax.plot(x_s, model.covariance(x_s), **kwargs) @@ -144,7 +159,10 @@ def plot_correlation( """Plot correlation function of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} correlation") ax.plot(x_s, model.correlation(x_s), **kwargs) @@ -159,7 +177,10 @@ def plot_vario_yadrenko( """Plot Yadrenko variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_scale, model.geo_scale * np.pi) + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko variogram") ax.plot(x_s, model.vario_yadrenko(x_s), **kwargs) @@ -174,7 +195,10 @@ def plot_cov_yadrenko( """Plot Yadrenko covariance of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_scale, model.geo_scale * np.pi) + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko covariance") ax.plot(x_s, model.cov_yadrenko(x_s), **kwargs) @@ -189,7 +213,10 @@ def plot_cor_yadrenko( """Plot Yadrenko correlation function of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_scale, model.geo_scale * np.pi) + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko correlation") ax.plot(x_s, model.cor_yadrenko(x_s), **kwargs) @@ -204,7 +231,10 @@ def plot_vario_axis( """Plot variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} variogram on axis {axis}") ax.plot(x_s, model.vario_axis(x_s, axis), **kwargs) @@ -219,7 +249,10 @@ def plot_cov_axis( """Plot variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} covariance on axis {axis}") ax.plot(x_s, model.cov_axis(x_s, axis), **kwargs) @@ -234,7 +267,10 @@ def plot_cor_axis( """Plot variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 * model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 * model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} correlation on axis {axis}") ax.plot(x_s, model.cor_axis(x_s, axis), **kwargs) @@ -249,7 +285,10 @@ def plot_spectrum( """Plot spectrum of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 / model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} {model.dim}D spectrum") ax.plot(x_s, model.spectrum(x_s), **kwargs) @@ -264,7 +303,10 @@ def plot_spectral_density( """Plot spectral density of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 / model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} {model.dim}D spectral-density") ax.plot(x_s, model.spectral_density(x_s), **kwargs) @@ -279,7 +321,10 @@ def plot_spectral_rad_pdf( """Plot radial spectral pdf of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = 3 / model.len_scale + if np.isclose(model.len_scale, 0): + x_max = 1.0 + else: + x_max = 3 / model.len_scale x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} {model.dim}D spectral-rad-pdf") ax.plot(x_s, model.spectral_rad_pdf(x_s), **kwargs) diff --git a/src/gstools/covmodel/sum_tools.py b/src/gstools/covmodel/sum_tools.py new file mode 100644 index 00000000..c4b9a90e --- /dev/null +++ b/src/gstools/covmodel/sum_tools.py @@ -0,0 +1,246 @@ +""" +GStools subpackage providing tools for sum-models. + +.. currentmodule:: gstools.covmodel.sum_tools + +The following classes and functions are provided + +.. autosummary:: + RatioError + ARG_DEF + default_mod_kwargs + sum_check + sum_default_arg_bounds + sum_default_opt_arg_bounds + sum_set_norm_var_ratios + sum_set_norm_len_ratios + sum_model_repr +""" + +# pylint: disable=W0212 +import numpy as np + +from gstools.covmodel.tools import check_arg_in_bounds +from gstools.tools import RADIAN_SCALE +from gstools.tools.misc import list_format + +__all__ = [ + "RatioError", + "ARG_DEF", + "default_mod_kwargs", + "sum_check", + "sum_default_arg_bounds", + "sum_default_opt_arg_bounds", + "sum_set_norm_var_ratios", + "sum_set_norm_len_ratios", + "sum_model_repr", +] + + +class RatioError(Exception): + """Error for invalid ratios in SumModel.""" + + +ARG_DEF = { + "dim": 3, + "latlon": False, + "temporal": False, + "geo_scale": RADIAN_SCALE, + "spatial_dim": None, + "hankel_kw": None, +} +"""dict: default model arguments""" + + +def default_mod_kwargs(kwargs): + """Generate default model keyword arguments.""" + mod_kw = {} + for arg, default in ARG_DEF.items(): + mod_kw[arg] = kwargs.get(arg, default) + return mod_kw + + +def sum_check(summod): + """Check consistency of contained models.""" + # prevent dim error in anis and angles + if any(mod.dim != summod.dim for mod in summod): + msg = "SumModel: models need to have same dimension." + raise ValueError(msg) + if any(mod.latlon != summod.latlon for mod in summod): + msg = "SumModel: models need to have same latlon config." + raise ValueError(msg) + if any(mod.temporal != summod.temporal for mod in summod): + msg = "SumModel: models need to have same temporal config." + raise ValueError(msg) + if not all(np.isclose(mod.nugget, 0) for mod in summod): + msg = "SumModel: models need to have 0 nugget." + raise ValueError(msg) + if not np.allclose([mod.geo_scale for mod in summod], summod.geo_scale): + msg = "SumModel: models need to have same geo_scale." + raise ValueError(msg) + if not all(np.allclose(mod.anis, summod.anis) for mod in summod): + msg = "SumModel: models need to have same anisotropy ratios." + raise ValueError(msg) + if not all(np.allclose(mod.angles, summod.angles) for mod in summod): + msg = "SumModel: models need to have same rotation angles." + raise ValueError(msg) + + +def sum_default_arg_bounds(summod): + """Default boundaries for arguments as dict.""" + var_bnds = [mod.var_bounds for mod in summod.models] + len_bnds = [mod.len_scale_bounds for mod in summod.models] + var_lo = sum((bnd[0] for bnd in var_bnds), start=0.0) + var_hi = sum((bnd[1] for bnd in var_bnds), start=0.0) + len_lo = min((bnd[0] for bnd in len_bnds), default=0.0) + len_hi = max((bnd[1] for bnd in len_bnds), default=0.0) + res = { + "var": (var_lo, var_hi), + "len_scale": (len_lo, len_hi), + "nugget": (0.0, np.inf, "co"), + "anis": (0.0, np.inf, "oo"), + } + return res + + +def sum_default_opt_arg_bounds(summod): + """Defaults boundaries for optional arguments as dict.""" + bounds = {} + for i, mod in enumerate(summod.models): + bounds.update( + {f"{opt}_{i}": bnd for opt, bnd in mod.opt_arg_bounds.items()} + ) + return bounds + + +def sum_set_norm_var_ratios(summod, ratios, skip=None, var=None): + """ + Set variances of contained models by normalized ratios in [0, 1]. + + Ratios are given as normalized ratios in [0, 1] as relative ratio of + variance to remaining difference to total variance of the Sum-Model. + + Parameters + ---------- + ratios : iterable + Ratios to set. Should have a length of len(models) - len(exclude) - 1 + skip : iterable, optional + Model indices to skip. Should have compatible lenth, by default None + var : float, optional + Desired variance, by default current variance + + Raises + ------ + ValueError + If number of ratios is not matching. + """ + skip = skip or set() + if len(summod) != len(ratios) + len(skip) + 1: + msg = "SumModel.set_norm_ratios: number of ratios not matching." + raise ValueError(msg) + ids = range(len(summod)) + if fail := set(skip) - set(ids): + msg = f"SumModel.set_norm_var_ratios: skip ids not valid: {fail}" + raise ValueError(msg) + var = summod.var if var is None else float(var) + check_arg_in_bounds(summod, "var", var, error=True) + var_sum = sum(summod.models[i].var for i in skip) + if var_sum > var: + msg = "SumModel.set_norm_var_ratios: skiped variances to big." + raise RatioError(msg) + j = 0 + for i in ids: + if i in skip: + continue + var_diff = var - var_sum + # last model gets remaining diff + var_set = var_diff * ratios[j] if j < len(ratios) else var_diff + summod[i].var = var_set + var_sum += var_set + j += 1 + + +def sum_set_norm_len_ratios(summod, ratios, skip=None, len_scale=None): + """ + Set length scales of contained models by normalized ratios in [0, 1]. + + Ratios are given as normalized ratios in [0, 1] as relative ratio of + len_scale * var / total_var to remaining difference to + total len_scale of the Sum-Model. + + Parameters + ---------- + ratios : iterable + Ratios to set. Should have a length of len(models) - len(exclude) - 1 + skip : iterable, optional + Model indices to skip. Should have compatible lenth, by default None + len_scale : float, optional + Desired len_scale, by default current len_scale + + Raises + ------ + ValueError + If number of ratios is not matching. + """ + skip = skip or set() + if len(summod) != len(ratios) + len(skip) + 1: + msg = "SumModel.set_norm_len_ratios: number of ratios not matching." + raise ValueError(msg) + ids = range(len(summod)) + if fail := set(skip) - set(ids): + msg = f"SumModel.set_norm_len_ratios: skip ids not valid: {fail}" + raise ValueError(msg) + len_scale = summod.len_scale if len_scale is None else float(len_scale) + check_arg_in_bounds(summod, "len_scale", len_scale, error=True) + len_sum = sum(summod[i].len_scale * summod.ratios[i] for i in skip) + if len_sum > len_scale: + msg = "SumModel.set_norm_len_ratios: skiped length scales to big." + raise RatioError(msg) + j = 0 + for i in ids: + if i in skip: + continue + len_diff = len_scale - len_sum + # last model gets remaining diff + len_set = len_diff * ratios[j] if j < len(ratios) else len_diff + summod[i].len_scale = ( + 0.0 + if np.isclose(summod.ratios[j], 0) + else len_set / summod.ratios[j] + ) + len_sum += len_set + j += 1 + + +def sum_model_repr(summod): # pragma: no cover + """ + Generate the sum-model string representation. + + Parameters + ---------- + model : :any:`SumModel` + The sum-model in use. + """ + m, p = summod, summod._prec + ani_str, ang_str, o_str, r_str, p_str = "", "", "", "", "" + m_str = ", ".join([mod.name for mod in m.models]) + t_str = ", temporal=True" if m.temporal else "" + d_str = f"latlon={m.latlon}" if m.latlon else f"dim={m.spatial_dim}" + if len(m) > 0: + m_str += ", " + p_str += f", vars={list_format(m.vars, p)}" + p_str += f", len_scales={list_format(m.len_scales, p)}" + p_str += "" if np.isclose(m.nugget, 0) else f", nugget={m.nugget:.{p}}" + for opt in m.opt_arg: + o_str += f", {opt}={getattr(m, opt):.{p}}" + if m.latlon: + if not m.is_isotropic and m.temporal: + ani_str = f", anis={m.anis[-1]:.{p}}" + if not np.isclose(m.geo_scale, 1): + r_str = f", geo_scale={m.geo_scale:.{p}}" + else: + if not m.is_isotropic: + ani_str = f", anis={list_format(m.anis, p)}" + if m.do_rotation: + ang_str = f", angles={list_format(m.angles, p)}" + return f"{m.name}({m_str}{d_str}{t_str}{p_str}{ani_str}{ang_str}{r_str}{o_str})" diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index dddeb441..69bf963b 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -297,18 +297,21 @@ def check_bounds(bounds): * "oo" : open - open * "oc" : open - close * "co" : close - open - * "cc" : close - close + * "cc" : close - close (default) """ + typ = bounds[2] if len(bounds) == 3 else "cc" if len(bounds) not in (2, 3): return False - if bounds[1] <= bounds[0]: + if (typ == "cc" and bounds[1] < bounds[0]) or ( + typ != "cc" and bounds[1] <= bounds[0] + ): return False if len(bounds) == 3 and bounds[2] not in ("oo", "oc", "co", "cc"): return False return True -def check_arg_in_bounds(model, arg, val=None): +def check_arg_in_bounds(model, arg, val=None, error=False): """Check if given argument value is in bounds of the given model.""" if arg not in model.arg_bounds: raise ValueError(f"check bounds: unknown argument: {arg}") @@ -317,6 +320,7 @@ def check_arg_in_bounds(model, arg, val=None): val = np.asarray(val) error_case = 0 if len(bnd) == 2: + bnd = list(bnd) bnd.append("cc") # use closed intervals by default if bnd[2][0] == "c": if np.any(val < bnd[0]): @@ -330,6 +334,15 @@ def check_arg_in_bounds(model, arg, val=None): else: if np.any(val >= bnd[1]): error_case = 4 + if error: + if error_case == 1: + raise ValueError(f"{arg} needs to be >= {bnd[0]}, got: {val}") + if error_case == 2: + raise ValueError(f"{arg} needs to be > {bnd[0]}, got: {val}") + if error_case == 3: + raise ValueError(f"{arg} needs to be <= {bnd[1]}, got: {val}") + if error_case == 4: + raise ValueError(f"{arg} needs to be < {bnd[1]}, got: {val}") return error_case @@ -452,24 +465,20 @@ def set_arg_bounds(model, check_args=True, **kwargs): is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ - # if variance needs to be resetted, do this at last - var_bnds = [] for arg, bounds in kwargs.items(): if not check_bounds(bounds): - raise ValueError( - f"Given bounds for '{arg}' are not valid, got: {bounds}" - ) + msg = f"Given bounds for '{arg}' are not valid, got: {bounds}" + raise ValueError(msg) if arg in model.opt_arg: model._opt_arg_bounds[arg] = bounds elif arg == "var": - var_bnds = bounds - continue + model._var_bounds = bounds elif arg == "len_scale": - model.len_scale_bounds = bounds + model._len_scale_bounds = bounds elif arg == "nugget": - model.nugget_bounds = bounds + model._nugget_bounds = bounds elif arg == "anis": - model.anis_bounds = bounds + model._anis_bounds = bounds else: raise ValueError(f"set_arg_bounds: unknown argument '{arg}'") if check_args and check_arg_in_bounds(model, arg) > 0: @@ -478,11 +487,6 @@ def set_arg_bounds(model, check_args=True, **kwargs): setattr(model, arg, [def_arg] * (model.dim - 1)) else: setattr(model, arg, def_arg) - # set var last like always - if var_bnds: - model.var_bounds = var_bnds - if check_args and check_arg_in_bounds(model, "var") > 0: - model.var = default_arg_from_bounds(var_bnds) def check_arg_bounds(model): @@ -501,19 +505,7 @@ def check_arg_bounds(model): """ # check var, len_scale, nugget and optional-arguments for arg in model.arg_bounds: - if not model.arg_bounds[arg]: - continue # no bounds given during init (called from self.dim) - bnd = list(model.arg_bounds[arg]) - val = getattr(model, arg) - error_case = check_arg_in_bounds(model, arg) - if error_case == 1: - raise ValueError(f"{arg} needs to be >= {bnd[0]}, got: {val}") - if error_case == 2: - raise ValueError(f"{arg} needs to be > {bnd[0]}, got: {val}") - if error_case == 3: - raise ValueError(f"{arg} needs to be <= {bnd[1]}, got: {val}") - if error_case == 4: - raise ValueError(f"{arg} needs to be < {bnd[1]}, got: {val}") + check_arg_in_bounds(model, arg, error=True) def set_dim(model, dim): @@ -557,16 +549,15 @@ def set_dim(model, dim): model._dim = int(dim) # create fourier transform just once (recreate for dim change) model._sft = SFT(ndim=model.dim, **model.hankel_kw) - # recalculate dimension related parameters - if model._anis is not None: - model._len_scale, model._anis = set_len_anis( - model.dim, model._len_scale, model._anis + # recalculate dimension related parameters (if model initialized) + if model._init: + model.len_scale, model.anis = set_len_anis( + model.dim, model.len_scale, model.anis ) - if model._angles is not None: - model._angles = set_model_angles( - model.dim, model._angles, model.latlon, model.temporal + model.angles = set_model_angles( + model.dim, model.angles, model.latlon, model.temporal ) - model.check_arg_bounds() + model.check_arg_bounds() def compare(this, that): @@ -587,12 +578,12 @@ def compare(this, that): equal = True equal &= this.name == that.name equal &= np.isclose(this.var, that.var) - equal &= np.isclose(this.var_raw, that.var_raw) # ?! needless? equal &= np.isclose(this.nugget, that.nugget) equal &= np.isclose(this.len_scale, that.len_scale) equal &= np.all(np.isclose(this.anis, that.anis)) equal &= np.all(np.isclose(this.angles, that.angles)) equal &= np.isclose(this.rescale, that.rescale) + equal &= np.isclose(this.geo_scale, that.geo_scale) equal &= this.latlon == that.latlon equal &= this.temporal == that.temporal for opt in this.opt_arg: @@ -609,39 +600,24 @@ def model_repr(model): # pragma: no cover model : :any:`CovModel` The covariance model in use. """ - m = model - p = model._prec - opt_str = "" + m, p = model, model._prec + ani_str, ang_str, o_str, r_str = "", "", "", "" t_str = ", temporal=True" if m.temporal else "" + d_str = f"latlon={m.latlon}" if m.latlon else f"dim={m.spatial_dim}" + p_str = f", var={m.var:.{p}}, len_scale={m.len_scale:.{p}}" + p_str += "" if np.isclose(m.nugget, 0) else f", nugget={m.nugget:.{p}}" if not np.isclose(m.rescale, m.default_rescale()): - opt_str += f", rescale={m.rescale:.{p}}" + o_str += f", rescale={m.rescale:.{p}}" for opt in m.opt_arg: - opt_str += f", {opt}={getattr(m, opt):.{p}}" + o_str += f", {opt}={getattr(m, opt):.{p}}" if m.latlon: - ani_str = ( - "" - if m.is_isotropic or not m.temporal - else f", anis={m.anis[-1]:.{p}}" - ) - r_str = ( - "" - if np.isclose(m.geo_scale, 1) - else f", geo_scale={m.geo_scale:.{p}}" - ) - repr_str = ( - f"{m.name}(latlon={m.latlon}{t_str}, var={m.var:.{p}}, " - f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" - f"{ani_str}{r_str}{opt_str})" - ) + if not m.is_isotropic and m.temporal: + ani_str = f", anis={m.anis[-1]:.{p}}" + if not np.isclose(m.geo_scale, 1): + r_str = f", geo_scale={m.geo_scale:.{p}}" else: - # only print anis and angles if model is anisotropic or rotated - ani_str = "" if m.is_isotropic else f", anis={list_format(m.anis, p)}" - ang_str = ( - f", angles={list_format(m.angles, p)}" if m.do_rotation else "" - ) - repr_str = ( - f"{m.name}(dim={m.spatial_dim}{t_str}, var={m.var:.{p}}, " - f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" - f"{ani_str}{ang_str}{opt_str})" - ) - return repr_str + if not m.is_isotropic: + ani_str = f", anis={list_format(m.anis, p)}" + if m.do_rotation: + ang_str = f", angles={list_format(m.angles, p)}" + return f"{m.name}({d_str}{t_str}{p_str}{ani_str}{ang_str}{r_str}{o_str})" diff --git a/src/gstools/covmodel/tpl_models.py b/src/gstools/covmodel/tpl_models.py index b728e7b9..40ed4ec8 100644 --- a/src/gstools/covmodel/tpl_models.py +++ b/src/gstools/covmodel/tpl_models.py @@ -31,6 +31,19 @@ class TPLCovModel(CovModel): """Truncated-Power-Law Covariance Model base class for super-position.""" + @property + def intensity(self): + """:class:`float`: Intensity of variation.""" + return self.var / self.intensity_scaling + + @property + def intensity_scaling(self): + """:class:`float`: Scaling of Intensity to result in variance.""" + return ( + self.len_up_rescaled ** (2 * self.hurst) + - self.len_low_rescaled ** (2 * self.hurst) + ) / (2 * self.hurst) + @property def len_up(self): """:class:`float`: Upper length scale truncation of the model. @@ -55,13 +68,6 @@ def len_low_rescaled(self): """ return self.len_low / self.rescale - def var_factor(self): - """Factor for C (intensity of variation) to result in variance.""" - return ( - self.len_up_rescaled ** (2 * self.hurst) - - self.len_low_rescaled ** (2 * self.hurst) - ) / (2 * self.hurst) - def cor(self, h): """TPL - normalized correlation function.""" @@ -125,7 +131,7 @@ class TPLGaussian(TPLCovModel): * :math:`C>0` : scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. - You can access C directly by ``model.var_raw`` + You can access C by ``model.intensity`` * :math:`00` : scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. - You can access C directly by ``model.var_raw`` + You can access C by ``model.intensity`` * :math:`00` : scaling factor from the Power-Law (intensity of variation) This parameter will be calculated internally by the given variance. - You can access C directly by ``model.var_raw`` + You can access C by ``model.intensity`` * :math:`0