Merge pull request #153 from nannau/feature/exact-values-option

Feature: exact values option

examples/exact_values_example_1D.py

@@ -0,0 +1,59 @@
+# -*- coding: utf-8 -*-
+"""
+Exact Values
+=================
+
+PyKrige demonstration and usage
+as a non-exact interpolator in 1D. 
+"""
+
+from pykrige.ok import OrdinaryKriging
+import matplotlib.pyplot as plt
+import numpy as np
+
+plt.style.use("ggplot")
+
+np.random.seed(42)
+
+x = np.linspace(0, 12.5, 50)
+xpred = np.linspace(0, 12.5, 393)
+y = np.sin(x) * np.exp(-0.25 * x) + np.random.normal(-0.25, 0.25, 50)
+
+# compare OrdinaryKriging as an exact and non exact interpolator
+uk = OrdinaryKriging(
+    x, np.zeros(x.shape), y, variogram_model="linear", exact_values=False
+)
+uk_exact = OrdinaryKriging(x, np.zeros(x.shape), y, variogram_model="linear")
+
+y_pred, y_std = uk.execute("grid", xpred, np.array([0.0]), backend="loop")
+y_pred_exact, y_std_exact = uk_exact.execute(
+    "grid", xpred, np.array([0.0]), backend="loop"
+)
+
+
+y_pred = np.squeeze(y_pred)
+y_std = np.squeeze(y_std)
+
+y_pred_exact = np.squeeze(y_pred_exact)
+y_std_exact = np.squeeze(y_std_exact)
+
+
+fig, ax = plt.subplots(1, 1, figsize=(10, 4))
+
+ax.scatter(x, y, label="Input Data")
+ax.plot(xpred, y_pred_exact, label="Exact Prediction")
+ax.plot(xpred, y_pred, label="Non Exact Prediction")
+
+ax.fill_between(
+    xpred,
+    y_pred - 3 * y_std,
+    y_pred + 3 * y_std,
+    alpha=0.3,
+    label="Confidence interval",
+)
+ax.legend(loc=9)
+ax.set_ylim(-1.8, 1.3)
+ax.legend(loc=9)
+plt.xlabel("X")
+plt.ylabel("Field")
+plt.show()

pykrige/lib/cok.pyx

