update kriging routines

gstools/krige/krigesum.pyx

@@ -21,21 +21,18 @@ def krigesum(double[:,:] krig_mat, double[:,:] krig_vecs, double[:] cond):
 
     cdef double[:] field = np.zeros(res_i)
     cdef double[:] error = np.zeros(res_i)
-    cdef double[:] krig_facs = np.zeros(mat_i)
+    cdef double krig_fac
 
     cdef int i, j, k
 
-    # krig_facs = cond * krig_mat
-    for j in prange(mat_i, nogil=True):
-        for i in range(mat_i):
-            krig_facs[j] += cond[i] * krig_mat[i,j]
-
     # error = krig_vecs * krig_mat * krig_vecs
     # field = krig_facs * krig_vecs
     for k in prange(res_i, nogil=True):
         for i in range(mat_i):
+            krig_fac = 0.0
             for j in range(mat_i):
-                error[k] += krig_mat[i,j] * krig_vecs[i,k] * krig_vecs[j,k]
-            field[k] += krig_facs[i] * krig_vecs[i,k]
+                krig_fac += krig_mat[i,j] * krig_vecs[j,k]
+            error[k] += krig_vecs[i,k] * krig_fac
+            field[k] += cond[i] * krig_fac
 
     return np.asarray(field), np.asarray(error)

gstools/krige/ordinary.py

@@ -33,7 +33,8 @@
 
 
 class Ordinary(object):
-    """A class for ordinary kriging.
+    """
+    A class for ordinary kriging.
 
     Parameters
     ----------
@@ -47,18 +48,21 @@ class Ordinary(object):
 
     def __init__(self, model, cond_pos, cond_val):
         # initialize private attributes
+        self.field = None
+        self.error = None
+
         self._model = None
         self._cond_pos = None
         self._cond_val = None
 
-        self.field = None
         self.model = model
         self.set_condition(cond_pos, cond_val)
 
         # initialize attributes
 
     def __call__(self, pos, mesh_type="unstructured"):
-        """Generate the ordinary kriging field.
+        """
+        Generate the ordinary kriging field.
 
         The field is saved as `self.field` and is also returned.
 
@@ -78,8 +82,8 @@ def __call__(self, pos, mesh_type="unstructured"):
             the kriging error
         """
         # internal conversation
-        x, y, z = pos2xyz(pos)
-        c_x, c_y, c_z = pos2xyz(self.cond_pos)
+        x, y, z = pos2xyz(pos, dtype=np.double)
+        c_x, c_y, c_z = pos2xyz(self.cond_pos, dtype=np.double)
         # format the positional arguments of the mesh
         check_mesh(self.model.dim, x, y, z, mesh_type)
         mesh_type_changed = False
@@ -99,17 +103,9 @@ def __call__(self, pos, mesh_type="unstructured"):
         c_y, c_z = make_isotropic(self.model.dim, self.model.anis, c_y, c_z)
 
         # set condtions
-        cond = np.array(np.concatenate((self.cond_val, [0])), dtype=np.double)
-        krig_mat = np.array(
-            inv(
-                self._get_krig_mat((c_x, c_y, c_z), (c_x, c_y, c_z), c_x.size)
-            ),
-            dtype=np.double,
-        )
-        krig_vecs = np.array(
-            self._get_vario_mat((c_x, c_y, c_z), (x, y, z), add=True),
-            dtype=np.double,
-        )
+        cond = np.concatenate((self.cond_val, [0]))
+        krig_mat = inv(self._get_krig_mat((c_x, c_y, c_z), (c_x, c_y, c_z)))
+        krig_vecs = self._get_vario_mat((c_x, c_y, c_z), (x, y, z), add=True)
         # generate the kriged field
         field, error = krigesum(krig_mat, krig_vecs, cond)
 
@@ -127,8 +123,9 @@ def __call__(self, pos, mesh_type="unstructured"):
         self.field = field
         return self.field, self.error
 
-    def _get_krig_mat(self, pos1, pos2, size):
-        res = np.empty((size + 1, size + 1), dtype=float)
+    def _get_krig_mat(self, pos1, pos2):
+        size = pos1[0].size
+        res = np.empty((size + 1, size + 1), dtype=np.double)
         res[:size, :size] = self._get_vario_mat(pos1, pos2)
         res[size, :] = 1
         res[:, size] = 1
@@ -157,18 +154,18 @@ def set_condition(self, cond_pos, cond_val):
             the values of the conditions
         """
         self._cond_pos = cond_pos
-        self._cond_val = cond_val
+        self._cond_val = np.array(cond_val, dtype=np.double)
 
     def structured(self, *args, **kwargs):
-        """Ordinary kriging on a structured mesh
+        """Ordinary kriging on a structured mesh.
 
         See :any:`Ordinary.__call__`
         """
         call = partial(self.__call__, mesh_type="structured")
         return call(*args, **kwargs)
 
     def unstructured(self, *args, **kwargs):
-        """Ordinary kriging on an unstructured mesh
+        """Ordinary kriging on an unstructured mesh.
 
         See :any:`Ordinary.__call__`
         """
@@ -177,20 +174,17 @@ def unstructured(self, *args, **kwargs):
 
     @property
     def cond_pos(self):
-        """:class:`list`: The position tuple of the conditions.
-        """
+        """:class:`list`: The position tuple of the conditions."""
         return self._cond_pos
 
     @property
     def cond_val(self):
