From 94b6c2988a32da57246b2adc4b264bec27573fb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 26 Jul 2024 23:11:40 +0200 Subject: [PATCH 01/28] plot: handle len_scale of 0 --- src/gstools/covmodel/plot.py | 75 ++++++++++++++++++++++++++++-------- 1 file changed, 60 insertions(+), 15 deletions(-) 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) From abbb062cffff08391bf64e6c3ed96ffee5996a44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 26 Jul 2024 23:12:34 +0200 Subject: [PATCH 02/28] covmodel: init len_scale with integral_scale if given; remove unused _integral_scale --- src/gstools/covmodel/base.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 23e19881..f5def43c 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -208,6 +208,8 @@ def __init__( 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 ) @@ -221,7 +223,6 @@ def __init__( 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: @@ -486,8 +487,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. @@ -1043,8 +1043,7 @@ 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): @@ -1054,7 +1053,7 @@ def integral_scale(self, integral_scale): integral_scale = self.len_scale # reset len_scale self.len_scale = 1.0 - int_tmp = self.calc_integral_scale() + int_tmp = self.integral_scale self.len_scale = integral_scale / int_tmp if not np.isclose(self.integral_scale, integral_scale, rtol=1e-3): raise ValueError( From c39866a90096a71db11092b746ca979e78c991a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 26 Jul 2024 23:13:14 +0200 Subject: [PATCH 03/28] bounds: allow interval with equal bounds if closed --- src/gstools/covmodel/tools.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index dddeb441..cf310b94 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -297,11 +297,14 @@ 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 From df9f2baa949201a6f54f41a2473f074909365f38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 27 Jul 2024 01:04:30 +0200 Subject: [PATCH 04/28] CovModel: add force_values logic --- src/gstools/covmodel/base.py | 65 ++++++++++++++++++++++++----------- src/gstools/covmodel/fit.py | 4 +++ src/gstools/covmodel/tools.py | 10 +++--- 3 files changed, 55 insertions(+), 24 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index f5def43c..3a61e229 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -139,6 +139,8 @@ class CovModel: used for the spectrum calculation. Use with caution (Better: Don't!). ``None`` is equivalent to ``{"a": -1, "b": 1, "N": 1000, "h": 0.001}``. Default: :any:`None` + fixed: :class:`set` or :any:`None`, optional + Names of fixed arguments. Default: :any:`None` **opt_arg Optional arguments are covered by these keyword arguments. If present, they are described in the section `Other Parameters`. @@ -161,6 +163,7 @@ def __init__( spatial_dim=None, var_raw=None, hankel_kw=None, + fixed=None, **opt_arg, ): # assert, that we use a subclass @@ -169,6 +172,16 @@ def __init__( if not hasattr(self, "variogram"): raise TypeError("Don't instantiate 'CovModel' directly!") + # indicate that arguments are fixed (True after __init__) + self._fix = False + # force input values + forced = self.force_values() + var = forced.get("var", var) + len_scale = forced.get("len_scale", len_scale) + nugget = forced.get("nugget", nugget) + integral_scale = forced.get("integral_scale", integral_scale) + for opt in opt_arg: + opt_arg[opt] = forced.get(opt, opt_arg[opt]) # prepare dim setting self._dim = None self._hankel_kw = None @@ -198,6 +211,16 @@ def __init__( # optional arguments for the variogram-model set_opt_args(self, opt_arg) + # check fixed + fixed = set(fixed) if fixed is not None else set() + fixed = fixed | self.always_fixed() + valid_fixed = set(self.iso_arg) # | set(["anis"]) + if not fixed <= valid_fixed: + raise ValueError( + f"CovModel: unknown names in 'fixed': {fixed - valid_fixed}" + ) + self._fixed = fixed + # set standard boundaries for variance, len_scale, nugget and opt_arg bounds = self.default_arg_bounds() bounds.update(self.default_opt_arg_bounds()) @@ -234,6 +257,11 @@ def __init__( self.check_arg_bounds() # additional checks for the optional arguments (provided by user) self.check_opt_arg() + # set fixed bounds after checking original bounds + for arg in self.fixed: + val = getattr(self, arg) + self.set_arg_bounds(check_args=False, **{arg: (val, val, "cc")}) + self._fix = True # precision for printing self._prec = 3 @@ -438,6 +466,14 @@ def pykrige_kwargs(self): # methods for optional/default arguments (can be overridden) + def always_fixed(self): + """Provide set of fixed arguments.""" + return set() + + def force_values(self): + """:class:`dict`: Forced values for arguments.""" + return {} + def default_opt_arg(self): """Provide default optional arguments by the user. @@ -795,6 +831,11 @@ def check_arg_bounds(self): # bounds properties + @property + def fixed(self): + """:class:`set`: Set with names of fixed arguments.""" + return self._fixed + @property def var_bounds(self): """:class:`list`: Bounds for the variance. @@ -809,11 +850,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 +866,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 +882,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 +898,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): diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 8b19f497..f19eaac3 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -165,10 +165,14 @@ def fit_variogram( The fitted parameters will be instantly set in the model. """ + for para in model.fixed: + para_select[para] = False # preprocess selected parameters 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 # check method diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index cf310b94..ad731aeb 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -458,6 +458,8 @@ def set_arg_bounds(model, check_args=True, **kwargs): # if variance needs to be resetted, do this at last var_bnds = [] for arg, bounds in kwargs.items(): + if model._fix and arg in model.fixed: + raise ValueError(f"Can't set bounds for fixed argument '{arg}'") if not check_bounds(bounds): raise ValueError( f"Given bounds for '{arg}' are not valid, got: {bounds}" @@ -468,11 +470,11 @@ def set_arg_bounds(model, check_args=True, **kwargs): var_bnds = bounds continue 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: @@ -483,7 +485,7 @@ def set_arg_bounds(model, check_args=True, **kwargs): setattr(model, arg, def_arg) # set var last like always if var_bnds: - model.var_bounds = 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) From a3f1ec6d897ba161f93abd133db7480ee70fa5d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sat, 27 Jul 2024 23:05:10 +0200 Subject: [PATCH 05/28] generator: shortcut for 0 variance --- src/gstools/field/generator.py | 75 ++++++++++++++++++++++------------ 1 file changed, 50 insertions(+), 25 deletions(-) diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index d4cb0dea..89914785 100755 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -152,6 +152,16 @@ def __call__(self, pos, add_nugget=True): the random modes """ + @property + @abstractmethod + def model(self): + """:any:`CovModel`: Covariance model of the spatial random field.""" + + @property + def zero_var(self): + """:class:`bool`: Whether Covariance model has zero variance.""" + return np.isclose(self.model.var, 0) + @property @abstractmethod def value_type(self): @@ -260,6 +270,10 @@ def __call__(self, pos, add_nugget=True): the random modes """ pos = np.asarray(pos, dtype=np.double) + if self.zero_var: + shp = pos.shape[1:] + return self.get_nugget(shp) if add_nugget else np.full(shp, 0.0) + # generate if var is not 0 summed_modes = _summate( self._cov_sample, self._z_1, self._z_2, pos, config.NUM_THREADS ) @@ -285,7 +299,7 @@ def get_nugget(self, shape): size=shape ) else: - nugget = 0.0 + nugget = np.full(shape, 0.0) if self.zero_var else 0.0 return nugget def update(self, model=None, seed=np.nan): @@ -362,23 +376,26 @@ def reset_seed(self, seed=np.nan): self._z_1 = self._rng.random.normal(size=self._mode_no) self._z_2 = self._rng.random.normal(size=self._mode_no) # sample uniform on a sphere - sphere_coord = self._rng.sample_sphere(self.model.dim, self._mode_no) - # sample radii according to radial spectral density of the model - if self.sampling == "inversion" or ( - self.sampling == "auto" and self.model.has_ppf - ): - pdf, cdf, ppf = self.model.dist_func - rad = self._rng.sample_dist( - size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0 - ) + if self.zero_var: + self._cov_sample = np.full((self.model.dim, self._mode_no), 0.0) else: - rad = self._rng.sample_ln_pdf( - ln_pdf=self.model.ln_spectral_rad_pdf, - size=self._mode_no, - sample_around=1.0 / self.model.len_rescaled, - ) - # get fully spatial samples by multiplying sphere samples and radii - self._cov_sample = rad * sphere_coord + sph_crd = self._rng.sample_sphere(self.model.dim, self._mode_no) + # sample radii according to radial spectral density of the model + if self.sampling == "inversion" or ( + self.sampling == "auto" and self.model.has_ppf + ): + pdf, cdf, ppf = self.model.dist_func + rad = self._rng.sample_dist( + size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0 + ) + else: + rad = self._rng.sample_ln_pdf( + ln_pdf=self.model.ln_spectral_rad_pdf, + size=self._mode_no, + sample_around=1.0 / self.model.len_rescaled, + ) + # get spatial samples by multiplying sphere samples and radii + self._cov_sample = rad * sph_crd @property def sampling(self): @@ -541,6 +558,10 @@ def __call__(self, pos, add_nugget=True): the random modes """ pos = np.asarray(pos, dtype=np.double) + nugget = self.get_nugget(pos.shape) if add_nugget else 0.0 + e1 = self._create_unit_vector(pos.shape) + if self.zero_var: + return self.mean_u * e1 + nugget summed_modes = _summate_incompr( self._cov_sample, self._z_1, @@ -548,8 +569,6 @@ def __call__(self, pos, add_nugget=True): pos, config.NUM_THREADS, ) - nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0 - e1 = self._create_unit_vector(summed_modes.shape) return ( self.mean_u * e1 + self.mean_u @@ -667,7 +686,10 @@ def __call__(self, pos, add_nugget=True): the random modes """ pos = np.asarray(pos, dtype=np.double) - + if self.zero_var: + shp = pos.shape[1:] + return self.get_nugget(shp) if add_nugget else np.full(shp, 0.0) + # generate if var is not 0 summed_modes = _summate_fourier( self._spectrum_factor, self._modes, @@ -698,7 +720,7 @@ def get_nugget(self, shape): size=shape ) else: - nugget = 0.0 + nugget = np.full(shape, 0.0) if self.zero_var else 0.0 return nugget def update(self, model=None, seed=np.nan, period=None, mode_no=None): @@ -800,10 +822,13 @@ def reset_seed(self, seed=np.nan): self._z_2 = self._rng.random.normal(size=np.prod(self._mode_no)) # pre calc. the spectrum for all wave numbers they are handed over to # Cython, which doesn't have access to the CovModel - k_norm = np.linalg.norm(self._modes, axis=0) - self._spectrum_factor = np.sqrt( - self._model.spectrum(k_norm) * np.prod(self._delta_k) - ) + if self.zero_var: + self._spectrum_factor = np.full(np.prod(self._mode_no), 0.0) + else: + k_norm = np.linalg.norm(self._modes, axis=0) + self._spectrum_factor = np.sqrt( + self._model.spectrum(k_norm) * np.prod(self._delta_k) + ) def _fill_to_dim( self, values, dim, dtype=float, default_value=None From 8b3ad819e3cf4431fa66f37bd11ace80d850db52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 29 Jul 2024 00:42:25 +0200 Subject: [PATCH 06/28] CovModel: safer usage of privat attr; anis set bug fix; simpler int scale --- src/gstools/covmodel/base.py | 51 +++++++++++++++++------------------- 1 file changed, 24 insertions(+), 27 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 3a61e229..dcbffe57 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -173,7 +173,7 @@ def __init__( raise TypeError("Don't instantiate 'CovModel' directly!") # indicate that arguments are fixed (True after __init__) - self._fix = False + self._init = False # force input values forced = self.force_values() var = forced.get("var", var) @@ -187,6 +187,7 @@ def __init__( self._hankel_kw = None self._sft = None # prepare parameters (they are checked in dim setting) + self._fixed = None self._rescale = None self._len_scale = None self._anis = None @@ -246,7 +247,8 @@ def __init__( self.var = var else: self._var = float(var_raw) - self.integral_scale = integral_scale + if integral_scale is not None: + self.integral_scale = integral_scale # set var again, if int_scale affects var_factor if var_raw is None: self._var = None @@ -261,7 +263,7 @@ def __init__( for arg in self.fixed: val = getattr(self, arg) self.set_arg_bounds(check_args=False, **{arg: (val, val, "cc")}) - self._fix = True + self._init = True # precision for printing self._prec = 3 @@ -801,8 +803,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"), } @@ -1014,10 +1016,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 @@ -1033,7 +1032,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): @@ -1072,19 +1071,17 @@ def integral_scale(self): @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.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, integral_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): @@ -1139,7 +1136,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): @@ -1152,7 +1149,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): @@ -1180,7 +1177,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 @@ -1226,7 +1223,7 @@ 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 __repr__(self): From a60b3d7e9e6066f986ea4baea60e0f78c7e1feb5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 29 Jul 2024 00:43:12 +0200 Subject: [PATCH 07/28] CovModel: check _init in bounds setter; also compare geo_scale; simpler repr --- src/gstools/covmodel/tools.py | 50 +++++++++++++---------------------- 1 file changed, 18 insertions(+), 32 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index ad731aeb..9d53de58 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -458,7 +458,7 @@ def set_arg_bounds(model, check_args=True, **kwargs): # if variance needs to be resetted, do this at last var_bnds = [] for arg, bounds in kwargs.items(): - if model._fix and arg in model.fixed: + if model._init and arg in model.fixed: raise ValueError(f"Can't set bounds for fixed argument '{arg}'") if not check_bounds(bounds): raise ValueError( @@ -598,6 +598,7 @@ def compare(this, that): 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: @@ -614,39 +615,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})" From 92210adabaae3c5052830d027d712bd71d1fce26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 29 Jul 2024 00:45:59 +0200 Subject: [PATCH 08/28] CovModel: remove force arguments mechanic --- src/gstools/covmodel/base.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index dcbffe57..d1d1cf0b 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -174,14 +174,6 @@ def __init__( # indicate that arguments are fixed (True after __init__) self._init = False - # force input values - forced = self.force_values() - var = forced.get("var", var) - len_scale = forced.get("len_scale", len_scale) - nugget = forced.get("nugget", nugget) - integral_scale = forced.get("integral_scale", integral_scale) - for opt in opt_arg: - opt_arg[opt] = forced.get(opt, opt_arg[opt]) # prepare dim setting self._dim = None self._hankel_kw = None @@ -472,10 +464,6 @@ def always_fixed(self): """Provide set of fixed arguments.""" return set() - def force_values(self): - """:class:`dict`: Forced values for arguments.""" - return {} - def default_opt_arg(self): """Provide default optional arguments by the user. From 8999bf98088446a91d728c4be376db474f7d06b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 29 Jul 2024 00:50:43 +0200 Subject: [PATCH 09/28] pylint fixes --- pyproject.toml | 2 +- src/gstools/covmodel/base.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) 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/covmodel/base.py b/src/gstools/covmodel/base.py index d1d1cf0b..56f776af 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -21,7 +21,6 @@ from gstools.covmodel.tools import ( _init_subclass, check_arg_bounds, - check_bounds, compare, default_arg_from_bounds, model_repr, From d643b450b86ff82f29f15c092c7d8fa271a1466d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 11 Aug 2024 12:16:13 +0200 Subject: [PATCH 10/28] TPL: remove var_raw and var_factor logic and only add intensity as property to TPL models --- src/gstools/covmodel/base.py | 44 ++++-------------------------- src/gstools/covmodel/fit.py | 32 ++++------------------ src/gstools/covmodel/tools.py | 1 - src/gstools/covmodel/tpl_models.py | 26 +++++++++++------- 4 files changed, 28 insertions(+), 75 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 56f776af..fce1e64b 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -126,12 +126,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` @@ -160,7 +154,6 @@ def __init__( geo_scale=RADIAN_SCALE, temporal=False, spatial_dim=None, - var_raw=None, hankel_kw=None, fixed=None, **opt_arg, @@ -232,20 +225,12 @@ 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._var = None + self.var = var + if integral_scale is not 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) + # final check for parameter bounds self.check_arg_bounds() # additional checks for the optional arguments (provided by user) @@ -500,10 +485,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 @@ -963,24 +944,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 diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index f19eaac3..e6d6b03d 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -19,6 +19,7 @@ __all__ = ["fit_variogram"] +# this should be: set(mod.iso_arg) - set(mod.opt_arg) DEFAULT_PARA = ["var", "len_scale", "nugget"] @@ -217,21 +218,13 @@ 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: 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 @@ -426,9 +419,9 @@ def curve(x, arg1, *args): para_skip = 0 opt_skip = 0 if para["var"]: - var_tmp = args[para_skip] + model.var = args[para_skip] if constrain_sill: - nugget_tmp = sill - var_tmp + 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) @@ -445,12 +438,6 @@ def curve(x, arg1, *args): 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 if is_dir_vario: if anis: model.anis = args[1 - model.dim :] @@ -469,13 +456,9 @@ def _post_fitting(model, para, popt, anis, is_dir_vario): fit_para = {} para_skip = 0 opt_skip = 0 - var_tmp = 0.0 # init value for par in DEFAULT_PARA: 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: @@ -491,9 +474,6 @@ def _post_fitting(model, para, popt, anis, 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/tools.py b/src/gstools/covmodel/tools.py index 9d53de58..1c3c9ca0 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -592,7 +592,6 @@ 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)) 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 Date: Sun, 11 Aug 2024 12:16:38 +0200 Subject: [PATCH 11/28] tests: remove tests for var_raw --- tests/test_covmodel.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index a2729dd6..c42e4774 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -201,15 +201,9 @@ def test_tpl_models(self): self.gamma_x, self.gamma_y, sill=1.1, nugget=False ) self.assertAlmostEqual(model.var, 1.1, delta=1e-5) - # check var_raw handling - model = Model(var_raw=1, len_low=0, integral_scale=10) - var_save = model.var - model.var_raw = 1.1 - self.assertAlmostEqual(model.var, var_save * 1.1) - self.assertAlmostEqual(model.integral_scale, 10) # integral scale is not setable when len_low is not 0 with self.assertRaises(ValueError): - Model(var_raw=1, len_low=5, integral_scale=10) + Model(len_low=5, integral_scale=10) def test_fitting(self): for Model in self.std_cov_models: From 669e9a203a5d40d7b35f11896e64674f3c57e0b2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Sun, 11 Aug 2024 16:58:40 +0200 Subject: [PATCH 12/28] CovModel: remove mechanism for fixed args again --- src/gstools/covmodel/base.py | 36 ++++------------------------------- src/gstools/covmodel/fit.py | 2 -- src/gstools/covmodel/tools.py | 2 -- 3 files changed, 4 insertions(+), 36 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index fce1e64b..4a1ee76c 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -132,8 +132,6 @@ class CovModel: used for the spectrum calculation. Use with caution (Better: Don't!). ``None`` is equivalent to ``{"a": -1, "b": 1, "N": 1000, "h": 0.001}``. Default: :any:`None` - fixed: :class:`set` or :any:`None`, optional - Names of fixed arguments. Default: :any:`None` **opt_arg Optional arguments are covered by these keyword arguments. If present, they are described in the section `Other Parameters`. @@ -155,7 +153,6 @@ def __init__( temporal=False, spatial_dim=None, hankel_kw=None, - fixed=None, **opt_arg, ): # assert, that we use a subclass @@ -164,14 +161,13 @@ def __init__( if not hasattr(self, "variogram"): raise TypeError("Don't instantiate 'CovModel' directly!") - # indicate that arguments are fixed (True after __init__) + # indicator for initialization status (True after __init__) self._init = False # prepare dim setting self._dim = None self._hankel_kw = None self._sft = None # prepare parameters (they are checked in dim setting) - self._fixed = None self._rescale = None self._len_scale = None self._anis = None @@ -196,16 +192,6 @@ def __init__( # optional arguments for the variogram-model set_opt_args(self, opt_arg) - # check fixed - fixed = set(fixed) if fixed is not None else set() - fixed = fixed | self.always_fixed() - valid_fixed = set(self.iso_arg) # | set(["anis"]) - if not fixed <= valid_fixed: - raise ValueError( - f"CovModel: unknown names in 'fixed': {fixed - valid_fixed}" - ) - self._fixed = fixed - # set standard boundaries for variance, len_scale, nugget and opt_arg bounds = self.default_arg_bounds() bounds.update(self.default_opt_arg_bounds()) @@ -213,6 +199,7 @@ 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 @@ -225,9 +212,6 @@ def __init__( self.dim, angles, self.latlon, self.temporal ) - self._var = None - self.var = var - if integral_scale is not None: self.integral_scale = integral_scale @@ -235,13 +219,10 @@ def __init__( self.check_arg_bounds() # additional checks for the optional arguments (provided by user) self.check_opt_arg() - # set fixed bounds after checking original bounds - for arg in self.fixed: - val = getattr(self, arg) - self.set_arg_bounds(check_args=False, **{arg: (val, val, "cc")}) - self._init = True # precision for printing self._prec = 3 + # initialized + self._init = True # one of these functions needs to be overridden def __init_subclass__(cls): @@ -444,10 +425,6 @@ def pykrige_kwargs(self): # methods for optional/default arguments (can be overridden) - def always_fixed(self): - """Provide set of fixed arguments.""" - return set() - def default_opt_arg(self): """Provide default optional arguments by the user. @@ -801,11 +778,6 @@ def check_arg_bounds(self): # bounds properties - @property - def fixed(self): - """:class:`set`: Set with names of fixed arguments.""" - return self._fixed - @property def var_bounds(self): """:class:`list`: Bounds for the variance. diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index e6d6b03d..b40fc845 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -166,8 +166,6 @@ def fit_variogram( The fitted parameters will be instantly set in the model. """ - for para in model.fixed: - para_select[para] = False # preprocess selected parameters para, sill, constrain_sill, anis = _pre_para( model, para_select, sill, anis diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 1c3c9ca0..78a90a7b 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -458,8 +458,6 @@ def set_arg_bounds(model, check_args=True, **kwargs): # if variance needs to be resetted, do this at last var_bnds = [] for arg, bounds in kwargs.items(): - if model._init and arg in model.fixed: - raise ValueError(f"Can't set bounds for fixed argument '{arg}'") if not check_bounds(bounds): raise ValueError( f"Given bounds for '{arg}' are not valid, got: {bounds}" From 01d1d8233ae44fe63ad84e170a0e68d1b296592b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Aug 2024 00:34:57 +0200 Subject: [PATCH 13/28] CovModel: simplify fitting routines --- src/gstools/covmodel/fit.py | 62 +++++++++---------------------------- 1 file changed, 14 insertions(+), 48 deletions(-) diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index b40fc845..b5836c0d 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -19,10 +19,6 @@ __all__ = ["fit_variogram"] -# this should be: set(mod.iso_arg) - set(mod.opt_arg) -DEFAULT_PARA = ["var", "len_scale", "nugget"] - - def fit_variogram( model, x_data, @@ -263,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 @@ -363,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 @@ -376,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]) @@ -408,34 +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"]: - model.var = args[para_skip] - 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 - 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 + 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 :] @@ -453,21 +427,13 @@ 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 - for par in DEFAULT_PARA: + for par in model.iso_arg: if para[par]: 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 :] From 418038f521a4f566b8016990db67547be1f39b69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Aug 2024 11:34:51 +0200 Subject: [PATCH 14/28] CovModel: fix setting integral scale as list of values --- src/gstools/covmodel/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 4a1ee76c..a02a4651 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -1004,7 +1004,7 @@ def integral_scale(self, integral_scale): 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, integral_scale, rtol=1e-3): + 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'!" From d8b037fe994c422fc29c147d3b01ef152e927df5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Aug 2024 12:38:40 +0200 Subject: [PATCH 15/28] no need to set var last anymore --- src/gstools/covmodel/tools.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 78a90a7b..448dc6c1 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -455,18 +455,14 @@ 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 elif arg == "nugget": @@ -481,11 +477,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): From 0b632253efcb6a1b69f76bf72b393ff2b1e94fe1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Aug 2024 12:38:57 +0200 Subject: [PATCH 16/28] fit: add comment --- src/gstools/covmodel/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index b5836c0d..5ba398ef 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -169,7 +169,7 @@ def fit_variogram( 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'") @@ -212,13 +212,13 @@ def fit_variogram( def _pre_para(model, para_select, sill, anis): """Preprocess selected parameters.""" + # 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): setattr(model, par, float(para_select[par])) para_select[par] = False - # remove those that were set to True para_select = {k: v for k, v in para_select.items() if not v} # handling the sill From 21d0e4a27d8b3dce26a8c2956ea398c0db601c50 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 14 Aug 2024 12:12:45 +0200 Subject: [PATCH 17/28] add ratio error class; better arg checking routines --- src/gstools/covmodel/tools.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 448dc6c1..50ef04dd 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -6,6 +6,7 @@ The following classes and functions are provided .. autosummary:: + RatioError AttributeWarning rad_fac set_opt_args @@ -34,6 +35,7 @@ from gstools.tools.misc import list_format __all__ = [ + "RatioError", "AttributeWarning", "rad_fac", "set_opt_args", @@ -52,6 +54,10 @@ ] +class RatioError(Exception): + """Error for invalid ratios in SumModel.""" + + class AttributeWarning(UserWarning): """Attribute warning for CovModel class.""" @@ -311,7 +317,7 @@ def check_bounds(bounds): 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}") @@ -320,6 +326,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]): @@ -333,6 +340,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 @@ -497,17 +513,7 @@ def check_arg_bounds(model): 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): From e6374db261076351920c3fb1857690c558a57b90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Aug 2024 23:26:21 +0200 Subject: [PATCH 18/28] CovModel: better dim setting --- src/gstools/covmodel/tools.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 50ef04dd..a70c8118 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -511,8 +511,6 @@ 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) check_arg_in_bounds(model, arg, error=True) @@ -557,16 +555,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): From dcb003bffec9fa2f4e0c1ebf0e174620a0183bd7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:04:00 +0200 Subject: [PATCH 19/28] add sum_tools submodule --- src/gstools/covmodel/sum_tools.py | 245 ++++++++++++++++++++++++++++++ src/gstools/covmodel/tools.py | 5 - 2 files changed, 245 insertions(+), 5 deletions(-) create mode 100644 src/gstools/covmodel/sum_tools.py diff --git a/src/gstools/covmodel/sum_tools.py b/src/gstools/covmodel/sum_tools.py new file mode 100644 index 00000000..a4a1647a --- /dev/null +++ b/src/gstools/covmodel/sum_tools.py @@ -0,0 +1,245 @@ +""" +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 +""" + +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 a70c8118..52a6e64c 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -6,7 +6,6 @@ The following classes and functions are provided .. autosummary:: - RatioError AttributeWarning rad_fac set_opt_args @@ -54,10 +53,6 @@ ] -class RatioError(Exception): - """Error for invalid ratios in SumModel.""" - - class AttributeWarning(UserWarning): """Attribute warning for CovModel class.""" From 057c8884791059fdf06de95108f8c12d7c3ebd8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:04:21 +0200 Subject: [PATCH 20/28] add SumModel class --- src/gstools/covmodel/base.py | 503 +++++++++++++++++++++++++++++++++++ 1 file changed, 503 insertions(+) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index a02a4651..8e694b26 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -18,9 +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_arg_in_bounds, compare, default_arg_from_bounds, model_repr, @@ -1156,3 +1166,496 @@ def __setattr__(self, name, value): def __repr__(self): """Return String representation.""" return model_repr(self) + + +class SumModel(CovModel): + r"""Class for sums of covariance models. + + Parameters + ---------- + *models + tuple of :any:`CovModel` instances of subclasses to sum. + 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`` + 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`` + 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`` + 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` + rescale : :class:`float` or :any:`None`, optional + Optional rescaling factor to divide the length scale with. + This could be used for unit conversion or rescaling the length scale + to coincide with e.g. the integral scale. + Will be set by each model individually. + 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 + 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 + 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_``. + """ + + def __init__(self, *models, **kwargs): + self._init = False + self._models = [] + add_nug = 0.0 + for mod in models: + if isinstance(mod, type) and issubclass(mod, SumModel): + continue # treat un-init sum-model as nugget model with 0 nugget + if isinstance(mod, SumModel): + self._models += mod.models + add_nug += mod.nugget + elif isinstance(mod, CovModel): + self._models.append(mod) + elif isinstance(mod, type) and issubclass(mod, CovModel): + 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.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) + self._rescale = 1.0 + 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 parameters 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) + # opt arg determining + self._prec = 3 + self._init = True + # set remaining args + for arg, val in kwargs.items(): + setattr(self, arg, val) + if var_set is not None: + self.var = var_set + if len_set is not None: + self.len_scale = len_set + # check for consistency + 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 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, vars): + if len(vars) != len(self): + msg = "SumModel: number of given variances not matching" + raise ValueError(msg) + for mod, var in zip(self.models, vars): + 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): + 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): + 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) From 656125a4c37345b7bd4cea12639e52bb538f31ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:09:43 +0200 Subject: [PATCH 21/28] add pure Nugget model --- src/gstools/__init__.py | 6 +++ src/gstools/covmodel/__init__.py | 7 +++- src/gstools/covmodel/models.py | 67 +++++++++++++++++++++++++++++++- 3 files changed, 78 insertions(+), 2 deletions(-) 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/models.py b/src/gstools/covmodel/models.py index b1a9d68e..80dbd705 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 @@ -27,11 +28,12 @@ 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. From dbdfa27fabcf5eb072ac0c43085c77d6bd008d74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:26:32 +0200 Subject: [PATCH 22/28] CovModel: let the sum magic happen --- src/gstools/covmodel/base.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 8e694b26..68af58d0 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -1163,6 +1163,14 @@ def __setattr__(self, name, value): 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) From 2858f063665ffb2b1d7b5a921ac25609f47c7835 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:35:53 +0200 Subject: [PATCH 23/28] pylint fixes --- src/gstools/covmodel/base.py | 10 ++++++---- src/gstools/covmodel/models.py | 2 +- src/gstools/covmodel/sum_tools.py | 9 +++++---- src/gstools/covmodel/tools.py | 1 - 4 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 68af58d0..a495d80b 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 @@ -1480,11 +1480,11 @@ def vars(self): return [mod.var for mod in self.models] @vars.setter - def vars(self, vars): - if len(vars) != len(self): + 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, vars): + 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) @@ -1555,6 +1555,7 @@ def angles(self, 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)) @@ -1593,6 +1594,7 @@ def spectral_density(self, k): ) def correlation(self, r): + """SumModel correlation function.""" return sum( ( mod.correlation(r) * rat diff --git a/src/gstools/covmodel/models.py b/src/gstools/covmodel/models.py index 80dbd705..7c1cd836 100644 --- a/src/gstools/covmodel/models.py +++ b/src/gstools/covmodel/models.py @@ -22,7 +22,7 @@ JBessel """ -# pylint: disable=C0103, E1101, R0201 +# pylint: disable=C0103, C0302, E1101, R0201 import warnings import numpy as np diff --git a/src/gstools/covmodel/sum_tools.py b/src/gstools/covmodel/sum_tools.py index a4a1647a..c4b9a90e 100644 --- a/src/gstools/covmodel/sum_tools.py +++ b/src/gstools/covmodel/sum_tools.py @@ -17,6 +17,7 @@ sum_model_repr """ +# pylint: disable=W0212 import numpy as np from gstools.covmodel.tools import check_arg_in_bounds @@ -89,10 +90,10 @@ 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) + 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), diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 52a6e64c..69bf963b 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -34,7 +34,6 @@ from gstools.tools.misc import list_format __all__ = [ - "RatioError", "AttributeWarning", "rad_fac", "set_opt_args", From 947619a0428e716e7b3ee78952b631be3c24824d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:40:47 +0200 Subject: [PATCH 24/28] fix doc-string of SumModel for docs --- src/gstools/covmodel/base.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index a495d80b..252bdeb3 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -1356,6 +1356,18 @@ def __init__(self, *models, **kwargs): # check for consistency self.check() + def __init_subclass__(cls): + """Initialize gstools sum-model.""" + _init_subclass(cls) + # overridden functions get standard doc if no new doc was created + ign = ["__", "variogram", "covariance", "cor"] + for att, attr_cls in cls.__dict__.items(): + if any(att.startswith(i) for i in ign) or att not in dir(CovModel): + continue + attr_doc = getattr(CovModel, att).__doc__ + if attr_cls.__doc__ is None: + attr_cls.__doc__ = attr_doc + def __iter__(self): return iter(self.models) From 480fa32e8d6d227ac9e8039e4dabc222951adea0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:57:54 +0200 Subject: [PATCH 25/28] CovModel: add switch to add doc string --- src/gstools/covmodel/base.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 252bdeb3..e080ba20 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -147,6 +147,8 @@ class CovModel: If present, they are described in the section `Other Parameters`. """ + _add_doc = True + def __init__( self, dim=3, @@ -242,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(): @@ -1259,6 +1262,8 @@ class SumModel(CovModel): Also covers ``var_`` and ``length_scale_``. """ + _add_doc = False + def __init__(self, *models, **kwargs): self._init = False self._models = [] @@ -1356,18 +1361,6 @@ def __init__(self, *models, **kwargs): # check for consistency self.check() - def __init_subclass__(cls): - """Initialize gstools sum-model.""" - _init_subclass(cls) - # overridden functions get standard doc if no new doc was created - ign = ["__", "variogram", "covariance", "cor"] - for att, attr_cls in cls.__dict__.items(): - if any(att.startswith(i) for i in ign) or att not in dir(CovModel): - continue - attr_doc = getattr(CovModel, att).__doc__ - if attr_cls.__doc__ is None: - attr_cls.__doc__ = attr_doc - def __iter__(self): return iter(self.models) From 171efbad0cbb1f55523be315b6db053ba8414a44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 16 Aug 2024 00:58:09 +0200 Subject: [PATCH 26/28] Fourier: fix doc --- src/gstools/field/generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py index 89914785..83d13816 100755 --- a/src/gstools/field/generator.py +++ b/src/gstools/field/generator.py @@ -736,7 +736,7 @@ def update(self, model=None, seed=np.nan, period=None, mode_no=None): the seed of the random number generator. If :any:`None`, a random seed is used. If :any:`numpy.nan`, the actual seed will be kept. Default: :any:`numpy.nan` - period : :class:`list` or :any:`None, optional + period : :class:`list` or :any:`None`, optional The spatial periodicity of the field, is often the domain size. mode_no : :class:`list` or :any:`None`, optional Number of Fourier modes per dimension. From e7838cf023ceff961054426e7b6eb0a169d32601 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 19 Aug 2024 11:30:24 +0200 Subject: [PATCH 27/28] SumModel: models need to either be all instances or all subclasses of CovModel --- src/gstools/covmodel/base.py | 57 +++++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 11 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index e080ba20..dd7ab2ee 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -1182,15 +1182,32 @@ def __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. + 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`` + 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 @@ -1206,14 +1223,14 @@ class SumModel(CovModel): 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`` + 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`` + 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. @@ -1234,20 +1251,20 @@ class SumModel(CovModel): 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 + 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` + 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 + 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 @@ -1259,7 +1276,8 @@ class SumModel(CovModel): 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_``. + Also covers ``var_`` and ``length_scale_`` but they should preferably be + set by ``vars`` and ``length_scales``. """ _add_doc = False @@ -1268,15 +1286,31 @@ 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." @@ -1315,7 +1349,7 @@ def __init__(self, *models, **kwargs): self._nugget = float( kwargs.pop( "nugget", - sum((mod.nugget for mod in self.models), 0.0) + add_nug, + sum((mod.nugget for mod in self.models), 0) + add_nug, ) ) for mod in self.models: @@ -1339,7 +1373,7 @@ def __init__(self, *models, **kwargs): self._angles = set_model_angles( self.dim, angles, self.latlon, self.temporal ) - # prepare parameters boundaries + # prepare parameter boundaries self._var_bounds = None self._len_scale_bounds = None self._nugget_bounds = None @@ -1348,17 +1382,18 @@ def __init__(self, *models, **kwargs): bounds = self.default_arg_bounds() bounds.update(self.default_opt_arg_bounds()) self.set_arg_bounds(check_args=False, **bounds) - # opt arg determining + # 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 for consistency + # check consistency of sub models self.check() def __iter__(self): From 0e35a222ef110af52f59ec2048b0be0d4c3cf2de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 19 Aug 2024 11:40:08 +0200 Subject: [PATCH 28/28] SumModel: sum models have a constant rescale factor of 1 since they are derived from correlation not cor --- src/gstools/covmodel/base.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index dd7ab2ee..38a47491 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -1193,6 +1193,7 @@ class SumModel(CovModel): 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 ---------- @@ -1235,12 +1236,6 @@ class SumModel(CovModel): If given, ``len_scale`` will be ignored and recalculated, so that the integral scale of the model matches the given one. Default: :any:`None` - rescale : :class:`float` or :any:`None`, optional - Optional rescaling factor to divide the length scale with. - This could be used for unit conversion or rescaling the length scale - to coincide with e.g. the integral scale. - Will be set by each model individually. - 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 @@ -1364,7 +1359,6 @@ def __init__(self, *models, **kwargs): 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) - self._rescale = 1.0 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 @@ -1477,6 +1471,11 @@ 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."""