diff --git a/MANIFEST.in b/MANIFEST.in index 24184482a..7d3d0bb9e 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,4 +1,4 @@ prune ** -recursive-include tests *.py +recursive-include tests *.py *.txt recursive-include src/gstools *.py *.pyx include AUTHORS.md LICENSE README.md pyproject.toml setup.py diff --git a/docs/source/conf.py b/docs/source/conf.py index 3af3c1eeb..aaddc29f3 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -199,7 +199,7 @@ def setup(app): # -- Options for LaTeX output --------------------------------------------- # latex_engine = 'lualatex' -# logo to big +# logo too big latex_logo = "pics/gstools_150.png" # latex_show_urls = 'footnote' @@ -301,6 +301,7 @@ def setup(app): "../../examples/09_spatio_temporal/", "../../examples/10_normalizer/", "../../examples/11_plurigaussian/", + "../../examples/12_sum_model/", ], # path where to save gallery generated examples "gallery_dirs": [ @@ -316,6 +317,7 @@ def setup(app): "examples/09_spatio_temporal/", "examples/10_normalizer/", "examples/11_plurigaussian/", + "examples/12_sum_model/", ], # Pattern to search for example files "filename_pattern": r"\.py", diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index 4c04d6f02..e3b20db1b 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -23,6 +23,7 @@ explore its whole beauty and power. examples/09_spatio_temporal/index examples/10_normalizer/index examples/11_plurigaussian/index + examples/12_sum_model/index examples/00_misc/index .. only:: html diff --git a/examples/12_sum_model/00_simple_sum_model.py b/examples/12_sum_model/00_simple_sum_model.py new file mode 100644 index 000000000..a9ded98ed --- /dev/null +++ b/examples/12_sum_model/00_simple_sum_model.py @@ -0,0 +1,61 @@ +r""" +Creating a Sum Model +-------------------- + +This tutorial demonstrates how to create and use sum models in GSTools. +We'll combine a Spherical and a Gaussian covariance model to construct +a sum model, visualize its variogram, and generate spatial random fields. + +Let's start with importing GSTools and setting up the domain size. +""" + +import gstools as gs + +x = y = range(100) + +############################################################################### +# First, we create two individual covariance models: a :any:`Spherical` model and a +# :any:`Gaussian` model. The Spherical model with its short length scale +# will emphasize small-scale variability, while the Gaussian model with a larger length scale +# captures larger-scale patterns. + +m0 = gs.Spherical(dim=2, var=2.0, len_scale=5.0) +m1 = gs.Gaussian(dim=2, var=1.0, len_scale=10.0) + +############################################################################### +# Next, we create a sum model by adding these two models together. +# Let's visualize the resulting variogram alongside the individual models. + +model = m0 + m1 +ax = model.plot(x_max=20) +m0.plot(x_max=20, ax=ax) +m1.plot(x_max=20, ax=ax) + +############################################################################### +# As shown, the Spherical model controls the behavior at shorter distances, +# while the Gaussian model dominates at longer distances. The ratio of influence +# is thereby controlled by the provided variances of the individual models. +# +# Using the sum model, we can generate a spatial random field. Let's visualize +# the field created by the sum model. + +srf = gs.SRF(model, seed=20250107) +srf.structured((x, y)) +srf.plot() + +############################################################################### +# For comparison, we generate random fields using the individual models +# to observe their contributions more clearly. + +srf0 = gs.SRF(m0, seed=20250107) +srf0.structured((x, y)) +srf0.plot() + +srf1 = gs.SRF(m1, seed=20250107) +srf1.structured((x, y)) +srf1.plot() + +############################################################################### +# As seen, the Gaussian model introduces large-scale structures, while the +# Spherical model influences the field's roughness. The sum model combines +# these effects, resulting in a field that reflects multi-scale variability. diff --git a/examples/12_sum_model/01_fitting_sum_model.py b/examples/12_sum_model/01_fitting_sum_model.py new file mode 100644 index 000000000..7a0dfa294 --- /dev/null +++ b/examples/12_sum_model/01_fitting_sum_model.py @@ -0,0 +1,80 @@ +r""" +Fitting a Sum Model +-------------------- + +In this tutorial, we demonstrate how to fit a sum model consisting of two +covariance models to an empirical variogram. + +We will generate synthetic data, compute an empirical variogram, and fit a +sum model combining a Spherical and Gaussian model to it. +""" + +import gstools as gs + +x = y = range(100) + +############################################################################### +# First, we create a synthetic random field based on a known sum model. +# This will serve as the ground truth for fitting. + +# Define the true sum model +m0 = gs.Spherical(dim=2, var=2.0, len_scale=5.0) +m1 = gs.Gaussian(dim=2, var=1.0, len_scale=10.0) +true_model = m0 + m1 + +# Generate synthetic field +srf = gs.SRF(true_model, seed=20250405) +field = srf.structured((x, y)) + +############################################################################### +# Next, we calculate the empirical variogram from the synthetic data. + +# Compute empirical variogram +bin_center, gamma = gs.vario_estimate((x, y), field) + +############################################################################### +# Now we define a sum model to fit to the empirical variogram. +# Initially, the parameters of the models are arbitrary. +# +# A sum model can also be created by a list of model classes together with +# the common arguments (like dim in this case). + +fit_model = gs.SumModel(gs.Spherical, gs.Gaussian, dim=2) + +############################################################################### +# We fit the sum model to the empirical variogram using GSTools' built-in +# fitting capabilities. We deactivate the nugget fitting to not overparameterize +# our model. + +fit_model.fit_variogram(bin_center, gamma, nugget=False) +print(f"{true_model=}") +print(f" {fit_model=}") + +############################################################################### +# The variance of a sum model is the sum of the sub variances (:any:`SumModel.vars`) +# from the contained models. The length scale is a weighted sum of the sub +# length scales (:any:`SumModel.len_scales`) where the weights are the ratios +# of the sub variances to the total variance of the sum model. + +print(f"{true_model.var=:.2}, {true_model.len_scale=:.2}") +print(f" {fit_model.var=:.2}, {fit_model.len_scale=:.2}") + +############################################################################### +# After fitting, we can visualize the empirical variogram alongside the +# fitted sum model and its components. A Sum Model is subscriptable to access +# the individual models its contains. + +ax = fit_model.plot(x_max=max(bin_center)) +ax.scatter(bin_center, gamma) +# Extract individual components +fit_model[0].plot(x_max=max(bin_center), ax=ax) +fit_model[1].plot(x_max=max(bin_center), ax=ax) +# True models +true_model.plot(x_max=max(bin_center), ax=ax, ls="--", c="C0", label="") +true_model[0].plot(x_max=max(bin_center), ax=ax, ls="--", c="C1", label="") +true_model[1].plot(x_max=max(bin_center), ax=ax, ls="--", c="C2", label="") + +############################################################################### +# As we can see, the fitted sum model closely matches the empirical variogram, +# demonstrating its ability to capture multi-scale variability effectively. +# The "true" variograms are shown with dashed lines for comparison. diff --git a/examples/12_sum_model/README.rst b/examples/12_sum_model/README.rst new file mode 100644 index 000000000..c4ccb2306 --- /dev/null +++ b/examples/12_sum_model/README.rst @@ -0,0 +1,30 @@ +Summing Covariance Models +========================= + +In geostatistics, the spatial relations of natural phenomena is often represented using covariance models, +which describe how values of a property correlate over distance. +A single covariance model may capture specific features of the spatial correlation, such as smoothness or the range of influence. +However, many real-world spatial processes are complex, involving multiple overlapping structures +that cannot be adequately described by a single covariance model. + +This is where **sum models** come into play. +A sum model combines multiple covariance models into a single representation, +allowing for a more flexible and comprehensive description of spatial variability. +By summing covariance models, we can: + +1. **Capture Multi-Scale Variability:** Many spatial phenomena exhibit variability at different scales. + For example, soil properties may have small-scale variation due to local heterogeneities and large-scale variation due to regional trends. + A sum model can combine short-range and long-range covariance models to reflect this behavior. +2. **Improve Model Fit and Prediction Accuracy:** By combining models, sum models can better match empirical variograms or other observed data, + leading to more accurate predictions in kriging or simulation tasks. +3. **Enhance Interpretability:** Each component of a sum model can be associated with a specific spatial process or scale, + providing insights into the underlying mechanisms driving spatial variability. + +The new :any:`SumModel` introduced in GSTools makes it straightforward to define and work with such composite covariance structures. +It allows users to combine any number of base models, each with its own parameters, in a way that is both intuitive and computationally efficient. + +In the following tutorials, we'll explore how to use the :any:`SumModel` in GSTools, +including practical examples that demonstrate its utility in different scenarios. + +Examples +-------- diff --git a/pyproject.toml b/pyproject.toml index 4e98c50d2..4a5aa5c0d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -157,7 +157,7 @@ target-version = [ max-locals = 50 max-branches = 30 max-statements = 85 - max-attributes = 25 + max-attributes = 30 max-public-methods = 80 max-positional-arguments=20 diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index f89e8e4da..c13f3e5a4 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -55,6 +55,7 @@ .. autosummary:: CovModel + SumModel Covariance Models ^^^^^^^^^^^^^^^^^ @@ -63,6 +64,7 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~ .. autosummary:: + Nugget Gaussian Exponential Matern @@ -154,9 +156,11 @@ JBessel, Linear, Matern, + Nugget, Rational, Spherical, Stable, + SumModel, SuperSpherical, TPLExponential, TPLGaussian, @@ -199,6 +203,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 28ab81f21..76920c61e 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 23e198812..53a925f09 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -7,9 +7,10 @@ .. autosummary:: CovModel + SumModel """ -# pylint: disable=C0103, R0201, E1101, C0302, W0613 +# pylint: disable=C0103, R0201, E1101, C0302, W0613, W0231 import copy import numpy as np @@ -18,10 +19,20 @@ from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram +from gstools.covmodel.sum_tools import ( + default_mod_kwargs, + sum_check, + sum_compare, + sum_default_arg_bounds, + sum_default_opt_arg_bounds, + sum_model_repr, + sum_set_len_weights, + sum_set_var_weights, +) 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 +138,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 +149,9 @@ class CovModel: If present, they are described in the section `Other Parameters`. """ + _add_doc = True + """Flag to append the above doc-string about parameters to the model implementation.""" + def __init__( self, dim=3, @@ -159,7 +167,6 @@ def __init__( geo_scale=RADIAN_SCALE, temporal=False, spatial_dim=None, - var_raw=None, hankel_kw=None, **opt_arg, ): @@ -169,6 +176,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 +214,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 +227,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 +247,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 +478,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 +486,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. @@ -499,6 +498,17 @@ def percentile_scale(self, per=0.9): # spectrum methods (can be overridden for speedup) + @property + def needs_fourier_transform(self): + """ + Flag whether the model needs a fourier transform to calculate + the spectral density from the covariance function or if + it implements an analytical spectral density. + """ + base_method = getattr(CovModel, "spectral_density") + instance_method = getattr(self.__class__, "spectral_density") + return base_method == instance_method + def spectrum(self, k): r""" Spectrum of the covariance model. @@ -765,8 +775,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 +819,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 +835,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 +851,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 +867,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 +943,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 +970,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 +986,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 +1021,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): @@ -1069,12 +1044,13 @@ def hankel_kw(self): @hankel_kw.setter def hankel_kw(self, hankel_kw): - if self._hankel_kw is None or hankel_kw is None: - self._hankel_kw = copy.copy(HANKEL_DEFAULT) - if hankel_kw is not None: - self._hankel_kw.update(hankel_kw) - if self.dim is not None: - self._sft = SFT(ndim=self.dim, **self.hankel_kw) + if self.needs_fourier_transform: + if self._hankel_kw is None or hankel_kw is None: + self._hankel_kw = copy.copy(HANKEL_DEFAULT) + if hankel_kw is not None: + self._hankel_kw.update(hankel_kw) + if self.dim is not None: + self._sft = SFT(ndim=self.dim, **self.hankel_kw) @property def dist_func(self): @@ -1115,7 +1091,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 +1104,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 +1132,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 @@ -1196,15 +1172,597 @@ def __eq__(self, other): """Compare CovModels.""" if not isinstance(other, CovModel): return False + if isinstance(other, SumModel): + return False return compare(self, other) 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 or 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 have 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 += copy.deepcopy(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(copy.deepcopy(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 = None + 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 self.size + + 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_var_weights(self, weights, skip=None, var=None): + """ + Set variances of contained models by weights. + + The variances of the sub-models act as ratios for the sum-model. + This function enables to keep the total variance of the sum-model + and reset the individual variances of the contained sub-models + by the given weights. + + Parameters + ---------- + weights : iterable + Weights to set. Should have a length of len(models) - len(skip) + skip : iterable, optional + Model indices to skip. Should have compatible length, by default None + var : float, optional + Desired variance, by default current variance + + Raises + ------ + ValueError + If number of weights is not matching. + """ + sum_set_var_weights(self, weights, skip, var) + + def set_len_weights(self, weights, skip=None, len_scale=None): + """ + Set length scales of contained models by weights. + + The variances of the sub-models act as ratios for the sub-model length + scales to determine the total length scale of the sum-model. + This function enables to keep the total length scale of the sum-model as + well as the variances of the sub-models and reset the individual length scales + of the contained sub-models by the given weights. + + Parameters + ---------- + weights : iterable + Weights to set. Should have a length of len(models) - len(skip) + skip : iterable, optional + Model indices to skip. Should have compatible length, by default None + len_scale : float, optional + Desired len_scale, by default current len_scale + + Raises + ------ + ValueError + If number of weights is not matching. + """ + sum_set_len_weights(self, weights, skip, len_scale) + + @property + def models(self): + """:class:`tuple`: The summed models.""" + return self._models + + @property + def size(self): + """:class:`int`: Number of summed models.""" + return len(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, + ) + + @property + def needs_fourier_transform(self): + """Whether the model doesn't implement an analytical spectral density.""" + return False + + 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 [ + f"{o}_{i}" for o in ["var", "len_scale"] for i in range(self.size) + ] + + @property + def sub_arg_bounds(self): + """:class:`dict`: Names of the sub-arguments for var and len_scale.""" + return { + f"{o}_{i}": mod.arg_bounds[o] + for o in ["var", "len_scale"] + for (i, mod) in enumerate(self.models) + } + + @property + def arg_bounds(self): + """:class:`dict`: Bounds for all parameters. + + Notes + ----- + Keys are the arg names and values are lists of 2 or 3 values: + ``[a, b]`` or ``[a, b, ]`` where + is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` + to define if the bounds are open ("o") or closed ("c"). + """ + res = { + "var": self.var_bounds, + "len_scale": self.len_scale_bounds, + "nugget": self.nugget_bounds, + "anis": self.anis_bounds, + } + res.update(self.opt_arg_bounds) + res.update(self.sub_arg_bounds) + return res + + 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 __eq__(self, other): + """Compare SumModels.""" + if not isinstance(other, SumModel): + return False + return sum_compare(self, other) + + 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 8b19f4977..8ad510de4 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, @@ -166,32 +163,33 @@ def fit_variogram( The fitted parameters will be instantly set in the model. """ # preprocess selected parameters - para, sill, constrain_sill, anis = _pre_para( - model, para_select, sill, anis - ) + para, sill, anis, sum_cfg = _pre_para(model, para_select, sill, anis) # 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'") # prepare variogram data # => concatenate directional variograms to have a 1D array for x and y x_data, y_data, is_dir_vario = _check_vario(model, x_data, y_data) + # only fit anisotropy if a directional variogram was given + anis &= is_dir_vario + sub_fitting = sum_cfg.get("var_size", 0) + sum_cfg.get("len_size", 0) > 0 + if not (any(para.values()) or anis or sub_fitting): + raise ValueError("fit: no parameters selected for fitting.") # prepare init guess dictionary init_guess = _pre_init_guess( model, init_guess, np.mean(x_data), np.mean(y_data) ) - # only fit anisotropy if a directional variogram was given - anis &= is_dir_vario # set weights _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario) # set the lower/upper boundaries for the variogram-parameters bounds, init_guess_list = _init_curve_fit_para( - model, para, init_guess, constrain_sill, sill, anis + model, para, init_guess, sill, anis, sum_cfg ) # create the fitting curve curve_fit_kwargs["f"] = _get_curve( - model, para, constrain_sill, sill, anis, is_dir_vario + model, para, sill, anis, is_dir_vario, sum_cfg ) # set the remaining kwargs for curve_fit curve_fit_kwargs["bounds"] = bounds @@ -204,7 +202,9 @@ def fit_variogram( # fit the variogram popt, pcov = curve_fit(**curve_fit_kwargs) # convert the results - fit_para = _post_fitting(model, para, popt, anis, is_dir_vario) + fit_para = _post_fitting( + model, para, popt, sill, anis, is_dir_vario, sum_cfg + ) # calculate the r2 score if wanted if return_r2: return fit_para, pcov, _r2_score(model, x_data, y_data, is_dir_vario) @@ -213,70 +213,189 @@ def fit_variogram( def _pre_para(model, para_select, sill, anis): """Preprocess selected parameters.""" - var_last = False - var_tmp = 0.0 # init value - for par in para_select: - if par not in model.arg_bounds: + sum_cfg = {"fix": {}} + is_sum = hasattr(model, "sub_arg") + sub_args = getattr(model, "sub_arg", []) + valid_args = model.iso_arg + sub_args + bnd = model.arg_bounds + # if values given, set them in the model, afterwards all entries are bool + for par, val in para_select.items(): + if par not in valid_args: 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]) + # if parameters given with value, set it and deselect from fitting + if not isinstance(val, bool): + # don't set sub-args, var or len_scale in sum model yet + if is_sum and par in sub_args + ["var", "len_scale"]: + sum_cfg["fix"][par] = float(val) else: - setattr(model, par, float(para_select[par])) + setattr(model, par, float(val)) 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 sum models + if is_sum: + _check_sum(model, para_select, sum_cfg, bnd) # handling the sill sill = None if (isinstance(sill, bool) and sill) else sill if sill is not None: sill = model.sill if isinstance(sill, bool) else float(sill) - constrain_sill = True - sill_low = model.arg_bounds["var"][0] + model.arg_bounds["nugget"][0] - sill_up = model.arg_bounds["var"][1] + model.arg_bounds["nugget"][1] - if not sill_low <= sill <= sill_up: - raise ValueError("fit: sill out of bounds.") - if "var" in para_select and "nugget" in para_select: - if model.var > sill: - model.nugget = model.arg_bounds["nugget"][0] - model.var = sill - model.nugget - else: - model.nugget = sill - model.var - elif "var" in para_select: - if model.var > sill: - raise ValueError( - "fit: if sill is fixed and variance deselected, " - "the set variance should be less than the given sill." - ) - para_select["nugget"] = False - model.nugget = sill - model.var - elif "nugget" in para_select: - if model.nugget > sill: - raise ValueError( - "fit: if sill is fixed and nugget deselected, " - "the set nugget should be less than the given sill." - ) - para_select["var"] = False - model.var = sill - model.nugget - else: - # deselect the nugget, to recalculate it accordingly - # nugget = sill - var - para_select["nugget"] = False - 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}) + _check_sill(model, para_select, sill, bnd, sum_cfg) + # select all parameters to be fitted (if bounds do not indicate fixed parameter) + para = {par: bnd[par][0] < bnd[par][1] for par in valid_args} # now deselect unwanted parameters para.update(para_select) # check if anisotropy should be fitted or set if not isinstance(anis, bool): model.anis = anis anis = False - return para, sill, constrain_sill, anis + return para, sill, anis, sum_cfg + + +def _check_sill(model, para_select, sill, bnd, sum_cfg): + """ + This functions checks if the selected values for + variance, nugget and sill are valid. + """ + is_sum = hasattr(model, "sub_arg") + sill_low = bnd["var"][0] + bnd["nugget"][0] + sill_up = bnd["var"][1] + bnd["nugget"][1] + if not sill_low <= sill <= sill_up: + raise ValueError("fit: sill out of bounds.") + if is_sum: + var_fixed = sum_cfg["var_fix"] + else: + var_fixed = "var" in para_select + if var_fixed and "nugget" in para_select: + if not np.isclose(model.var + model.nugget, sill): + msg = "fit: if sill, var and nugget are fixed, var + nugget should match the given sill" + raise ValueError(msg) + elif var_fixed: + if model.var > sill: + raise ValueError( + "fit: if sill is fixed and variance deselected, " + "the set variance should be less than the given sill." + ) + para_select["nugget"] = False + # this also works for a pure nugget model + model.nugget = sill - model.var + elif "nugget" in para_select: + if model.nugget > sill: + raise ValueError( + "fit: if sill is fixed and nugget deselected, " + "the set nugget should be less than the given sill." + ) + para_select["var"] = False + var = sill - model.nugget + if is_sum: + model.set_var_weights( + np.ones_like(sum_cfg["var_fit"]), sum_cfg["var_skip"], var + ) + # if only 1 sub-var was to fit, this is now also fixed + if len(sum_cfg["var_fit"]) == 1: + i = sum_cfg["var_fit"][0] + sum_cfg["fix"].setdefault(f"var_{i}", model.vars[i]) + sum_cfg["var_fit"] = [] + sum_cfg["var_skip"] = list(range(model.size)) + sum_cfg["fix"].setdefault("var", var) + sum_cfg["var_fix"] = True + # number or sub-var parameters + sum_cfg["var_size"] = max( + model.size - len(sum_cfg["var_skip"]) - 1, 0 + ) + else: + # in case of a nugget model, this should raise an error + model.var = var + else: + # deselect the nugget, to recalculate it accordingly + # nugget = sill - var + para_select["nugget"] = False + + +def _check_sum(model, para_select, sum_cfg, bnd): + """Check for consistent parameter selection in case of a SumModel.""" + # check len_scale + if "len_scale" in para_select: + for par in para_select: + if par.startswith("len_scale_"): + msg = ( + "fit: for sum-models you can only fix " + "'len_scale' or the sub-arguments 'len_scale_', not both." + ) + raise ValueError(msg) + sum_cfg["fix"].setdefault("len_scale", model.len_scale) + # use len_scale_ for fitting if len_scale not fixed + # use weights for fitting if len_scale fixed + # either way: len_scale not used for fitting in sum_model + para_select["len_scale"] = False + # check variance + if "var" in para_select: + sum_cfg["fix"].setdefault("var", model.var) + # use var_ for fitting if var not fixed + # use weights for fitting if var fixed + # either way: var not used for fitting in sum_model + para_select["var"] = False + # whether var and len_scale are fixed + var_fix = "var" in sum_cfg["fix"] + len_fix = "len_scale" in sum_cfg["fix"] + size = model.size + # check for fixed bounds + for i in range(size): + for par in [f"var_{i}", f"len_scale_{i}"]: + if not bnd[par][0] < bnd[par][1]: + para_select[par] = False + # check sub arguments (var_ and len_scale_) + var_skip = [] + len_skip = [] + remove = [] + for par in para_select: + if par.startswith("var_"): + var_skip.append(int(par[4:])) + # for fixed var, fit by weights + if var_fix: + remove.append(par) + if par.startswith("len_scale_"): + len_skip.append(int(par[10:])) + for par in remove: + para_select.pop(par) + var_skip.sort() + len_skip.sort() + var_fit = sorted(set(range(size)) - set(var_skip)) + len_fit = sorted(set(range(size)) - set(len_skip)) + # if all sub-vars fixed, total variance is fixed + if not var_fit: + para_select["var"] = False + # set values related to var and len_scale in sum-model + for i in var_skip: + sum_cfg["fix"].setdefault(f"var_{i}", model.vars[i]) + setattr(model, f"var_{i}", sum_cfg["fix"][f"var_{i}"]) + for i in len_skip: + sum_cfg["fix"].setdefault(f"len_scale_{i}", model.len_scales[i]) + setattr(model, f"len_scale_{i}", sum_cfg["fix"][f"len_scale_{i}"]) + if var_fix: + var_min = sum(sum_cfg["fix"][f"var_{i}"] for i in var_skip) + if var_min > sum_cfg["fix"]["var"]: + msg = "fit: fixed sub-variances greater than fixed total variance." + raise ValueError(msg) + model.set_var_weights( + np.ones_like(var_fit), var_skip, sum_cfg["fix"]["var"] + ) + # if all sub_vars except one is fixed and total var as well, all are fixed + if len(var_fit) == 1: + var_fit = [] + var_skip = list(range(size)) + if len_fix: + model.len_scale = sum_cfg["fix"]["len_scale"] + # update config for the sum-model fitting + sum_cfg["var_skip"] = var_skip # sub vars to skip + sum_cfg["len_skip"] = len_skip # sub lens to skip + sum_cfg["var_fit"] = var_fit # sub vars to fit + sum_cfg["len_fit"] = len_fit # sub lens to fit + sum_cfg["var_fix"] = var_fix # total variance fixed + sum_cfg["len_fix"] = len_fix # total len-scale fixed + # number or sub-var parameters + sum_cfg["var_size"] = max(size - len(var_skip) - int(var_fix), 0) + # number or sub-len parameters + sum_cfg["len_size"] = max(size - len(len_skip) - int(len_fix), 0) def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0): @@ -288,15 +407,17 @@ def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0): if default_guess not in ["default", "current"]: raise ValueError(f"fit_variogram: unknown def. guess: {default_guess}") default = default_guess == "default" + is_sum = hasattr(model, "sub_arg") + sub_args = getattr(model, "sub_arg", []) + valid_args = model.iso_arg + sub_args + ["anis"] # check invalid names for given init guesses - invalid_para = set(init_guess) - set(model.iso_arg + ["anis"]) + invalid_para = set(init_guess) - set(valid_args) if invalid_para: raise ValueError(f"fit_variogram: unknown init guess: {invalid_para}") bnd = model.arg_bounds # default length scale is mean of given bin centers (respecting "rescale") - init_guess.setdefault( - "len_scale", mean_x * model.rescale if default else model.len_scale - ) + len_def = mean_x * model.rescale + init_guess.setdefault("len_scale", len_def if default else model.len_scale) # init guess for variance and nugget is mean of given variogram for par in ["var", "nugget"]: init_guess.setdefault(par, mean_y if default else getattr(model, par)) @@ -316,8 +437,17 @@ def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0): else getattr(model, opt) ), ) + # SumModel: check for var_ and len_scale_ and set defaults + if is_sum: + for i in range(model.size): + init_guess.setdefault( + f"len_scale_{i}", len_def if default else model.len_scales[i] + ) + init_guess.setdefault( + f"var_{i}", mean_y / model.size if default else model.vars[i] + ) # convert all init guesses to float (except "anis") - for arg in model.iso_arg: + for arg in model.iso_arg + sub_args: init_guess[arg] = float(init_guess[arg]) return init_guess @@ -363,34 +493,59 @@ def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario): curve_fit_kwargs["absolute_sigma"] = True -def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): +def _init_curve_fit_para(model, para, init_guess, sill, anis, sum_cfg): """Create initial guess and bounds for fitting.""" + is_sum = hasattr(model, "sub_arg") low_bounds = [] top_bounds = [] init_guess_list = [] - for par in DEFAULT_PARA: + bnd = model.arg_bounds + 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 + low_bounds.append(bnd[par][0]) + if par == "var" and sill is not None: # var <= sill in this case top_bounds.append(sill) else: - top_bounds.append(model.arg_bounds[par][1]) + top_bounds.append(bnd[par][1]) init_guess_list.append( _init_guess( bounds=[low_bounds[-1], top_bounds[-1]], 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 is_sum: + if sum_cfg["var_fix"]: + for _ in range(sum_cfg["var_size"]): + low_bounds.append(0.0) + top_bounds.append(1.0) + init_guess_list.append(0.5) + else: + for i in sum_cfg["var_fit"]: + par = f"var_{i}" + low_bounds.append(bnd[par][0]) + top_bounds.append(bnd[par][1]) + init_guess_list.append( + _init_guess( + bounds=[low_bounds[-1], top_bounds[-1]], + default=init_guess[par], + ) + ) + if sum_cfg["len_fix"]: + for _ in range(sum_cfg["len_size"]): + low_bounds.append(0.0) + top_bounds.append(1.0) + init_guess_list.append(0.5) + else: + for i in sum_cfg["len_fit"]: + par = f"len_scale_{i}" + low_bounds.append(bnd[par][0]) + top_bounds.append(bnd[par][1]) + init_guess_list.append( + _init_guess( + bounds=[low_bounds[-1], top_bounds[-1]], + default=init_guess[par], + ) ) - ) if anis: for i in range(model.dim - 1): low_bounds.append(model.anis_bounds[0]) @@ -411,42 +566,53 @@ def _init_guess(bounds, default): return default_arg_from_bounds(bounds) -def _get_curve(model, para, constrain_sill, sill, anis, is_dir_vario): +def _get_curve(model, para, sill, anis, is_dir_vario, sum_cfg): """Create the curve for scipys curve_fit.""" - var_save = model.var + is_sum = hasattr(model, "sub_arg") # 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 + # set var and len-scale ratios in sum-models + if is_sum: + if sum_cfg["var_size"] > 0: + var_vals = args[para_skip : para_skip + sum_cfg["var_size"]] + para_skip += sum_cfg["var_size"] + if sum_cfg["var_fix"]: + model.set_var_weights( + stick_breaking_uniform(var_vals), + sum_cfg["var_skip"], + sum_cfg["fix"]["var"], + ) + else: + for i, val in zip(sum_cfg["var_fit"], var_vals): + setattr(model, f"var_{i}", val) + if sum_cfg["len_size"] > 0: + len_vals = args[para_skip : para_skip + sum_cfg["len_size"]] + para_skip += sum_cfg["len_size"] + if sum_cfg["len_fix"]: + model.set_len_weights( + stick_breaking_uniform(len_vals), + sum_cfg["len_skip"], + sum_cfg["fix"]["len_scale"], + ) + else: + for i, val in zip(sum_cfg["len_fit"], len_vals): + setattr(model, f"len_scale_{i}", val) + # handle sill + if sill is not None and para["var"]: + 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 :] @@ -460,36 +626,56 @@ def curve(x, arg1, *args): return curve -def _post_fitting(model, para, popt, anis, is_dir_vario): +def _post_fitting(model, para, popt, sill, anis, is_dir_vario, sum_cfg): """Postprocess fitting results and application to model.""" + is_sum = hasattr(model, "sub_arg") 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) + # set var and len-scale ratios in sum-models + if is_sum: + if sum_cfg["var_size"] > 0: + var_vals = popt[para_skip : para_skip + sum_cfg["var_size"]] + para_skip += sum_cfg["var_size"] + if sum_cfg["var_fix"]: + model.set_var_weights( + stick_breaking_uniform(var_vals), + sum_cfg["var_skip"], + sum_cfg["fix"]["var"], + ) + else: + for i, val in zip(sum_cfg["var_fit"], var_vals): + setattr(model, f"var_{i}", val) + if sum_cfg["len_size"] > 0: + len_vals = popt[para_skip : para_skip + sum_cfg["len_size"]] + para_skip += sum_cfg["len_size"] + if sum_cfg["len_fix"]: + model.set_len_weights( + stick_breaking_uniform(len_vals), + sum_cfg["len_skip"], + sum_cfg["fix"]["len_scale"], + ) + else: + for i, val in zip(sum_cfg["len_fit"], len_vals): + setattr(model, f"len_scale_{i}", val) + for i in range(model.size): + fit_para[f"var_{i}"] = model.vars[i] + fit_para[f"len_scale_{i}"] = model.len_scales[i] + # handle sill + if sill is not None and para["var"]: + nugget = sill - model.var + fit_para["nugget"] = nugget + model.nugget = nugget 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 @@ -538,3 +724,34 @@ def func(x_data): return 1.0 / (1.0 + np.exp(growth * (x_mean - x_data))) return func + + +def stick_breaking_uniform(u): + """ + Generate a single sample (x_1, ..., x_n) uniformly from the (n-1)-simplex. + + This is using Beta transforms of uniform samples. The x_i will sum to 1. + + Parameters + ---------- + u : array-like of shape (n-1,) + Uniform(0,1) random values + + Returns + ------- + x : ndarray of shape (n,) + A random sample in the (n-1)-simplex. + """ + n = len(u) + 1 + x = np.zeros(n) + leftover = 1.0 + for i in range(n - 1): + # 2) Compute the inverse CDF of Beta(1, b) = Beta(1, n-1-i) + fraction = 1.0 - (1.0 - u[i]) ** (1.0 / (n - 1 - i)) + # 3) Break that fraction of the current leftover + x[i] = leftover * fraction + # 4) Subtract that from leftover + leftover -= x[i] + # Last coordinate + x[-1] = leftover + return x diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index b1a9d68ec..cb495751f 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,72 @@ ] +class Nugget(SumModel): + r"""Pure nugget model. + + This model has no correlated variability and represents pure noise. + The length scale of the model will be zero. + + 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 32148c14a..f7fa58db5 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 000000000..253ee5a4b --- /dev/null +++ b/src/gstools/covmodel/sum_tools.py @@ -0,0 +1,266 @@ +""" +GStools subpackage providing tools for sum-models. + +.. currentmodule:: gstools.covmodel.sum_tools + +The following classes and functions are provided + +.. autosummary:: + ARG_DEF + default_mod_kwargs + sum_check + sum_compare + 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.tools import RADIAN_SCALE +from gstools.tools.misc import list_format + +__all__ = [ + "ARG_DEF", + "default_mod_kwargs", + "sum_check", + "sum_compare", + "sum_default_arg_bounds", + "sum_default_opt_arg_bounds", + "sum_set_var_weights", + "sum_set_len_weights", + "sum_model_repr", +] + + +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_var_weights(summod, weights, skip=None, var=None): + """ + Set variances of contained models by weights. + + Parameters + ---------- + weights : iterable + Weights to set. Should have a length of len(models) - len(skip) + skip : iterable, optional + Model indices to skip. Should have compatible length, by default None + var : float, optional + Desired variance, by default current variance + + Raises + ------ + ValueError + If number of weights is not matching. + """ + skip = set() if skip is None else set(skip) + if len(summod) != len(weights) + len(skip): + msg = "SumModel.set_var_weights: number of ratios not matching." + raise ValueError(msg) + ids = range(len(summod)) + if fail := set(skip) - set(ids): + msg = ( + f"SumModel.set_var_weights: ids given by 'skip' not valid: {fail}" + ) + raise ValueError(msg) + var = summod.var if var is None else float(var) + var_sum = sum(summod.models[i].var for i in skip) + var_diff = var - var_sum + if var_diff < 0: + msg = ( + "SumModel.set_var_weights: summed variances selected " + "with 'skip' already too big to keep total variance." + ) + raise ValueError(msg) + weights_sum = sum(weights) + var_list = summod.vars + j = 0 + for i in ids: + if i in skip: + continue + var_list[i] = var_diff * weights[j] / weights_sum + j += 1 + summod.vars = var_list + + +def sum_set_len_weights(summod, weights, skip=None, len_scale=None): + """ + Set length scales of contained models by weights. + + Parameters + ---------- + weights : iterable + Weights to set. Should have a length of len(models) - len(skip) + skip : iterable, optional + Model indices to skip. Should have compatible length, by default None + len_scale : float, optional + Desired len_scale, by default current len_scale + + Raises + ------ + ValueError + If number of weights is not matching. + """ + skip = set() if skip is None else set(skip) + if len(summod) != len(weights) + len(skip): + msg = "SumModel.set_len_weights: number of weights not matching." + raise ValueError(msg) + ids = range(len(summod)) + if fail := set(skip) - set(ids): + msg = ( + f"SumModel.set_len_weights: ids given by 'skip' not valid: {fail}" + ) + raise ValueError(msg) + len_scale = summod.len_scale if len_scale is None else float(len_scale) + # also skip models with no variance (not contributing to total len scale) + j = 0 + wei = [] + for i in ids: + if i in skip: + continue + if np.isclose(summod.ratios[i], 0): + skip.add(i) + else: + wei.append(weights[j]) + j += 1 + weights = wei + len_sum = sum(summod[i].len_scale * summod.ratios[i] for i in skip) + len_diff = len_scale - len_sum + if len_diff < 0: + msg = ( + "SumModel.set_len_weights: summed length scales " + "selected with 'skip' already too big to keep total length scale." + ) + raise ValueError(msg) + weights_sum = sum(weights) + len_scales = summod.len_scales + j = 0 + for i in ids: + if i in skip: + continue + len_scales[i] = len_diff * weights[j] / weights_sum / summod.ratios[j] + j += 1 + summod.len_scales = len_scales + + +def sum_compare(this, that): + """ + Compare SumModels. + + Parameters + ---------- + this / that : :any:`SumModel` + The sum models to compare. + """ + if len(this) != len(that): + return False + if not np.isclose(this.nugget, that.nugget): + return False + return all(mod1 == mod2 for (mod1, mod2) in zip(this, that)) + + +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 dddeb4413..cfba8259a 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,18 @@ 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}" - ) - if arg in model.opt_arg: + msg = f"Given bounds for '{arg}' are not valid, got: {bounds}" + raise ValueError(msg) + if arg in getattr(model, "sub_arg", []): + # var_ and len_scale_ + name, i = _split_sub(arg) + setattr(model[i], f"_{name}_bounds", bounds) + elif arg in model.opt_arg: model._opt_arg_bounds[arg] = bounds - elif arg == "var": - var_bnds = bounds - continue - elif arg == "len_scale": - model.len_scale_bounds = bounds - elif arg == "nugget": - model.nugget_bounds = bounds - elif arg == "anis": - model.anis_bounds = bounds + elif arg in model.arg_bounds: + setattr(model, f"_{arg}_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 +485,15 @@ 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 _split_sub(name): + if name.startswith("var_"): + return "var", int(name[4:]) + if name.startswith("len_scale_"): + return "len_scale", int(name[10:]) + msg = f"Unknown sub variable: {name}" + raise ValueError(msg) def check_arg_bounds(model): @@ -501,19 +512,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): @@ -556,17 +555,17 @@ 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 + if model.needs_fourier_transform: + model._sft = SFT(ndim=model.dim, **model.hankel_kw) + # 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 +586,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 +608,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 b728e7b98..40ed4ec8e 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