-        """:class:`list`: The values of the conditions.
-        """
+        """:class:`list`: The values of the conditions."""
         return self._cond_val
 
     @property
     def model(self):
-        """:any:`CovModel`: The covariance model used for kriging.
-        """
+        """:any:`CovModel`: The covariance model used for kriging."""
         return self._model
 
     @model.setter
@@ -205,15 +199,15 @@ def model(self, model):
 
     @property
     def do_rotation(self):
-        """:any:`bool`: State if a rotation should be performed
-        depending on the model.
-        """
+        """:any:`bool`: State if a rotation should be performed."""
         return not np.all(np.isclose(self.model.angles, 0.0))
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "Ordinary(model={0}, cond_pos={1}, cond_val={2}".format(
             self.model, self.cond_pos, self.cond_val
         )

gstools/krige/simple.py

@@ -33,7 +33,8 @@
 
 
 class Simple(object):
-    """A class for simple kriging.
+    """
+    A class for simple kriging.
 
     Parameters
     ----------
@@ -63,7 +64,8 @@ def __init__(self, model, mean, cond_pos, cond_val):
         # initialize attributes
 
     def __call__(self, pos, mesh_type="unstructured"):
-        """Generate the simple kriging field.
+        """
+        Generate the simple kriging field.
 
         The field is saved as `self.field` and is also returned.
 
@@ -83,8 +85,8 @@ def __call__(self, pos, mesh_type="unstructured"):
             the kriging error
         """
         # internal conversation
-        x, y, z = pos2xyz(pos)
-        c_x, c_y, c_z = pos2xyz(self.cond_pos)
+        x, y, z = pos2xyz(pos, dtype=np.double)
+        c_x, c_y, c_z = pos2xyz(self.cond_pos, dtype=np.double)
         # format the positional arguments of the mesh
         check_mesh(self.model.dim, x, y, z, mesh_type)
         mesh_type_changed = False
@@ -104,14 +106,9 @@ def __call__(self, pos, mesh_type="unstructured"):
         c_y, c_z = make_isotropic(self.model.dim, self.model.anis, c_y, c_z)
 
         # set condtions to zero mean
-        cond = np.array(self.cond_val - self.mean, dtype=np.double)
-        krig_mat = np.array(
-            inv(self._get_cov_mat((c_x, c_y, c_z), (c_x, c_y, c_z))),
-            dtype=np.double,
-        )
-        krig_vecs = np.array(
-            self._get_cov_mat((c_x, c_y, c_z), (x, y, z)), dtype=np.double
-        )
+        cond = self.cond_val - self.mean
+        krig_mat = inv(self._get_cov_mat((c_x, c_y, c_z), (c_x, c_y, c_z)))
+        krig_vecs = self._get_cov_mat((c_x, c_y, c_z), (x, y, z))
         # generate the kriged field
         field, error = krigesum(krig_mat, krig_vecs, cond)
 
@@ -152,15 +149,15 @@ def set_condition(self, cond_pos, cond_val):
         self._cond_val = np.array(cond_val, dtype=np.double)
 
     def structured(self, *args, **kwargs):
-        """Simple kriging on a structured mesh
+        """Simple kriging on a structured mesh.
 
         See :any:`Simple.__call__`
         """
         call = partial(self.__call__, mesh_type="structured")
         return call(*args, **kwargs)
 
     def unstructured(self, *args, **kwargs):
-        """Simple kriging on an unstructured mesh
+        """Simple kriging on an unstructured mesh.
 
         See :any:`Simple.__call__`
         """
@@ -169,20 +166,17 @@ def unstructured(self, *args, **kwargs):
 
     @property
     def cond_pos(self):
-        """:class:`list`: The position tuple of the conditions.
-        """
+        """:class:`list`: The position tuple of the conditions."""
         return self._cond_pos
 
     @property
     def cond_val(self):
-        """:class:`list`: The values of the conditions.
-        """
+        """:class:`list`: The values of the conditions."""
         return self._cond_val
 
     @property
     def model(self):
-        """:any:`CovModel`: The covariance model used for kriging.
-        """
+        """:any:`CovModel`: The covariance model used for kriging."""
         return self._model
 
     @model.setter
@@ -198,15 +192,15 @@ def model(self, model):
 
     @property
     def do_rotation(self):
-        """:any:`bool`: State if a rotation should be performed
-        depending on the model.
-        """
+        """:any:`bool`: State if a rotation should be performed."""
         return not np.all(np.isclose(self.model.angles, 0.0))
 
     def __str__(self):
+        """Return String representation."""
         return self.__repr__()
 
     def __repr__(self):
+        """Return String representation."""
         return "Simple(model={0}, mean={1}, cond_pos={2}, cond_val={3}".format(
             self.model, self.mean, self.cond_pos, self.cond_val
         )
@@ -216,3 +210,11 @@ def __repr__(self):
     import doctest
 
     doctest.testmod()
+    # from gstools import Gaussian
+    # gridx = gridy = np.linspace(0, 10, 1000)
+    # data = np.random.random((100, 3))
+    # model = Gaussian(dim=2, len_scale=1.5, anis=0.2, angles=0.1, var=0.5)
+    # krige = Simple(
+    #     model, mean=1, cond_pos=(data[:, 0], data[:, 1]), cond_val=data[:, 2]
+    # )
+    # krige((gridx, gridy), mesh_type="structured")