@@ -48,7 +48,8 @@ cpdef _c_exec_loop(double [:, ::1] a_all,
             b[k] *= -1
         b[n] = 1.0
 
-        check_b_vect(n, bd, b, eps)
+        if not pars['exact_values']:
+            check_b_vect(n, bd, b, eps)
 
 
         # Do the BLAS matrix-vector multiplication call

pykrige/ok.py

@@ -15,6 +15,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
+.. [2] N. Cressie, Statistics for spatial data, 
+   (Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
 """
@@ -140,11 +142,21 @@ class OrdinaryKriging:
         coordinates, x is interpreted as longitude and y as latitude
         coordinates, both given in degree. Longitudes are expected in
         [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'.
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     eps = 1.0e-10  # Cutoff for comparison to zero
@@ -173,11 +185,17 @@ def __init__(
         enable_plotting=False,
         enable_statistics=False,
         coordinates_type="euclidean",
+        exact_values=True,
     ):
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -585,11 +603,11 @@ def _get_kriging_matrix(self, n):
             )
         a = np.zeros((n + 1, n + 1))
         a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d)
+
         np.fill_diagonal(a, 0.0)
         a[n, :] = 1.0
         a[:, n] = 1.0
         a[n, n] = 0.0
-
         return a
 
     def _exec_vector(self, a, bd, mask):
@@ -609,7 +627,7 @@ def _exec_vector(self, a, bd, mask):
 
         b = np.zeros((npt, n + 1, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
         b[:, n, 0] = 1.0
 
@@ -647,7 +665,7 @@ def _exec_loop(self, a, bd_all, mask):
 
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
             x = np.dot(a_inv, b)
@@ -683,7 +701,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx):
                 zero_value = False
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 
@@ -705,25 +723,6 @@ def execute(
     ):
         """Calculates a kriged grid and the associated variance.
 
-        This is now the method that performs the main kriging calculation.
-        Note that currently measurements (i.e., z values) are considered
-        'exact'. This means that, when a specified coordinate for interpolation
-        is exactly the same as one of the data points, the variogram evaluated
-        at the point is forced to be zero. Also, the diagonal of the kriging
-        matrix is also always forced to be zero. In forcing the variogram
-        evaluated at data points to be zero, we are effectively saying that
-        there is no variance at that point (no uncertainty,
-        so the value is 'exact').
-
-        In the future, the code may include an extra 'exact_values' boolean
-        flag that can be adjusted to specify whether to treat the measurements
-        as 'exact'. Setting the flag to false would indicate that the
-        variogram should not be forced to be zero at zero distance
-        (i.e., when evaluated at data points). Instead, the uncertainty in
-        the point will be equal to the nugget. This would mean that the
-        diagonal of the kriging matrix would be set to
-        the nugget instead of to zero.
-
         Parameters
         ----------
         style : str
@@ -876,6 +875,7 @@ def execute(
                     "eps",
                     "variogram_model_parameters",
                     "variogram_function",
+                    "exact_values",
                 ]
             }
 

pykrige/ok3d.py

@@ -14,6 +14,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
+.. [2] N. Cressie, Statistics for spatial data, 
+(Wiley Series in Probability and Statistics, 1993) 137 p.
 
 Copyright (c) 2015-2020, PyKrige Developers
 """
@@ -151,11 +153,21 @@ class OrdinaryKriging3D:
         Default is False (off).
     enable_plotting : bool, optional
         Enables plotting to display variogram. Default is False (off).
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     eps = 1.0e-10  # Cutoff for comparison to zero
@@ -186,11 +198,17 @@ def __init__(
         anisotropy_angle_z=0.0,
         verbose=False,
         enable_plotting=False,
+        exact_values=True,
     ):
 
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -594,7 +612,7 @@ def _exec_vector(self, a, bd, mask):
 
         b = np.zeros((npt, n + 1, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
         b[:, n, 0] = 1.0
 
@@ -632,7 +650,7 @@ def _exec_loop(self, a, bd_all, mask):
 
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 
@@ -669,7 +687,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx):
                 zero_index = None
             b = np.zeros((n + 1, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
             b[n, 0] = 1.0
 

pykrige/uk.py

@@ -15,7 +15,8 @@
 ----------
 .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
     Hydrogeology, (Cambridge University Press, 1997) 272 p.
-
+.. [2] N. Cressie, Statistics for spatial data, 
+   (Wiley Series in Probability and Statistics, 1993) 137 p.
 Copyright (c) 2015-2020, PyKrige Developers
 """
 import numpy as np
@@ -172,11 +173,21 @@ class UniversalKriging:
         Default is False (off).
     enable_plotting : boolean, optional
         Enables plotting to display variogram. Default is False (off).
+    exact_values : bool, optional
+        If True, interpolation provides input values at input locations. 
+        If False, interpolation accounts for variance/nugget within input 
+        values at input locations and does not behave as an 
+        exact-interpolator [2]. Note that this only has an effect if
+        there is variance/nugget present within the input data since it is
+        interpreted as measurement error. If the nugget is zero, the kriged
+        field will behave as an exact interpolator.
 
     References
     ----------
     .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in
        Hydrogeology, (Cambridge University Press, 1997) 272 p.
+    .. [2] N. Cressie, Statistics for spatial data, 
+       (Wiley Series in Probability and Statistics, 1993) 137 p.
     """
 
     UNBIAS = True  # This can be changed to remove the unbiasedness condition
@@ -212,6 +223,7 @@ def __init__(
         functional_drift=None,
         verbose=False,
         enable_plotting=False,
+        exact_values=True,
     ):
 
         # Deal with mutable default argument
@@ -225,6 +237,11 @@ def __init__(
         # set up variogram model and parameters...
         self.variogram_model = variogram_model
         self.model = None
+
+        if not isinstance(exact_values, bool):
+            raise ValueError("exact_values has to be boolean True or False")
+        self.exact_values = exact_values
+
         # check if a GSTools covariance model is given
         if hasattr(self.variogram_model, "pykrige_kwargs"):
             # save the model in the class
@@ -820,6 +837,7 @@ def _get_kriging_matrix(self, n, n_withdrifts):
         else:
             a = np.zeros((n_withdrifts, n_withdrifts))
         a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d)
+
         np.fill_diagonal(a, 0.0)
 
         i = n
@@ -888,7 +906,7 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids)
         else:
             b = np.zeros((npt, n_withdrifts, 1))
         b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-        if zero_value:
+        if zero_value and self.exact_values:
             b[zero_index[0], zero_index[1], 0] = 0.0
 
         i = n
@@ -979,7 +997,7 @@ def _exec_loop(self, a, bd_all, xy, xy_orig, mask, n_withdrifts, spec_drift_grid
             else:
                 b = np.zeros((n_withdrifts, 1))
             b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd)
-            if zero_value:
+            if zero_value and self.exact_values:
                 b[zero_index[0], 0] = 0.0
 
             i = n
@@ -1039,24 +1057,6 @@ def execute(
         """Calculates a kriged grid and the associated variance.
         Includes drift terms.
 
-        This is now the method that performs the main kriging calculation.
-        Note that currently measurements (i.e., z values) are considered
-        'exact'. This means that, when a specified coordinate for interpolation
-        is exactly the same as one of the data points, the variogram evaluated
-        at the point is forced to be zero. Also, the diagonal of the kriging
-        matrix is also always forced to be zero. In forcing the variogram
-        evaluated at data points to be zero, we are effectively saying that
-        there is no variance at that point (no uncertainty,
-        so the value is 'exact').
-
-        In the future, the code may include an extra 'exact_values' boolean
-        flag that can be adjusted to specify whether to treat the measurements
-        as 'exact'. Setting the flag to false would indicate that the variogram
-        should not be forced to be zero at zero distance (i.e., when evaluated
-        at data points). Instead, the uncertainty in the point will be equal to
-        the nugget. This would mean that the diagonal of the kriging matrix
-        would be set to the nugget instead of to zero.
-
         Parameters
         ----------
         style : str