diff --git a/examples/05_kriging/10_simple_collocated_cokriging.py b/examples/05_kriging/10_simple_collocated_cokriging.py new file mode 100644 index 000000000..0e739d1fd --- /dev/null +++ b/examples/05_kriging/10_simple_collocated_cokriging.py @@ -0,0 +1,86 @@ +r""" +Simple Collocated Cokriging +---------------------------- + +Simple collocated cokriging uses secondary data at the estimation location +to improve the primary variable estimate. + +This example demonstrates the new correlogram-based API using MarkovModel1, +which encapsulates the Markov Model I (MM1) cross-covariance structure. + +Example +^^^^^^^ + +Here we compare Simple Kriging with Simple Collocated Cokriging using the +new MarkovModel1 correlogram. +""" + +import matplotlib.pyplot as plt +import numpy as np + +from gstools import Gaussian, MarkovModel1, krige +from gstools.cokriging import SimpleCollocated + +# condtions +np.random.seed(4) +cond_pos = np.array([0.5, 2.1, 3.8, 6.2, 13.5]) +cond_val = np.array([0.8, 1.2, 1.8, 2.1, 1.4]) +gridx = np.linspace(0.0, 15.0, 151) +model = Gaussian(dim=1, var=0.5, len_scale=2.0) + +############################################################################### +# Generate correlated secondary data + +sec_pos = np.linspace(0, 15, 31) +primary_trend = np.interp(sec_pos, cond_pos, cond_val) +gap_feature = -1.6 * np.exp(-(((sec_pos - 10.0) / 2.0) ** 2)) +gap_feature2 = -0.95 * np.exp(-(((sec_pos - 4.0) / 2.0) ** 2)) +sec_val = 0.99 * primary_trend + gap_feature + gap_feature2 + +sec_grid = np.interp(gridx, sec_pos, sec_val) +sec_at_primary = np.interp(cond_pos, sec_pos, sec_val) + +############################################################################### +# Simple Kriging and Simple Collocated Cokriging + +sk = krige.Simple(model, cond_pos=cond_pos, cond_val=cond_val, mean=1.0) +sk_field, sk_var = sk(gridx, return_var=True) + +# Compute cross-correlation from data +cross_corr = np.corrcoef(cond_val, sec_at_primary)[0, 1] + +# Create MarkovModel1 correlogram (NEW API) +correlogram = MarkovModel1( + primary_model=model, + cross_corr=cross_corr, + secondary_var=np.var(sec_val), + primary_mean=1.0, + secondary_mean=np.mean(sec_val), +) + +# Simple Collocated Cokriging with new API +scck = SimpleCollocated(correlogram, cond_pos=cond_pos, cond_val=cond_val) +scck_field, scck_var = scck(gridx, secondary_data=sec_grid, return_var=True) + +############################################################################### + +fig, ax = plt.subplots(1, 2, figsize=(10, 3.5)) + +ax[0].scatter(cond_pos, cond_val, color="red", label="Primary data") +ax[0].scatter( + cond_pos, + sec_at_primary, + color="blue", + marker="s", + label="Secondary at primary", +) +ax[0].plot(sec_pos, sec_val, "b-", alpha=0.6, label="Secondary data") +ax[0].legend() + +ax[1].plot(gridx, sk_field, label="Simple Kriging") +ax[1].plot(gridx, scck_field, label="Simple Collocated Cokriging") +ax[1].scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") +ax[1].legend() + +plt.tight_layout() +plt.show() diff --git a/examples/05_kriging/11_intrinsic_collocated_cokriging.py b/examples/05_kriging/11_intrinsic_collocated_cokriging.py new file mode 100644 index 000000000..fc6c2b2f5 --- /dev/null +++ b/examples/05_kriging/11_intrinsic_collocated_cokriging.py @@ -0,0 +1,95 @@ +r""" +Intrinsic Collocated Cokriging +------------------------------- + +Intrinsic Collocated Cokriging (ICCK) improves variance estimation +compared to Simple Collocated Cokriging. + +This example demonstrates the new correlogram-based API using MarkovModel1. + +The variance formula is: + +.. math:: \sigma^2_{ICCK} = (1 - \rho_0^2) \cdot \sigma^2_{SK} + +Example +^^^^^^^ + +Here we compare Simple Kriging with Intrinsic Collocated Cokriging using the +new MarkovModel1 correlogram. +""" + +import matplotlib.pyplot as plt +import numpy as np + +from gstools import Gaussian, MarkovModel1, krige +from gstools.cokriging import IntrinsicCollocated + +# condtions +np.random.seed(4) +cond_pos = np.array([0.5, 2.1, 3.8, 6.2, 13.5]) +cond_val = np.array([0.8, 1.2, 1.8, 2.1, 1.4]) +gridx = np.linspace(0.0, 15.0, 151) +model = Gaussian(dim=1, var=0.5, len_scale=2.0) + +############################################################################### +# Generate correlated secondary data + +sec_pos = np.linspace(0, 15, 31) +primary_trend = np.interp(sec_pos, cond_pos, cond_val) +gap_feature = -1.6 * np.exp(-(((sec_pos - 10.0) / 2.0) ** 2)) +gap_feature2 = -0.95 * np.exp(-(((sec_pos - 4.0) / 2.0) ** 2)) +sec_val = 0.99 * primary_trend + gap_feature + gap_feature2 + +sec_grid = np.interp(gridx, sec_pos, sec_val) +sec_at_primary = np.interp(cond_pos, sec_pos, sec_val) + +############################################################################### +# Simple Kriging and Intrinsic Collocated Cokriging + +sk = krige.Simple(model, cond_pos=cond_pos, cond_val=cond_val, mean=1.0) +sk_field, sk_var = sk(gridx, return_var=True) + +# Compute cross-correlation from data +cross_corr = np.corrcoef(cond_val, sec_at_primary)[0, 1] + +# Create MarkovModel1 correlogram (NEW API) +correlogram = MarkovModel1( + primary_model=model, + cross_corr=cross_corr, + secondary_var=np.var(sec_val), + primary_mean=1.0, + secondary_mean=np.mean(sec_val), +) + +# Intrinsic Collocated Cokriging with new API +icck = IntrinsicCollocated( + correlogram, + cond_pos=cond_pos, + cond_val=cond_val, + secondary_cond_pos=cond_pos, + secondary_cond_val=sec_at_primary, +) +icck_field, icck_var = icck(gridx, secondary_data=sec_grid, return_var=True) + +############################################################################### + +fig, ax = plt.subplots(1, 2, figsize=(10, 3.5)) + +ax[0].scatter(cond_pos, cond_val, color="red", label="Primary data") +ax[0].scatter( + cond_pos, + sec_at_primary, + color="blue", + marker="s", + label="Secondary at primary", +) +ax[0].plot(sec_pos, sec_val, "b-", alpha=0.6, label="Secondary data") +ax[0].legend() + +ax[1].plot(gridx, sk_field, label="Simple Kriging") +ax[1].plot(gridx, icck_field, label="Intrinsic Collocated Cokriging") +ax[1].scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") +ax[1].legend() + +plt.tight_layout() +plt.show() diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 4d12007c9..7cc5e53f8 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -19,6 +19,7 @@ field variogram krige + cokriging random tools transform @@ -36,6 +37,17 @@ .. autosummary:: Krige +Cokriging +^^^^^^^^^ +Collocated cokriging methods for multivariate estimation + +.. currentmodule:: gstools.cokriging + +.. autosummary:: + SimpleCollocated + IntrinsicCollocated + MarkovModel1 + Spatial Random Field ^^^^^^^^^^^^^^^^^^^^ Classes for (conditioned) random field generation @@ -134,17 +146,23 @@ """ # Hooray! -from gstools import ( +from gstools import ( # noqa: I001 config, covmodel, field, krige, + cokriging, normalizer, random, tools, transform, variogram, ) +from gstools.cokriging import ( + IntrinsicCollocated, + MarkovModel1, + SimpleCollocated, +) from gstools.covmodel import ( Circular, CovModel, @@ -199,7 +217,15 @@ __version__ = "0.0.0.dev0" __all__ = ["__version__"] -__all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] +__all__ += [ + "covmodel", + "field", + "variogram", + "krige", + "cokriging", + "random", + "tools", +] __all__ += ["transform", "normalizer", "config"] __all__ += [ "CovModel", @@ -234,6 +260,9 @@ __all__ += [ "Krige", + "SimpleCollocated", + "IntrinsicCollocated", + "MarkovModel1", "SRF", "CondSRF", "PGS", diff --git a/src/gstools/cokriging/__init__.py b/src/gstools/cokriging/__init__.py new file mode 100644 index 000000000..ba037ba6c --- /dev/null +++ b/src/gstools/cokriging/__init__.py @@ -0,0 +1,36 @@ +""" +GStools subpackage providing cokriging. + +.. currentmodule:: gstools.cokriging + +Cokriging Classes +^^^^^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: + + CollocatedCokriging + SimpleCollocated + IntrinsicCollocated + +Correlogram Models +^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: + + Correlogram + MarkovModel1 +""" + +from gstools.cokriging.base import CollocatedCokriging +from gstools.cokriging.correlogram import Correlogram, MarkovModel1 +from gstools.cokriging.methods import IntrinsicCollocated, SimpleCollocated + +__all__ = [ + "CollocatedCokriging", + "SimpleCollocated", + "IntrinsicCollocated", + "Correlogram", + "MarkovModel1", +] diff --git a/src/gstools/cokriging/base.py b/src/gstools/cokriging/base.py new file mode 100644 index 000000000..44eb9697e --- /dev/null +++ b/src/gstools/cokriging/base.py @@ -0,0 +1,330 @@ +""" +GStools subpackage providing collocated cokriging. + +.. currentmodule:: gstools.cokriging.base + +The following classes are provided + +.. autosummary:: + CollocatedCokriging +""" + +import numpy as np + +from gstools.cokriging.correlogram import Correlogram +from gstools.krige.base import Krige + +__all__ = ["CollocatedCokriging"] + + +class CollocatedCokriging(Krige): + """ + Collocated cokriging base class using Correlogram models. + + Collocated cokriging uses secondary data at the estimation location + to improve the primary variable estimate. The cross-covariance structure + is defined by a :any:`Correlogram` object (e.g., :any:`MarkovModel1`). + + Two algorithms are supported: Simple Collocated ("simple") uses only + collocated secondary at the estimation point, while Intrinsic Collocated + ("intrinsic") additionally uses secondary data at all primary locations + for more accurate variance estimation. + + Parameters + ---------- + correlogram : :any:`Correlogram` + Correlogram object defining the cross-covariance structure between + primary and secondary variables (e.g., :any:`MarkovModel1`). + cond_pos : :class:`list` + tuple, containing the given condition positions (x, [y, z]) + cond_val : :class:`numpy.ndarray` + the values of the primary variable conditions (nan values will be ignored) + algorithm : :class:`str` + Cokriging algorithm to use. Either "simple" (SCCK) or "intrinsic" (ICCK). + secondary_cond_pos : :class:`list`, optional + tuple, containing secondary variable condition positions (only for ICCK) + secondary_cond_val : :class:`numpy.ndarray`, optional + values of secondary variable at primary locations (only for ICCK) + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the input data to gain normality. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signature: f(x, [y, z, ...]) + This is used for detrended kriging, where the trend is subtracted + from the conditions before kriging is applied. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. + exact : :class:`bool`, optional + Whether the interpolator should reproduce the exact input values. + If `False`, `cond_err` is interpreted as measurement error + at the conditioning points and the result will be more smooth. + Default: False + cond_err : :class:`str`, :class:`float` or :class:`list`, optional + The measurement error at the conditioning points. + Either "nugget" to apply the model-nugget, a single value applied to + all points or an array with individual values for each point. + The measurement error has to be <= nugget. + The "exact=True" variant only works with "cond_err='nugget'". + Default: "nugget" + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: True + pseudo_inv_type : :class:`str` or :any:`callable`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `"pinv"`: use `pinv` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values + + If you want to use another routine to invert the kriging matrix, + you can pass a callable which takes a matrix and returns the inverse. + Default: `"pinv"` + fit_normalizer : :class:`bool`, optional + Whether to fit the data-normalizer to the given conditioning data. + Default: False + fit_variogram : :class:`bool`, optional + Whether to fit the given variogram model to the data. + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with + standard bins provided by the :any:`standard_bins` routine. + Default: False + + References + ---------- + .. [Samson2020] Samson, M., & Deutsch, C. V. (2020). Collocated Cokriging. + In J. L. Deutsch (Ed.), Geostatistics Lessons. Retrieved from + http://geostatisticslessons.com/lessons/collocatedcokriging + .. [Wackernagel2003] Wackernagel, H. Multivariate Geostatistics, + Springer, Berlin, 2003. + """ + + def __init__( + self, + correlogram, + cond_pos, + cond_val, + algorithm, + secondary_cond_pos=None, + secondary_cond_val=None, + normalizer=None, + trend=None, + exact=False, + cond_err="nugget", + pseudo_inv=True, + pseudo_inv_type="pinv", + fit_normalizer=False, + fit_variogram=False, + ): + # Validate correlogram + if not isinstance(correlogram, Correlogram): + raise TypeError( + f"correlogram must be a Correlogram instance, got {type(correlogram)}" + ) + self.correlogram = correlogram + + # validate algorithm parameter + if algorithm not in ["simple", "intrinsic"]: + raise ValueError("algorithm must be 'simple' or 'intrinsic'") + self.algorithm = algorithm + + # handle secondary conditioning data (required for intrinsic) + if algorithm == "intrinsic": + if secondary_cond_pos is None or secondary_cond_val is None: + raise ValueError( + "secondary_cond_pos and secondary_cond_val required for ICCK" + ) + self.secondary_cond_pos = secondary_cond_pos + self.secondary_cond_val = np.asarray( + secondary_cond_val, dtype=np.double + ) + + if len(self.secondary_cond_val) != len(cond_val): + raise ValueError( + "secondary_cond_val must have same length as primary cond_val" + ) + else: + self.secondary_cond_pos = None + self.secondary_cond_val = None + + # initialize as simple kriging (unbiased=False) + super().__init__( + model=correlogram.primary_model, + cond_pos=cond_pos, + cond_val=cond_val, + mean=correlogram.primary_mean, + unbiased=False, # Simple kriging base + normalizer=normalizer, + trend=trend, + exact=exact, + cond_err=cond_err, + pseudo_inv=pseudo_inv, + pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, + ) + + def __call__(self, pos=None, secondary_data=None, **kwargs): + """ + Generate the collocated cokriging field. + + The field is saved as `self.field` and is also returned. + The error variance is saved as `self.krige_var` and is also returned. + + Parameters + ---------- + pos : :class:`list` + the position tuple, containing main direction and transversal + directions (x, [y, z]) + secondary_data : :class:`numpy.ndarray` + Secondary variable values at the given positions. + **kwargs + Keyword arguments passed to Krige.__call__. + + Returns + ------- + field : :class:`numpy.ndarray` + the collocated cokriging field + krige_var : :class:`numpy.ndarray`, optional + the collocated cokriging error variance + (if return_var is True) + """ + if secondary_data is None: + raise ValueError( + "secondary_data required for collocated cokriging" + ) + + user_return_var = kwargs.get("return_var", True) + # always get variance for weight calculation + kwargs_with_var = kwargs.copy() + kwargs_with_var["return_var"] = True + # get simple kriging results + sk_field, sk_var = super().__call__(pos=pos, **kwargs_with_var) + secondary_data = np.asarray(secondary_data, dtype=np.double) + + # apply algorithm-specific post-processing + if self.algorithm == "simple": + cokriging_field, cokriging_var = self._apply_simple_collocated( + sk_field, sk_var, secondary_data, user_return_var + ) + elif self.algorithm == "intrinsic": + cokriging_field, cokriging_var = self._apply_intrinsic_collocated( + sk_field, sk_var, secondary_data, user_return_var + ) + else: + raise ValueError(f"Unknown algorithm: {self.algorithm}") + + if user_return_var: + return cokriging_field, cokriging_var + return cokriging_field + + def _apply_simple_collocated( + self, sk_field, sk_var, secondary_data, return_var + ): + """Apply simple collocated cokriging.""" + C_Z0, C_Y0, C_YZ0 = self._compute_covariances() + k = C_YZ0 / C_Z0 + + # compute collocated weight + numerator = k * sk_var + denominator = C_Y0 - (k**2) * (C_Z0 - sk_var) + collocated_weights = np.where( + np.abs(denominator) < 1e-15, 0.0, numerator / denominator + ) + + # apply collocated cokriging estimator + scck_field = ( + sk_field * (1 - k * collocated_weights) + + collocated_weights + * (secondary_data - self.correlogram.secondary_mean) + + k * collocated_weights * self.mean + ) + + if return_var: + # simple collocated variance + scck_variance = sk_var * (1 - collocated_weights * k) + scck_variance = np.maximum(0.0, scck_variance) + else: + scck_variance = None + return scck_field, scck_variance + + def _apply_intrinsic_collocated( + self, sk_field, sk_var, secondary_data, return_var + ): + """ + Apply intrinsic collocated cokriging. + + Adds the collocated secondary contribution at estimation locations + and computes ICCK variance. + + Note: The secondary-at-primary contribution is already added during + the kriging solve in _summate(). + """ + # apply collocated secondary contribution + collocated_contribution = self._lambda_Y0 * ( + secondary_data - self.correlogram.secondary_mean + ) + icck_field = sk_field + collocated_contribution + + # compute intrinsic variance + if return_var: + C_Z0, C_Y0, C_YZ0 = self._compute_covariances() + if C_Y0 * C_Z0 < 1e-15: + rho_squared = 0.0 + else: + rho_squared = (C_YZ0**2) / (C_Y0 * C_Z0) + icck_var = (1.0 - rho_squared) * sk_var + icck_var = np.maximum(0.0, icck_var) + else: + icck_var = None + return icck_field, icck_var + + def _summate(self, field, krige_var, c_slice, k_vec, return_var): + """Apply intrinsic collocated cokriging during kriging solve.""" + if self.algorithm == "simple": + super()._summate(field, krige_var, c_slice, k_vec, return_var) + return + + elif self.algorithm == "intrinsic": + sk_weights = self._krige_mat @ k_vec + C_Z0, C_Y0, C_YZ0 = self._compute_covariances() + + if abs(C_YZ0) < 1e-15: + self._lambda_Y0 = 0.0 + self._secondary_at_primary = 0.0 + super()._summate(field, krige_var, c_slice, k_vec, return_var) + return + + lambda_weights = sk_weights[: self.cond_no] + mu_weights = -(C_YZ0 / C_Y0) * lambda_weights + lambda_Y0 = C_YZ0 / C_Y0 + + secondary_residuals = ( + self.secondary_cond_val - self.correlogram.secondary_mean + ) + if sk_weights.ndim == 1: + secondary_at_primary = np.sum(mu_weights * secondary_residuals) + else: + secondary_at_primary = np.sum( + mu_weights * secondary_residuals[:, None], axis=0 + ) + + self._lambda_Y0 = lambda_Y0 + self._secondary_at_primary = secondary_at_primary + + super()._summate(field, krige_var, c_slice, k_vec, return_var) + field[c_slice] += secondary_at_primary + else: + raise ValueError(f"Unknown algorithm: {self.algorithm}") + + def _compute_covariances(self): + """ + Compute covariances at zero lag. + + Delegates to the correlogram object. + """ + return self.correlogram.compute_covariances() diff --git a/src/gstools/cokriging/correlogram/__init__.py b/src/gstools/cokriging/correlogram/__init__.py new file mode 100644 index 000000000..5685c1a42 --- /dev/null +++ b/src/gstools/cokriging/correlogram/__init__.py @@ -0,0 +1,36 @@ +""" +GStools subpackage providing correlogram models for collocated cokriging. + +.. currentmodule:: gstools.cokriging.correlogram + +Correlogram models define the cross-covariance structure between primary +and secondary variables in collocated cokriging. + +Base Class +^^^^^^^^^^ + +.. autosummary:: + :toctree: + + Correlogram + +Markov Models +^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: + + MarkovModel1 + +Future Models +^^^^^^^^^^^^^ + +Planned implementations: + - MarkovModel2: Uses secondary variable's spatial structure + - LinearModelCoregionalization: Full multivariate model +""" + +from gstools.cokriging.correlogram.base import Correlogram +from gstools.cokriging.correlogram.markov import MarkovModel1 + +__all__ = ["Correlogram", "MarkovModel1"] diff --git a/src/gstools/cokriging/correlogram/base.py b/src/gstools/cokriging/correlogram/base.py new file mode 100644 index 000000000..55e464c74 --- /dev/null +++ b/src/gstools/cokriging/correlogram/base.py @@ -0,0 +1,130 @@ +from abc import ABC, abstractmethod + +__all__ = ["Correlogram"] + + +class Correlogram(ABC): + """ + Abstract base class for cross-covariance models in collocated cokriging. + + A correlogram encapsulates the spatial relationship between primary and + secondary variables, including their cross-covariance structure and + statistical parameters (means, variances). + + This design allows for different cross-covariance models (MM1, MM2, etc.) + to be implemented as separate classes, making the cokriging framework + extensible and future-proof. + + Parameters + ---------- + primary_model : :any:`CovModel` + Covariance model for the primary variable. + cross_corr : :class:`float` + Cross-correlation coefficient between primary and secondary variables + at zero lag (collocated). Must be in [-1, 1]. + secondary_var : :class:`float` + Variance of the secondary variable. Must be positive. + primary_mean : :class:`float`, optional + Mean value of the primary variable. Default: 0.0 + secondary_mean : :class:`float`, optional + Mean value of the secondary variable. Default: 0.0 + + Attributes + ---------- + primary_model : :any:`CovModel` + The primary variable's covariance model. + cross_corr : :class:`float` + Cross-correlation at zero lag. + secondary_var : :class:`float` + Secondary variable variance. + primary_mean : :class:`float` + Primary variable mean. + secondary_mean : :class:`float` + Secondary variable mean. + + Notes + ----- + Subclasses must implement :any:`compute_covariances` and + :any:`cross_covariance` to define the cross-covariance structure. + """ + + def __init__( + self, + primary_model, + cross_corr, + secondary_var, + primary_mean=0.0, + secondary_mean=0.0, + ): + """Initialize the correlogram with spatial and statistical parameters.""" + self.primary_model = primary_model + self.cross_corr = float(cross_corr) + self.secondary_var = float(secondary_var) + self.primary_mean = float(primary_mean) + self.secondary_mean = float(secondary_mean) + + # Validate parameters + self._validate() + + def _validate(self): + """ + Validate correlogram parameters. + + Raises + ------ + ValueError + If cross_corr is not in [-1, 1] or secondary_var is not positive. + """ + if not -1.0 <= self.cross_corr <= 1.0: + raise ValueError( + f"cross_corr must be in [-1, 1], got {self.cross_corr}" + ) + + if self.secondary_var <= 0: + raise ValueError( + f"secondary_var must be positive, got {self.secondary_var}" + ) + + @abstractmethod + def compute_covariances(self): + """ + Compute covariances at zero lag. + + Returns + ------- + C_Z0 : :class:`float` + Primary variable variance :math:`C_Z(0)`. + C_Y0 : :class:`float` + Secondary variable variance :math:`C_Y(0)`. + C_YZ0 : :class:`float` + Cross-covariance between primary and secondary at zero lag + :math:`C_{YZ}(0)`. + + Notes + ----- + This method defines how the cross-covariance at zero lag is computed + from the cross-correlation and variances. Different correlogram models + may use different formulas. + """ + + @abstractmethod + def cross_covariance(self, h): + """ + Compute cross-covariance :math:`C_{YZ}(h)` at distance :math:`h`. + + Parameters + ---------- + h : :class:`float` or :class:`numpy.ndarray` + Distance(s) at which to compute cross-covariance. + + Returns + ------- + C_YZ_h : :class:`float` or :class:`numpy.ndarray` + Cross-covariance at distance :math:`h`. + + Notes + ----- + This is the key method that differentiates correlogram models. + For example, MM1 uses the primary variable's spatial structure + while MM2 would use the secondary variable's structure. + """ diff --git a/src/gstools/cokriging/correlogram/markov.py b/src/gstools/cokriging/correlogram/markov.py new file mode 100644 index 000000000..099b2bb2c --- /dev/null +++ b/src/gstools/cokriging/correlogram/markov.py @@ -0,0 +1,118 @@ +""" +GStools subpackage providing Markov model correlograms. + +.. currentmodule:: gstools.cokriging.correlogram.markov + +The following classes are provided + +.. autosummary:: + MarkovModel1 +""" + +import numpy as np + +from gstools.cokriging.correlogram.base import Correlogram + +__all__ = ["MarkovModel1"] + + +class MarkovModel1(Correlogram): + """ + Markov Model I (MM1) correlogram for collocated cokriging. + + The Markov Model I assumes that the cross-covariance between primary + and secondary variables follows the primary variable's spatial structure: + + .. math:: + C_{YZ}(h) = \\frac{C_{YZ}(0)}{C_Z(0)} \\cdot C_Z(h) + + where :math:`C_{YZ}(h)` is the cross-covariance at distance h, + :math:`C_{YZ}(0)` is the cross-covariance at zero lag, + :math:`C_Z(h)` is the primary variable's covariance at distance h, + and :math:`C_Z(0)` is the primary variable's variance. + + This implies that both variables share the same spatial correlation + structure: :math:`\\rho_Y(h) = \\rho_Z(h)`. + + Parameters + ---------- + primary_model : :any:`CovModel` + Covariance model for the primary variable (Z). This defines the + spatial structure that both variables are assumed to share. + cross_corr : :class:`float` + Cross-correlation coefficient :math:`\\rho_{YZ}(0)` at zero lag. + Must be in [-1, 1]. Computed as: + :math:`\\rho_{YZ}(0) = C_{YZ}(0) / \\sqrt{C_Y(0) \\cdot C_Z(0)}` + secondary_var : :class:`float` + Variance of the secondary variable :math:`C_Y(0)`. Must be positive. + primary_mean : :class:`float`, optional + Mean value of the primary variable :math:`m_Z`. Default: 0.0 + secondary_mean : :class:`float`, optional + Mean value of the secondary variable :math:`m_Y`. Default: 0.0 + + Attributes + ---------- + primary_model : :any:`CovModel` + The primary variable's covariance model. + cross_corr : :class:`float` + Cross-correlation at zero lag. + secondary_var : :class:`float` + Secondary variable variance. + primary_mean : :class:`float` + Primary variable mean. + secondary_mean : :class:`float` + Secondary variable mean. + + References + ---------- + .. [Samson2020] Samson, M., & Deutsch, C. V. (2020). Collocated Cokriging. + In J. L. Deutsch (Ed.), Geostatistics Lessons. Retrieved from + http://geostatisticslessons.com/lessons/collocatedcokriging + .. [Wackernagel2003] Wackernagel, H. Multivariate Geostatistics, + Springer, Berlin, 2003. + """ + + def compute_covariances(self): + """ + Compute covariances at zero lag using MM1 formula. + + Returns + ------- + C_Z0 : :class:`float` + Primary variable variance (sill of primary model). + C_Y0 : :class:`float` + Secondary variable variance (as specified). + C_YZ0 : :class:`float` + Cross-covariance at zero lag, computed as: + :math:`C_{YZ}(0) = \\rho_{YZ}(0) \\cdot \\sqrt{C_Y(0) \\cdot C_Z(0)}` + """ + C_Z0 = self.primary_model.sill + C_Y0 = self.secondary_var + C_YZ0 = self.cross_corr * np.sqrt(C_Z0 * C_Y0) + return C_Z0, C_Y0, C_YZ0 + + def cross_covariance(self, h): + """ + Compute cross-covariance at distance h using MM1 formula. + + Parameters + ---------- + h : :class:`float` or :class:`numpy.ndarray` + Distance(s) at which to compute cross-covariance. + + Returns + ------- + C_YZ_h : :class:`float` or :class:`numpy.ndarray` + Cross-covariance at distance h, computed using MM1: + :math:`C_{YZ}(h) = \\frac{C_{YZ}(0)}{C_Z(0)} \\cdot C_Z(h)` + """ + C_Z0, C_Y0, C_YZ0 = self.compute_covariances() + + # Handle edge case: zero primary variance + if C_Z0 < 1e-15: + return np.zeros_like(h) if isinstance(h, np.ndarray) else 0.0 + + # MM1 formula: C_YZ(h) = (C_YZ(0) / C_Z(0)) * C_Z(h) + k = C_YZ0 / C_Z0 + C_Z_h = self.primary_model.covariance(h) + return k * C_Z_h diff --git a/src/gstools/cokriging/methods.py b/src/gstools/cokriging/methods.py new file mode 100644 index 000000000..81f111eff --- /dev/null +++ b/src/gstools/cokriging/methods.py @@ -0,0 +1,246 @@ +""" +GStools subpackage providing cokriging methods. + +.. currentmodule:: gstools.cokriging.methods + +The following classes are provided + +.. autosummary:: + SimpleCollocated + IntrinsicCollocated +""" + +from gstools.cokriging.base import CollocatedCokriging +from gstools.cokriging.correlogram import Correlogram + +__all__ = ["SimpleCollocated", "IntrinsicCollocated"] + + +class SimpleCollocated(CollocatedCokriging): + """ + Simple collocated cokriging. + + Simple collocated cokriging extends simple kriging by incorporating + secondary variable data at the estimation location only. + + This method can produce variance inflation in some cases. + For accurate variance estimation, use :any:`IntrinsicCollocated` instead. + + Parameters + ---------- + correlogram : :any:`Correlogram` + Correlogram object defining the cross-covariance structure. + Typically a :any:`MarkovModel1` instance. + cond_pos : :class:`list` + tuple, containing the given condition positions (x, [y, z]) + cond_val : :class:`numpy.ndarray` + the values of the conditions (nan values will be ignored) + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the input data to gain normality. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signature: f(x, [y, z, ...]) + This is used for detrended kriging, where the trended is subtracted + from the conditions before kriging is applied. + This can be used for regression kriging, where the trend function + is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. + exact : :class:`bool`, optional + Whether the interpolator should reproduce the exact input values. + If `False`, `cond_err` is interpreted as measurement error + at the conditioning points and the result will be more smooth. + Default: False + cond_err : :class:`str`, :class :class:`float` or :class:`list`, optional + The measurement error at the conditioning points. + Either "nugget" to apply the model-nugget, a single value applied to + all points or an array with individual values for each point. + The measurement error has to be <= nugget. + The "exact=True" variant only works with "cond_err='nugget'". + Default: "nugget" + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: True + pseudo_inv_type : :class:`str` or :any:`callable`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `"pinv"`: use `pinv` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values + + If you want to use another routine to invert the kriging matrix, + you can pass a callable which takes a matrix and returns the inverse. + Default: `"pinv"` + fit_normalizer : :class:`bool`, optional + Whether to fit the data-normalizer to the given conditioning data. + Default: False + fit_variogram : :class:`bool`, optional + Whether to fit the given variogram model to the data. + Default: False + + References + ---------- + .. [Samson2020] Samson, M., & Deutsch, C. V. (2020). Collocated Cokriging. + In J. L. Deutsch (Ed.), Geostatistics Lessons. Retrieved from + http://geostatisticslessons.com/lessons/collocatedcokriging + .. [Wackernagel2003] Wackernagel, H. Multivariate Geostatistics, + Springer, Berlin, 2003. + """ + + def __init__( + self, + correlogram, + cond_pos, + cond_val, + normalizer=None, + trend=None, + exact=False, + cond_err="nugget", + pseudo_inv=True, + pseudo_inv_type="pinv", + fit_normalizer=False, + fit_variogram=False, + ): + # Check if correlogram is actually a Correlogram object + if not isinstance(correlogram, Correlogram): + raise TypeError( + f"First argument must be a Correlogram instance. " + f"Got {type(correlogram).__name__}." + ) + + # Initialize using base class with simple collocated algorithm + super().__init__( + correlogram=correlogram, + cond_pos=cond_pos, + cond_val=cond_val, + algorithm="simple", + normalizer=normalizer, + trend=trend, + exact=exact, + cond_err=cond_err, + pseudo_inv=pseudo_inv, + pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, + ) + + +class IntrinsicCollocated(CollocatedCokriging): + """ + Intrinsic collocated cokriging. + + Intrinsic collocated cokriging extends simple kriging by incorporating + secondary variable data at both the estimation location and at all + primary conditioning locations. + + This method provides accurate variance estimation that eliminates the + variance inflation issue of :any:`SimpleCollocated`, at the cost of + requiring secondary data at all primary locations. + + Parameters + ---------- + correlogram : :any:`Correlogram` + Correlogram object defining the cross-covariance structure. + Typically a :any:`MarkovModel1` instance. + cond_pos : :class:`list` + tuple, containing the given condition positions (x, [y, z]) + cond_val : :class:`numpy.ndarray` + the values of the primary variable conditions (nan values will be ignored) + secondary_cond_pos : :class:`list` + tuple, containing the secondary variable condition positions (x, [y, z]) + secondary_cond_val : :class:`numpy.ndarray` + the values of the secondary variable conditions at primary locations + normalizer : :any:`None` or :any:`Normalizer`, optional + Normalizer to be applied to the input data to gain normality. + The default is None. + trend : :any:`None` or :class:`float` or :any:`callable`, optional + A callable trend function. Should have the signature: f(x, [y, z, ...]) + This is used for detrended kriging, where the trended is subtracted + from the conditions before kriging is applied. + This can be used for regression kriging, where the trend function + is determined by an external regression algorithm. + If no normalizer is applied, this behaves equal to 'mean'. + The default is None. + exact : :class:`bool`, optional + Whether the interpolator should reproduce the exact input values. + If `False`, `cond_err` is interpreted as measurement error + at the conditioning points and the result will be more smooth. + Default: False + cond_err : :class:`str`, :class :class:`float` or :class:`list`, optional + The measurement error at the conditioning points. + Either "nugget" to apply the model-nugget, a single value applied to + all points or an array with individual values for each point. + The measurement error has to be <= nugget. + The "exact=True" variant only works with "cond_err='nugget'". + Default: "nugget" + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: True + pseudo_inv_type : :class:`str` or :any:`callable`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `"pinv"`: use `pinv` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values + + If you want to use another routine to invert the kriging matrix, + you can pass a callable which takes a matrix and returns the inverse. + Default: `"pinv"` + fit_normalizer : :class:`bool`, optional + Whether to fit the data-normalizer to the given conditioning data. + Default: False + fit_variogram : :class:`bool`, optional + Whether to fit the given variogram model to the data. + Default: False + + References + ---------- + .. [Samson2020] Samson, M., & Deutsch, C. V. (2020). Collocated Cokriging. + In J. L. Deutsch (Ed.), Geostatistics Lessons. Retrieved from + http://geostatisticslessons.com/lessons/collocatedcokriging + .. [Wackernagel2003] Wackernagel, H. Multivariate Geostatistics, + Springer, Berlin, 2003. + """ + + def __init__( + self, + correlogram, + cond_pos, + cond_val, + secondary_cond_pos, + secondary_cond_val, + normalizer=None, + trend=None, + exact=False, + cond_err="nugget", + pseudo_inv=True, + pseudo_inv_type="pinv", + fit_normalizer=False, + fit_variogram=False, + ): + # Check if correlogram is actually a Correlogram object + if not isinstance(correlogram, Correlogram): + raise TypeError( + f"First argument must be a Correlogram instance. " + f"Got {type(correlogram).__name__}. " + ) + + # Initialize using base class with intrinsic algorithm + super().__init__( + correlogram=correlogram, + cond_pos=cond_pos, + cond_val=cond_val, + algorithm="intrinsic", + secondary_cond_pos=secondary_cond_pos, + secondary_cond_val=secondary_cond_val, + normalizer=normalizer, + trend=trend, + exact=exact, + cond_err=cond_err, + pseudo_inv=pseudo_inv, + pseudo_inv_type=pseudo_inv_type, + fit_normalizer=fit_normalizer, + fit_variogram=fit_variogram, + ) diff --git a/tests/test_cokriging.py b/tests/test_cokriging.py new file mode 100644 index 000000000..2af942e95 --- /dev/null +++ b/tests/test_cokriging.py @@ -0,0 +1,232 @@ +""" +This is the unittest of the cokriging module. +""" + +import unittest + +import numpy as np + +import gstools as gs + + +class TestCokriging(unittest.TestCase): + def setUp(self): + # Simple 1D test case + self.model = gs.Gaussian(dim=1, var=2, len_scale=2) + self.cond_pos = ([0.3, 1.9, 1.1, 3.3, 4.7],) + self.cond_val = np.array([0.47, 0.56, 0.74, 1.47, 1.74]) + self.sec_cond_val = np.array([1.8, 1.2, 2.1, 2.9, 2.4]) + self.pos = np.linspace(0, 5, 51) + # Dummy secondary data + self.sec_data = np.random.RandomState(42).rand(len(self.pos)) + + def test_secondary_data_required(self): + """Test that secondary_data is required on call.""" + correlogram = gs.MarkovModel1( + self.model, cross_corr=0.5, secondary_var=1.0 + ) + scck = gs.cokriging.SimpleCollocated( + correlogram, self.cond_pos, self.cond_val + ) + with self.assertRaises(ValueError): + scck(self.pos) + + def test_correlogram_type_required(self): + """Test that first argument must be a Correlogram.""" + with self.assertRaises(TypeError): + gs.cokriging.SimpleCollocated( + self.model, self.cond_pos, self.cond_val + ) + + def test_icck_secondary_cond_required(self): + """Test ICCK requires secondary conditioning data.""" + correlogram = gs.MarkovModel1( + self.model, cross_corr=0.5, secondary_var=1.0 + ) + with self.assertRaises(ValueError): + gs.cokriging.IntrinsicCollocated( + correlogram, + self.cond_pos, + self.cond_val, + secondary_cond_pos=None, + secondary_cond_val=None, + ) + + def test_icck_secondary_cond_length(self): + """Test ICCK secondary conditioning data length validation.""" + correlogram = gs.MarkovModel1( + self.model, cross_corr=0.5, secondary_var=1.0 + ) + with self.assertRaises(ValueError): + gs.cokriging.IntrinsicCollocated( + correlogram, + self.cond_pos, + self.cond_val, + self.cond_pos, + self.sec_cond_val[:3], # Wrong length + ) + + def test_zero_correlation_equals_sk(self): + """Test that ρ=0 gives Simple Kriging results.""" + # Reference: Simple Kriging + sk = gs.krige.Simple( + self.model, self.cond_pos, self.cond_val, mean=0.0 + ) + sk_field, sk_var = sk(self.pos, return_var=True) + + # SCCK with ρ=0 + correlogram_scck = gs.MarkovModel1( + self.model, cross_corr=0.0, secondary_var=1.5 + ) + scck = gs.cokriging.SimpleCollocated( + correlogram_scck, self.cond_pos, self.cond_val + ) + scck_field, scck_var = scck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + np.testing.assert_allclose(scck_field, sk_field, rtol=1e-6, atol=1e-9) + np.testing.assert_allclose(scck_var, sk_var, rtol=1e-6, atol=1e-9) + + # ICCK with ρ=0 + correlogram_icck = gs.MarkovModel1( + self.model, cross_corr=0.0, secondary_var=1.5 + ) + icck = gs.cokriging.IntrinsicCollocated( + correlogram_icck, + self.cond_pos, + self.cond_val, + self.cond_pos, + self.sec_cond_val, + ) + icck_field, icck_var = icck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + np.testing.assert_allclose(icck_field, sk_field, rtol=1e-6, atol=1e-9) + np.testing.assert_allclose(icck_var, sk_var, rtol=1e-6, atol=1e-9) + + def test_scck_variance_formula(self): + """Test SCCK variance formula.""" + cross_corr = 0.7 + secondary_var = 1.5 + + # Get SK variance + sk = gs.krige.Simple( + self.model, self.cond_pos, self.cond_val, mean=0.0 + ) + _, sk_var = sk(self.pos, return_var=True) + + # Calculate expected SCCK variance + C_Z0 = self.model.sill + C_Y0 = secondary_var + C_YZ0 = cross_corr * np.sqrt(C_Z0 * C_Y0) + k = C_YZ0 / C_Z0 + + numerator = k * sk_var + denominator = C_Y0 - (k**2) * (C_Z0 - sk_var) + lambda_Y0 = np.where( + np.abs(denominator) < 1e-15, 0.0, numerator / denominator + ) + expected_var = sk_var * (1.0 - lambda_Y0 * k) + expected_var = np.maximum(0.0, expected_var) + + # Actual SCCK variance + correlogram = gs.MarkovModel1( + self.model, cross_corr=cross_corr, secondary_var=secondary_var + ) + scck = gs.cokriging.SimpleCollocated( + correlogram, self.cond_pos, self.cond_val + ) + _, actual_var = scck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + np.testing.assert_allclose( + actual_var, expected_var, rtol=1e-6, atol=1e-9 + ) + + def test_icck_variance_formula(self): + """Test ICCK variance formula.""" + cross_corr = 0.7 + secondary_var = 1.5 + + # Get SK variance + sk = gs.krige.Simple( + self.model, self.cond_pos, self.cond_val, mean=0.0 + ) + _, sk_var = sk(self.pos, return_var=True) + + # Expected ICCK variance + C_Z0 = self.model.sill + C_Y0 = secondary_var + C_YZ0 = cross_corr * np.sqrt(C_Z0 * C_Y0) + rho_squared = (C_YZ0**2) / (C_Y0 * C_Z0) + expected_var = (1.0 - rho_squared) * sk_var + + # Actual ICCK variance + correlogram = gs.MarkovModel1( + self.model, cross_corr=cross_corr, secondary_var=secondary_var + ) + icck = gs.cokriging.IntrinsicCollocated( + correlogram, + self.cond_pos, + self.cond_val, + self.cond_pos, + self.sec_cond_val, + ) + _, actual_var = icck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + np.testing.assert_allclose( + actual_var, expected_var, rtol=1e-6, atol=1e-9 + ) + + def test_perfect_correlation_variance(self): + """Test that ρ=±1 gives near-zero variance for ICCK.""" + for rho in [-1.0, 1.0]: + correlogram = gs.MarkovModel1( + self.model, cross_corr=rho, secondary_var=1.5 + ) + icck = gs.cokriging.IntrinsicCollocated( + correlogram, + self.cond_pos, + self.cond_val, + self.cond_pos, + self.sec_cond_val, + ) + _, icck_var = icck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + self.assertTrue(np.allclose(icck_var, 0.0, atol=1e-12)) + + def test_variance_reduction(self): + """Test that cokriging reduces variance compared to simple kriging.""" + cross_corr = 0.8 + secondary_var = 1.5 + + # Get SK variance + sk = gs.krige.Simple( + self.model, self.cond_pos, self.cond_val, mean=0.0 + ) + _, sk_var = sk(self.pos, return_var=True) + + # Get ICCK variance + correlogram = gs.MarkovModel1( + self.model, cross_corr=cross_corr, secondary_var=secondary_var + ) + icck = gs.cokriging.IntrinsicCollocated( + correlogram, + self.cond_pos, + self.cond_val, + self.cond_pos, + self.sec_cond_val, + ) + _, icck_var = icck( + self.pos, secondary_data=self.sec_data, return_var=True + ) + + # ICCK variance ≤ SK variance + self.assertTrue(np.all(icck_var <= sk_var + 1e-8)) + self.assertTrue(np.mean(icck_var) < np.mean(sk_var)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_correlogram.py b/tests/test_correlogram.py new file mode 100644 index 000000000..a3a583aa0 --- /dev/null +++ b/tests/test_correlogram.py @@ -0,0 +1,116 @@ +""" +This is the unittest of the correlogram module. +""" + +import unittest + +import numpy as np + +import gstools as gs + + +class TestCorrelogram(unittest.TestCase): + def setUp(self): + self.model = gs.Gaussian(dim=1, var=0.5, len_scale=2.0) + self.cross_corr = 0.8 + self.secondary_var = 1.5 + self.primary_mean = 1.0 + self.secondary_mean = 0.5 + + def test_markov_model1_covariances(self): + """Test MM1 covariance computation at zero lag.""" + mm1 = gs.MarkovModel1( + primary_model=self.model, + cross_corr=self.cross_corr, + secondary_var=self.secondary_var, + primary_mean=self.primary_mean, + secondary_mean=self.secondary_mean, + ) + + C_Z0, C_Y0, C_YZ0 = mm1.compute_covariances() + + # Check primary variance + self.assertAlmostEqual(C_Z0, self.model.sill) + # Check secondary variance + self.assertAlmostEqual(C_Y0, self.secondary_var) + # Check cross-covariance formula + expected_C_YZ0 = self.cross_corr * np.sqrt(C_Z0 * C_Y0) + self.assertAlmostEqual(C_YZ0, expected_C_YZ0) + + def test_markov_model1_cross_covariance(self): + """Test MM1 cross-covariance formula at distance h.""" + mm1 = gs.MarkovModel1( + primary_model=self.model, + cross_corr=self.cross_corr, + secondary_var=self.secondary_var, + ) + + # Test at zero lag + C_YZ_0 = mm1.cross_covariance(0.0) + _, _, C_YZ0_expected = mm1.compute_covariances() + self.assertAlmostEqual(C_YZ_0, C_YZ0_expected) + + # Test MM1 formula: C_YZ(h) = (C_YZ(0) / C_Z(0)) * C_Z(h) + h = 1.0 + C_YZ_h = mm1.cross_covariance(h) + C_Z0, _, C_YZ0 = mm1.compute_covariances() + C_Z_h = self.model.covariance(h) + expected = (C_YZ0 / C_Z0) * C_Z_h + self.assertAlmostEqual(C_YZ_h, expected) + + def test_markov_model1_cross_covariance_array(self): + """Test MM1 cross-covariance with array input.""" + mm1 = gs.MarkovModel1( + primary_model=self.model, + cross_corr=self.cross_corr, + secondary_var=self.secondary_var, + ) + + h_array = np.array([0.0, 0.5, 1.0, 2.0]) + C_YZ_array = mm1.cross_covariance(h_array) + + # Check array shape + self.assertEqual(C_YZ_array.shape, h_array.shape) + + # Verify each element matches scalar computation + for i, h in enumerate(h_array): + C_YZ_single = mm1.cross_covariance(h) + self.assertAlmostEqual(C_YZ_array[i], C_YZ_single) + + def test_validation_cross_corr(self): + """Test parameter validation for cross_corr.""" + # cross_corr too large + with self.assertRaises(ValueError): + gs.MarkovModel1( + primary_model=self.model, + cross_corr=1.5, + secondary_var=self.secondary_var, + ) + # cross_corr too small + with self.assertRaises(ValueError): + gs.MarkovModel1( + primary_model=self.model, + cross_corr=-1.5, + secondary_var=self.secondary_var, + ) + + def test_validation_secondary_var(self): + """Test parameter validation for secondary_var.""" + # negative variance + with self.assertRaises(ValueError): + gs.MarkovModel1( + primary_model=self.model, + cross_corr=self.cross_corr, + secondary_var=-1.0, + ) + # zero variance + with self.assertRaises(ValueError): + gs.MarkovModel1( + primary_model=self.model, + cross_corr=self.cross_corr, + secondary_var=0.0, + ) + + +if __name__ == "__main__": + unittest.main()