From 6020e1e1a7fd0d444573ea28399d811b5a219d6f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 00:34:09 +0100 Subject: [PATCH 01/28] Generator: remove structured option from generators --- gstools/field/generator.py | 125 ++++---------------------- gstools/field/summator.pyx | 180 +------------------------------------ 2 files changed, 18 insertions(+), 287 deletions(-) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index a5a9c2ab..c0a48862 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -16,12 +16,8 @@ import numpy as np from gstools.covmodel.base import CovModel from gstools.random.rng import RNG -from gstools.field.summator import ( - summate_unstruct, - summate_struct, - summate_incompr_unstruct, - summate_incompr_struct, -) +from gstools.field.summator import summate_unstruct, summate_incompr_unstruct + __all__ = ["RandMeth", "IncomprRandMeth"] @@ -85,7 +81,7 @@ def __init__( # set model and seed self.update(model, seed) - def __call__(self, x, y=None, z=None, mesh_type="unstructured"): + def __call__(self, pos): """Calculate the random modes for the randomization method. This method calls the `summate_*` Cython methods, which are the @@ -93,32 +89,18 @@ def __call__(self, x, y=None, z=None, mesh_type="unstructured"): Parameters ---------- - x : :class:`float`, :class:`numpy.ndarray` - The x components of the pos. tuple. - y : :class:`float`, :class:`numpy.ndarray`, optional - The y components of the pos. tuple. - z : :class:`float`, :class:`numpy.ndarray`, optional - The z components of the pos. tuple. - mesh_type : :class:`str`, optional - 'structured' / 'unstructured' + pos : (d, n), :class:`numpy.ndarray` + the position tuple with d dimensions and n points. Returns ------- :class:`numpy.ndarray` the random modes """ - if mesh_type == "unstructured": - pos = _reshape_pos(x, y, z, dtype=np.double) - - summed_modes = summate_unstruct( - self._cov_sample, self._z_1, self._z_2, pos - ) - else: - x, y, z = _set_dtype(x, y, z, dtype=np.double) - summed_modes = summate_struct( - self._cov_sample, self._z_1, self._z_2, x, y, z - ) - + pos = np.array(pos, dtype=np.double) + summed_modes = summate_unstruct( + self._cov_sample, self._z_1, self._z_2, pos + ) nugget = self._set_nugget(summed_modes.shape) return np.sqrt(self.model.var / self._mode_no) * summed_modes + nugget @@ -367,7 +349,7 @@ def __init__( self.mean_u = mean_velocity self._value_type = "vector" - def __call__(self, x, y=None, z=None, mesh_type="unstructured"): + def __call__(self, pos): """Calculate the random modes for the randomization method. This method calls the `summate_incompr_*` Cython methods, @@ -377,34 +359,18 @@ def __call__(self, x, y=None, z=None, mesh_type="unstructured"): Parameters ---------- - x : :class:`float`, :class:`numpy.ndarray` - the x components of the position tuple, the shape has to be - (len(x), 1, 1) for 3d and accordingly shorter for lower - dimensions - y : :class:`float`, :class:`numpy.ndarray`, optional - the y components of the pos. tuples. Default: ``None`` - z : :class:`float`, :class:`numpy.ndarray`, optional - the z components of the pos. tuple. Default: ``None`` - mesh_type : :class:`str`, optional - 'structured' / 'unstructured' + pos : (d, n), :class:`numpy.ndarray` + the position tuple with d dimensions and n points. Returns ------- :class:`numpy.ndarray` the random modes """ - if mesh_type == "unstructured": - pos = _reshape_pos(x, y, z, dtype=np.double) - - summed_modes = summate_incompr_unstruct( - self._cov_sample, self._z_1, self._z_2, pos - ) - else: - x, y, z = _set_dtype(x, y, z, dtype=np.double) - summed_modes = summate_incompr_struct( - self._cov_sample, self._z_1, self._z_2, x, y, z - ) - + pos = np.array(pos, dtype=np.double) + summed_modes = summate_incompr_unstruct( + self._cov_sample, self._z_1, self._z_2, pos + ) nugget = self._set_nugget(summed_modes.shape) e1 = self._create_unit_vector(summed_modes.shape) @@ -441,62 +407,3 @@ def _create_unit_vector(self, broadcast_shape, axis=0): e1 = np.zeros(shape) e1[axis] = 1.0 return e1 - - -def _reshape_pos(x, y=None, z=None, dtype=np.double): - """ - Reshape the 1d x, y, z positions to a 2d position array. - - Parameters - ---------- - x : :class:`float`, :class:`numpy.ndarray` - the x components of the position tuple, the shape has to be - (len(x), 1, 1) for 3d and accordingly shorter for lower - dimensions - y : :class:`float`, :class:`numpy.ndarray`, optional - the y components of the pos. tuple - z : :class:`float`, :class:`numpy.ndarray`, optional - the z components of the pos. tuple - dtype : :class:`numpy.dtype`, optional - the numpy dtype to which the elements should be converted - - Returns - ------- - :class:`numpy.ndarray` - the positions in one convinient data structure - """ - if y is None and z is None: - pos = np.array(x.reshape(1, len(x)), dtype=dtype) - elif z is None: - pos = np.array(np.vstack((x, y)), dtype=dtype) - else: - pos = np.array(np.vstack((x, y, z)), dtype=dtype) - return pos - - -def _set_dtype(x, y=None, z=None, dtype=np.double): - """ - Convert the dtypes of the input arrays to given dtype. - - Parameters - ---------- - x : :class:`float`, :class:`numpy.ndarray` - The array to be converted. - y : :class:`float`, :class:`numpy.ndarray`, optional - The array to be converted. - z : :class:`float`, :class:`numpy.ndarray`, optional - The array to be converted. - dtype : :class:`numpy.dtype`, optional - The numpy dtype to which the elements should be converted. - - Returns - ------- - :class:`numpy.ndarray` - The input lists/ arrays as numpy arrays with given dtype. - """ - x = x.astype(dtype, copy=False) - if y is not None: - y = y.astype(dtype, copy=False) - if z is not None: - z = z.astype(dtype, copy=False) - return x, y, z diff --git a/gstools/field/summator.pyx b/gstools/field/summator.pyx index 490c7e63..c3da4459 100644 --- a/gstools/field/summator.pyx +++ b/gstools/field/summator.pyx @@ -1,7 +1,7 @@ #cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True # -*- coding: utf-8 -*- """ -This is the variogram estimater, implemented in cython. +This is the randomization method summator, implemented in cython. """ import numpy as np @@ -41,97 +41,6 @@ def summate_unstruct( return np.asarray(summed_modes) -def summate_struct( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y=None, - const double[:] z=None, -): - if y == None and z == None: - return summate_struct_1d(cov_samples, z_1, z_2, x) - elif z == None: - return summate_struct_2d(cov_samples, z_1, z_2, x, y) - else: - return summate_struct_3d(cov_samples, z_1, z_2, x, y, z) - -def summate_struct_1d( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - ): - - cdef int i, j, X_len, N - cdef double phase - - X_len = x.shape[0] - N = cov_samples.shape[1] - - cdef double[:] summed_modes = np.zeros(X_len, dtype=DTYPE) - - for i in prange(X_len, nogil=True): - for j in range(N): - phase = cov_samples[0,j] * x[i] - summed_modes[i] += z_1[j] * cos(phase) + z_2[j] * sin(phase) - - return np.asarray(summed_modes) - -def summate_struct_2d( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y, - ): - cdef int i, j, k, X_len, Y_len, N - cdef double phase - - X_len = x.shape[0] - Y_len = y.shape[0] - N = cov_samples.shape[1] - - cdef double[:,:] summed_modes = np.zeros((X_len, Y_len), dtype=DTYPE) - - for i in prange(X_len, nogil=True): - for j in range(Y_len): - for k in range(N): - phase = cov_samples[0,k] * x[i] + cov_samples[1,k] * y[j] - summed_modes[i,j] += z_1[k] * cos(phase) + z_2[k] * sin(phase) - - return np.asarray(summed_modes) - -def summate_struct_3d( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y, - const double[:] z, - ): - cdef int i, j, k, l, X_len, Y_len, Z_len, N - cdef double phase - - X_len = x.shape[0] - Y_len = y.shape[0] - Z_len = z.shape[0] - N = cov_samples.shape[1] - - cdef double[:,:,:] summed_modes = np.zeros((X_len, Y_len, Z_len), dtype=DTYPE) - - for i in prange(X_len, nogil=True): - for j in range(Y_len): - for k in range(Z_len): - for l in range(N): - phase = ( - cov_samples[0,l] * x[i] + - cov_samples[1,l] * y[j] + - cov_samples[2,l] * z[k] - ) - summed_modes[i,j,k] += z_1[l] * cos(phase) + z_2[l] * sin(phase) - - return np.asarray(summed_modes) cdef (double) abs_square(const double[:] vec) nogil: cdef int i @@ -142,6 +51,7 @@ cdef (double) abs_square(const double[:] vec) nogil: return r + def summate_incompr_unstruct( const double[:,:] cov_samples, const double[:] z_1, @@ -174,89 +84,3 @@ def summate_incompr_unstruct( summed_modes[d,i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase)) return np.asarray(summed_modes) - -def summate_incompr_struct( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y=None, - const double[:] z=None, -): - if z == None: - return summate_incompr_struct_2d(cov_samples, z_1, z_2, x, y) - else: - return summate_incompr_struct_3d(cov_samples, z_1, z_2, x, y, z) - -def summate_incompr_struct_2d( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y, - ): - cdef int i, j, k, d, X_len, Y_len, N - cdef double phase - cdef int dim = 2 - cdef double k_2 - - cdef double[:] e1 = np.zeros(dim, dtype=DTYPE) - e1[0] = 1. - cdef double[:] proj = np.empty(dim, dtype=DTYPE) - - X_len = x.shape[0] - Y_len = y.shape[0] - N = cov_samples.shape[1] - - cdef double[:,:,:] summed_modes = np.zeros((dim, X_len, Y_len), dtype=DTYPE) - - for i in range(X_len): - for j in range(Y_len): - for k in range(N): - k_2 = abs_square(cov_samples[:,k]) - phase = cov_samples[0,k] * x[i] + cov_samples[1,k] * y[j] - for d in range(dim): - proj[d] = e1[d] - cov_samples[d,k] * cov_samples[0,k] / k_2 - summed_modes[d,i,j] += proj[d] * (z_1[k] * cos(phase) + z_2[k] * sin(phase)) - - return np.asarray(summed_modes) - -def summate_incompr_struct_3d( - const double[:,:] cov_samples, - const double[:] z_1, - const double[:] z_2, - const double[:] x, - const double[:] y, - const double[:] z, - ): - cdef int i, j, k, l, d, X_len, Y_len, Z_len, N - cdef double phase - cdef int dim = 3 - cdef double k_2 - - cdef double[:] e1 = np.zeros(dim, dtype=DTYPE) - e1[0] = 1. - cdef double[:] proj = np.empty(dim, dtype=DTYPE) - - X_len = x.shape[0] - Y_len = y.shape[0] - Z_len = z.shape[0] - N = cov_samples.shape[1] - - cdef double[:,:,:,:] summed_modes = np.zeros((dim, X_len, Y_len, Z_len), dtype=DTYPE) - - for i in range(X_len): - for j in range(Y_len): - for k in range(Z_len): - for l in range(N): - k_2 = abs_square(cov_samples[:,l]) - phase = ( - cov_samples[0,l] * x[i] + - cov_samples[1,l] * y[j] + - cov_samples[2,l] * z[k] - ) - for d in range(dim): - proj[d] = e1[d] - cov_samples[d,l] * cov_samples[0,l] / k_2 - summed_modes[d,i,j,k] += proj[d] * (z_1[l] * cos(phase) + z_2[l] * sin(phase)) - - return np.asarray(summed_modes) From 3eb6c5813be80ae88db13d4648ef0ecb531abb61 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 00:42:13 +0100 Subject: [PATCH 02/28] Generator: explicitly assert dim=2;3 when generating vector fields --- gstools/field/generator.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index c0a48862..77858f68 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -339,10 +339,9 @@ def __init__( verbose=False, **kwargs ): - if model.dim < 2: + if model.dim < 2 or model.dim > 3: raise ValueError( - "Only 2- and 3-dimensional incompressible fields " - + "can be generated." + "Only 2D and 3D incompressible fields can be generated." ) super().__init__(model, mode_no, seed, verbose, **kwargs) From 74b12f72d53ffe02f46d4a4b0061e73cd972de37 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 01:27:27 +0100 Subject: [PATCH 03/28] Field: simplify pre_pos routine; use new CovModel.isometrize routine there --- gstools/field/base.py | 108 +++++++++-------- gstools/field/srf.py | 14 +-- gstools/field/tools.py | 256 ----------------------------------------- 3 files changed, 66 insertions(+), 312 deletions(-) delete mode 100644 gstools/field/tools.py diff --git a/gstools/field/base.py b/gstools/field/base.py index c3349e0a..182dc952 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -17,14 +17,7 @@ from gstools.covmodel.base import CovModel from gstools.tools.export import to_vtk, vtk_export -from gstools.field.tools import ( - _get_select, - check_mesh, - make_isotropic, - unrotate_mesh, - reshape_axis_from_struct_to_unstruct, -) -from gstools.tools.geometric import pos2xyz, xyz2pos + __all__ = ["Field"] @@ -154,7 +147,7 @@ def mesh( mesh.point_data[name] = field return out - def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False): + def pre_pos(self, pos, mesh_type="unstructured"): """ Preprocessing positions and mesh_type. @@ -165,46 +158,39 @@ def _pre_pos(self, pos, mesh_type="unstructured", make_unstruct=False): directions mesh_type : :class:`str` 'structured' / 'unstructured' - make_unstruct: :class:`bool` - State if mesh_type should be made unstructured. Returns ------- - x : :class:`numpy.ndarray` - first components of unrotated and isotropic position vectors - y : :class:`numpy.ndarray` or None - analog to x - z : :class:`numpy.ndarray` or None - analog to x - pos : :class:`tuple` of :class:`numpy.ndarray` - the normalized position tuple - mesh_type_gen : :class:`str` - 'structured' / 'unstructured' for the generator - mesh_type_changed : :class:`bool` - State if the mesh_type was changed. - axis_lens : :class:`tuple` or :any:`None` - axis lengths of the structured mesh if mesh type was changed. + iso_pos : (d, n), :class:`numpy.ndarray` + the isometrize position tuple + shape : :class:`tuple` + Shape of the resulting field. """ - x, y, z = pos2xyz(pos, dtype=np.double, max_dim=self.model.dim) - pos = xyz2pos(x, y, z) - mesh_type_gen = mesh_type - # format the positional arguments of the mesh - check_mesh(self.model.dim, x, y, z, mesh_type) - mesh_type_changed = False - axis_lens = None - if ( - self.model.do_rotation or make_unstruct - ) and mesh_type == "structured": - mesh_type_changed = True - mesh_type_gen = "unstructured" - x, y, z, axis_lens = reshape_axis_from_struct_to_unstruct( - self.model.dim, x, y, z - ) - if self.model.do_rotation: - x, y, z = unrotate_mesh(self.model.dim, self.model.angles, x, y, z) - if not self.model.is_isotropic: - y, z = make_isotropic(self.model.dim, self.model.anis, y, z) - return x, y, z, pos, mesh_type_gen, mesh_type_changed, axis_lens + # save mesh-type + self.mesh_type = mesh_type + # save pos tuple + if mesh_type != "unstructured": + if self.model.dim == 1: + self.pos = (np.array(pos, dtype=np.double).reshape(-1),) + else: + if len(pos) != self.model.dim: + raise ValueError( + "Field: position tuple doesn't match dimension." + ) + self.pos = tuple( + [ + np.array(pos_i, dtype=np.double).reshape(-1) + for pos_i in pos + ] + ) + pos = np.meshgrid(*self.pos, indexing="ij") + shape = tuple([len(pos_i) for pos_i in self.pos]) + else: + pos = np.array(pos, dtype=np.double).reshape(self.model.dim, -1) + self.pos = pos + shape = np.shape(self.pos[0]) + # return isometrized pos tuple and resulting field shape + return self.model.isometrize(pos), shape def _to_vtk_helper( self, filename=None, field_select="field", fieldname="field" @@ -396,3 +382,35 @@ def __repr__(self): return "Field(model={0}, mean={1:.{p}})".format( self.model, self.mean, p=self.model._prec ) + + +def _get_select(direction): + select = [] + if not (0 < len(direction) < 4): + raise ValueError( + "Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction) + ) + for axis in direction: + if axis == "x": + if 0 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(0) + elif axis == "y": + if 1 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(1) + elif axis == "z": + if 2 in select: + raise ValueError( + "Field.mesh: got duplicate directions {}".format(direction) + ) + select.append(2) + else: + raise ValueError( + "Field.mesh: got unknown direction {}".format(axis) + ) + return select diff --git a/gstools/field/srf.py b/gstools/field/srf.py index befba22c..caa902cd 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -132,20 +132,12 @@ def __call__( field : :class:`numpy.ndarray` the SRF """ - self.mesh_type = mesh_type # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) - # internal conversation - x, y, z, self.pos, mt_gen, mt_changed, axis_lens = self._pre_pos( - pos, mesh_type - ) + # get isometrized positions and the resulting field-shape + iso_pos, shape = self.pre_pos(pos, mesh_type) # generate the field - self.raw_field = self.generator.__call__(x, y, z, mt_gen) - # reshape field if we got an unstructured mesh - if mt_changed: - self.raw_field = reshape_field_from_unstruct_to_struct( - self.model.dim, self.raw_field, axis_lens - ) + self.raw_field = np.reshape(self.generator(iso_pos), shape) # apply given conditions to the field if self.condition: ( diff --git a/gstools/field/tools.py b/gstools/field/tools.py deleted file mode 100644 index 00beb8fd..00000000 --- a/gstools/field/tools.py +++ /dev/null @@ -1,256 +0,0 @@ -# -*- coding: utf-8 -*- -""" -GStools subpackage providing tools for the spatial random field. - -.. currentmodule:: gstools.field.tools - -The following classes and functions are provided - -.. autosummary:: - reshape_input - reshape_input_axis_from_unstruct - reshape_input_axis_from_struct - check_mesh - make_isotropic - make_anisotropic - unrotate_mesh - rotate_mesh - reshape_axis_from_struct_to_unstruct - reshape_field_from_unstruct_to_struct -""" -# pylint: disable=C0103 - -import numpy as np -from gstools.tools.geometric import matrix_rotate, matrix_derotate - -__all__ = [ - "reshape_input", - "reshape_input_axis_from_unstruct", - "reshape_input_axis_from_struct", - "check_mesh", - "make_isotropic", - "make_anisotropic", - "unrotate_mesh", - "rotate_mesh", - "reshape_axis_from_struct_to_unstruct", - "reshape_field_from_unstruct_to_struct", -] - - -# Geometric functions ######################################################### - - -def reshape_input(x, y=None, z=None, mesh_type="unstructured"): - """Reshape given axes, depending on the mesh type.""" - if mesh_type == "unstructured": - x, y, z = reshape_input_axis_from_unstruct(x, y, z) - elif mesh_type == "structured": - x, y, z = reshape_input_axis_from_struct(x, y, z) - return x, y, z - - -def reshape_input_axis_from_unstruct(x, y=None, z=None): - """Reshape given axes for vectorisation on unstructured grid.""" - x = np.atleast_1d(x) - y = np.atleast_1d(y) - z = np.atleast_1d(z) - x = np.reshape(x, (len(x), 1)) - y = np.reshape(y, (len(y), 1)) - z = np.reshape(z, (len(z), 1)) - return (x, y, z) - - -def reshape_input_axis_from_struct(x, y=None, z=None): - """Reshape given axes for vectorisation on unstructured grid.""" - x = np.atleast_1d(x) - y = np.atleast_1d(y) - z = np.atleast_1d(z) - x = np.reshape(x, (len(x), 1, 1, 1)) - y = np.reshape(y, (1, len(y), 1, 1)) - z = np.reshape(z, (1, 1, len(z), 1)) - return (x, y, z) - - -# SRF helpers ################################################################# - - -def check_mesh(dim, x, y, z, mesh_type): - """Do a basic check of the shapes of the input arrays.""" - if dim >= 2: - if y is None: - raise ValueError( - "The y-component is missing for " "{0} dimensions".format(dim) - ) - if dim == 3: - if z is None: - raise ValueError( - "The z-component is missing for " "{0} dimensions".format(dim) - ) - if mesh_type == "unstructured": - if dim >= 2: - try: - if len(x) != len(y): - raise ValueError( - "len(x) = {0} != len(y) = {1} " - "for unstructured grids".format(len(x), len(y)) - ) - except TypeError: - pass - if dim == 3: - try: - if len(x) != len(z): - raise ValueError( - "len(x) = {0} != len(z) = {1} " - "for unstructured grids".format(len(x), len(z)) - ) - except TypeError: - pass - elif mesh_type == "structured": - pass - else: - raise ValueError("Unknown mesh type {0}".format(mesh_type)) - - -def make_isotropic(dim, anis, y, z): - """Stretch given axes in order to implement anisotropy.""" - if dim == 1: - return y, z - if dim == 2: - return y / anis[0], z - if dim == 3: - return y / anis[0], z / anis[1] - return None - - -def make_anisotropic(dim, anis, y, z): - """Re-stretch given axes.""" - if dim == 1: - return y, z - if dim == 2: - return y * anis[0], z - if dim == 3: - return y * anis[0], z * anis[1] - return None - - -def unrotate_mesh(dim, angles, x, y, z): - """Rotate axes in order to implement rotation. - - for 3d: yaw, pitch, and roll angles are alpha, beta, and gamma, - of intrinsic rotation whose Tait-Bryan angles are - alpha, beta, gamma about axes x, y, z. - """ - if dim == 1: - return x, y, z - if dim == 2: - # extract 2d rotation matrix - rot_mat = matrix_derotate(dim, angles) - pos_tuple = np.vstack((x, y)) - pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 2) - x = pos_tuple[0].reshape(np.shape(x)) - y = pos_tuple[1].reshape(np.shape(y)) - return x, y, z - if dim == 3: - rot_mat = matrix_derotate(dim, angles) - pos_tuple = np.vstack((x, y, z)) - pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 3) - x = pos_tuple[0].reshape(np.shape(x)) - y = pos_tuple[1].reshape(np.shape(y)) - z = pos_tuple[2].reshape(np.shape(z)) - return x, y, z - return None - - -def rotate_mesh(dim, angles, x, y, z): - """Rotate axes. - - for 3d: yaw, pitch, and roll angles are alpha, beta, and gamma, - of intrinsic rotation whose Tait-Bryan angles are - alpha, beta, gamma about axes x, y, z. - """ - if dim == 1: - return x, y, z - if dim == 2: - # extract 2d rotation matrix - rot_mat = matrix_rotate(dim, angles) - pos_tuple = np.vstack((x, y)) - pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 2) - x = pos_tuple[0].reshape(np.shape(x)) - y = pos_tuple[1].reshape(np.shape(y)) - return x, y, z - if dim == 3: - rot_mat = matrix_rotate(dim, angles) - pos_tuple = np.vstack((x, y, z)) - pos_tuple = np.vsplit(np.dot(rot_mat, pos_tuple), 3) - x = pos_tuple[0].reshape(np.shape(x)) - y = pos_tuple[1].reshape(np.shape(y)) - z = pos_tuple[2].reshape(np.shape(z)) - return x, y, z - return None - - -def reshape_axis_from_struct_to_unstruct( - dim, x, y=None, z=None, indexing="ij" -): - """Reshape given axes from struct to unstruct for rotation.""" - if dim == 1: - return x, y, z, (len(x),) - if dim == 2: - x_u, y_u = np.meshgrid(x, y, indexing=indexing) - len_unstruct = len(x) * len(y) - x_u = np.reshape(x_u, len_unstruct) - y_u = np.reshape(y_u, len_unstruct) - return x_u, y_u, z, (len(x), len(y)) - if dim == 3: - x_u, y_u, z_u = np.meshgrid(x, y, z, indexing=indexing) - len_unstruct = len(x) * len(y) * len(z) - x_u = np.reshape(x_u, len_unstruct) - y_u = np.reshape(y_u, len_unstruct) - z_u = np.reshape(z_u, len_unstruct) - return x_u, y_u, z_u, (len(x), len(y), len(z)) - return None - - -def reshape_field_from_unstruct_to_struct(dim, field, axis_lens): - """Reshape the rotated field back to struct.""" - if dim == 1: - return field - if dim == 2: - field = np.reshape(field, axis_lens) - return field - if dim == 3: - field = np.reshape(field, axis_lens) - return field - return None - - -def _get_select(direction): - select = [] - if not (0 < len(direction) < 4): - raise ValueError( - "Field.mesh: need 1 to 3 direction(s), got '{}'".format(direction) - ) - for axis in direction: - if axis == "x": - if 0 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(0) - elif axis == "y": - if 1 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(1) - elif axis == "z": - if 2 in select: - raise ValueError( - "Field.mesh: got duplicate directions {}".format(direction) - ) - select.append(2) - else: - raise ValueError( - "Field.mesh: got unknown direction {}".format(axis) - ) - return select From 0ca920cca145fc3db7db57aafd70ef75fbde968d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 15:35:29 +0100 Subject: [PATCH 04/28] Tools: add routines to format pos tuple with given dim or shape --- gstools/tools/geometric.py | 267 +++++++++++++++++++++++++++---------- 1 file changed, 199 insertions(+), 68 deletions(-) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index c11ad3f4..8d062d2f 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -19,8 +19,10 @@ matrix_isometrize matrix_anisometrize rotated_main_axes - pos2xyz - xyz2pos + gen_mesh + format_struct_pos_dim + format_struct_pos_shape + format_unstruct_pos_shape ang2dir """ # pylint: disable=C0103 @@ -40,8 +42,10 @@ "matrix_isometrize", "matrix_anisometrize", "rotated_main_axes", - "pos2xyz", - "xyz2pos", + "gen_mesh", + "format_struct_pos_dim", + "format_struct_pos_shape", + "format_unstruct_pos_shape", "ang2dir", ] @@ -326,87 +330,214 @@ def rotated_main_axes(dim, angles): # conversion ################################################################## -def pos2xyz(pos, dtype=None, calc_dim=False, max_dim=3): - """Convert postional arguments to x, y, z. +def gen_mesh(pos): + """ + Generate mesh from a structured position tuple. + + Parameters + ---------- + pos : :class:`tuple` of :class:`numpy.ndarray` + The structured position tuple. + + Returns + ------- + :class:`numpy.ndarray` + Unstructured position tuple. + """ + return np.array(np.meshgrid(*pos, indexing="ij"), dtype=np.double).reshape( + (len(pos), -1) + ) + + +def format_struct_pos_dim(pos, dim): + """ + Format a structured position tuple with given dimension. Parameters ---------- pos : :any:`iterable` - the position tuple, containing main direction and transversal - directions - dtype : data-type, optional - The desired data-type for the array. - If not given, then the type will be determined as the minimum type - required to hold the objects in the sequence. Default: None - calc_dim : :class:`bool`, optional - State if the dimension should be returned. Default: False - max_dim : :class:`int`, optional - Cut of information above the given dimension. Default: 3 + Position tuple, containing main direction and transversal directions. + dim : :class:`int` + Spatial dimension. + + Raises + ------ + ValueError + When position tuple doesn't match the given dimension. Returns ------- - x : :class:`numpy.ndarray` - first components of position vectors - y : :class:`numpy.ndarray` or None - analog to x - z : :class:`numpy.ndarray` or None - analog to x - dim : :class:`int`, optional - dimension (only if calc_dim is True) + pos : :class:`tuple` of :class:`numpy.ndarray` + The formatted structured position tuple. + shape : :class:`tuple` + Shape of the resulting field. + """ + if dim == 1: + pos = (np.array(pos, dtype=np.double).reshape(-1),) + elif len(pos) != dim: + raise ValueError("Formatting: position tuple doesn't match dimension.") + else: + pos = tuple(np.array(p_i, dtype=np.double).reshape(-1) for p_i in pos) + shape = tuple(len(p_i) for p_i in pos) + return pos, shape - Notes - ----- - If len(pos) > 3, everything after pos[2] will be ignored. + +def format_struct_pos_shape(pos, shape, check_stacked_shape=False): """ - if max_dim == 1: # sanity check - pos = np.array(pos, ndmin=2) - x = np.array(pos[0], dtype=dtype).reshape(-1) - dim = 1 - y = z = None - if len(pos) > 1 and max_dim > 1: - dim = 2 - y = np.array(pos[1], dtype=dtype).reshape(-1) - if len(pos) > 2 and max_dim > 2: - dim = 3 - z = np.array(pos[2], dtype=dtype).reshape(-1) - if calc_dim: - return x, y, z, dim - return x, y, z - - -def xyz2pos(x, y=None, z=None, dtype=None, max_dim=3): - """Convert x, y, z to postional arguments. + Format a structured position tuple with given shape. + + Shape could be stacked, when multiple fields were given. Parameters ---------- - x : :class:`numpy.ndarray` - grid axis in x-direction if structured, or first components of - position vectors if unstructured - y : :class:`numpy.ndarray`, optional - analog to x - z : :class:`numpy.ndarray`, optional - analog to x - dtype : data-type, optional - The desired data-type for the array. - If not given, then the type will be determined as the minimum type - required to hold the objects in the sequence. Default: None - max_dim : :class:`int`, optional - Cut of information above the given dimension. Default: 3 + pos : :any:`iterable` + Position tuple, containing main direction and transversal directions. + shape : :class:`tuple` + Shape of the input field. + check_stacked_shape : :class:`bool`, optional + Whether to check if given shape comes from stacked fields. + Default: False. + + Raises + ------ + ValueError + When position tuple doesn't match the given dimension. + + Returns + ------- + pos : :class:`tuple` of :class:`numpy.ndarray` + The formatted structured position tuple. + shape : :class:`tuple` + Shape of the resulting field. + dim : :class:`int` + Spatial dimension. + """ + # some help from the given shape + shape_size = np.prod(shape) + stacked_shape_size = np.prod(shape[1:]) + wrong_shape = False + # now we try to be smart + try: + # if this works we have either: + # - a 1D array + # - nD array where all axes have same length (corner case) + check_pos = np.array(pos, dtype=np.double, ndmin=2) + except ValueError: + # if it doesn't work, we have a tuple of differently sized axes (easy) + dim = len(pos) + pos, pos_shape = format_struct_pos_dim(pos, dim) + # determine if we have a stacked field if wanted + if check_stacked_shape and stacked_shape_size == np.prod(pos_shape): + shape = (shape[0],) + pos_shape + # check if we have a single field with matching size + elif shape_size == np.prod(pos_shape): + shape = (1,) + pos_shape if check_stacked_shape else pos_shape + # if nothing works, we raise an error + else: + wrong_shape = True + else: + struct_size = np.prod([p_i.size for p_i in check_pos]) + # case: 1D unstacked + if check_pos.size == shape_size: + dim = 1 + pos, pos_shape = format_struct_pos_dim(check_pos, dim) + shape = (1,) + pos_shape if check_stacked_shape else pos_shape + # case: 1D and stacked + elif check_pos.size == stacked_shape_size: + dim = 1 + pos, pos_shape = format_struct_pos_dim(check_pos, dim) + cnt = shape[0] + shape = (cnt,) + pos_shape + wrong_shape = not check_stacked_shape + # case: nD unstacked + elif struct_size == shape_size: + dim = len(check_pos) + pos, pos_shape = format_struct_pos_dim(pos, dim) + shape = (1,) + pos_shape if check_stacked_shape else pos_shape + # case: nD and stacked + elif struct_size == stacked_shape_size: + dim = len(check_pos) + pos, pos_shape = format_struct_pos_dim(pos, dim) + cnt = shape[0] + shape = (cnt,) + pos_shape + wrong_shape = not check_stacked_shape + # if nothing works, we raise an error + else: + wrong_shape = True + + # if shape was wrong at one point we raise an error + if wrong_shape: + raise ValueError("Formatting: position tuple doesn't match dimension.") + + return pos, shape, dim + + +def format_unstruct_pos_shape(pos, shape, check_stacked_shape=False): + """ + Format an unstructured position tuple with given shape. + + Shape could be stacked, when multiple fields were given. + + Parameters + ---------- + pos : :any:`iterable` + Position tuple, containing point coordinates. + shape : :class:`tuple` + Shape of the input field. + check_stacked_shape : :class:`bool`, optional + Whether to check if given shape comes from stacked fields. + Default: False. + + Raises + ------ + ValueError + When position tuple doesn't match the given dimension. Returns ------- pos : :class:`tuple` of :class:`numpy.ndarray` - the position tuple + The formatted structured position tuple. + shape : :class:`tuple` + Shape of the resulting field. + dim : :class:`int` + Spatial dimension. """ - if y is None and z is not None: - raise ValueError("gstools.tools.xyz2pos: if z is given, y is needed!") - pos = [] - pos.append(np.array(x, dtype=dtype).reshape(-1)) - if y is not None and max_dim > 1: - pos.append(np.array(y, dtype=dtype).reshape(-1)) - if z is not None and max_dim > 2: - pos.append(np.array(z, dtype=dtype).reshape(-1)) - return tuple(pos) + # some help from the given shape + shape_size = np.prod(shape) + stacked_shape_size = np.prod(shape[1:]) + wrong_shape = False + # now we try to be smart + pos = np.array(pos, dtype=np.double, ndmin=2) + pnt_cnt = pos[0].size + # case: 1D unstacked + if pos.size == shape_size: + dim = 1 + shape = (1, pos.size) if check_stacked_shape else (pos.size,) + # case: 1D and stacked + elif pos.size == stacked_shape_size: + dim = 1 + shape = (shape[0], pos.size) + wrong_shape = not check_stacked_shape + # case: nD unstacked + elif pnt_cnt == shape_size: + dim = len(pos) + shape = (1, pnt_cnt) if check_stacked_shape else pnt_cnt + # case: nD and stacked + elif pnt_cnt == stacked_shape_size: + dim = len(pos) + shape = (shape[0], pnt_cnt) + wrong_shape = not check_stacked_shape + # if nothing works, we raise an error + else: + wrong_shape = True + + # if shape was wrong at one point we raise an error + if wrong_shape: + raise ValueError("Formatting: position tuple doesn't match dimension.") + + pos = pos.reshape((dim, -1)) + + return pos, shape, dim def ang2dir(angles, dtype=np.double, dim=None): From 88fbb79bb049824db2d468fecc10763004bd2e14 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 15:37:49 +0100 Subject: [PATCH 05/28] Variogram: better care about pos tuple with given field shape --- gstools/variogram/variogram.py | 35 ++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index e0af0d97..db005ce9 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -15,7 +15,12 @@ import numpy as np from gstools.field.tools import reshape_axis_from_struct_to_unstruct -from gstools.tools.geometric import pos2xyz, xyz2pos, ang2dir +from gstools.tools.geometric import ( + gen_mesh, + format_struct_pos_shape, + format_unstruct_pos_shape, + ang2dir, +) from gstools.variogram.estimator import ( unstructured, structured, @@ -209,17 +214,19 @@ def vario_estimate( return bin_centres, estimates if not masked: field = field.filled() - # check_mesh shape - x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) + # check mesh shape if mesh_type != "unstructured": - x, y, z, __ = reshape_axis_from_struct_to_unstruct(dim, x, y, z) - if len(field.shape) > 2 or field.shape[1] != len(x): - try: - field = field.reshape((-1, len(x))) - except ValueError as exc: - raise ValueError("'field' has wrong shape") from exc - # prepare positions - pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) + pos, __, dim = format_struct_pos_shape( + pos, field.shape, check_stacked_shape=True + ) + pos = gen_mesh(pos) + else: + pos, __, dim = format_unstruct_pos_shape( + pos, field.shape, check_stacked_shape=True + ) + # prepare the field + pnt_cnt = len(pos[0]) + field = field.reshape((-1, pnt_cnt)) # apply mask if wanted if masked: # if fields have different masks, take the minimal common mask @@ -228,7 +235,7 @@ def vario_estimate( if np.size(mask) > 1: # not only np.ma.nomask select = np.invert( np.logical_or( - np.reshape(mask, len(x)), np.all(field.mask, axis=0) + np.reshape(mask, pnt_cnt), np.all(field.mask, axis=0) ) ) else: @@ -264,9 +271,9 @@ def vario_estimate( bandwidth = float(bandwidth) if bandwidth is not None else -1.0 angles_tol = float(angles_tol) # prepare sampled variogram - if sampling_size is not None and sampling_size < len(x): + if sampling_size is not None and sampling_size < pnt_cnt: sampled_idx = np.random.RandomState(sampling_seed).choice( - np.arange(len(x)), sampling_size, replace=False + np.arange(pnt_cnt), sampling_size, replace=False ) field = field[:, sampled_idx] pos = pos[:, sampled_idx] From 3ffdf39f3f50b7ea4b5dc92a797600589e5ae9d5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 15:39:45 +0100 Subject: [PATCH 06/28] Field: update pre_pos routine: name of attribute to set selectable --- gstools/field/base.py | 31 ++++++++++++------------------- gstools/field/srf.py | 1 - 2 files changed, 12 insertions(+), 20 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 182dc952..6bf318eb 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -16,6 +16,7 @@ import numpy as np from gstools.covmodel.base import CovModel +from gstools.tools.geometric import format_struct_pos_dim, gen_mesh from gstools.tools.export import to_vtk, vtk_export @@ -147,7 +148,7 @@ def mesh( mesh.point_data[name] = field return out - def pre_pos(self, pos, mesh_type="unstructured"): + def pre_pos(self, pos, mesh_type="unstructured", name="pos"): """ Preprocessing positions and mesh_type. @@ -156,8 +157,12 @@ def pre_pos(self, pos, mesh_type="unstructured"): pos : :any:`iterable` the position tuple, containing main direction and transversal directions - mesh_type : :class:`str` + mesh_type : :class:`str`, optional 'structured' / 'unstructured' + Default: `"unstructured"` + name : :class:`str`, optional + Name of the attribute to be set in the class. + Default: `"pos"` Returns ------- @@ -170,25 +175,13 @@ def pre_pos(self, pos, mesh_type="unstructured"): self.mesh_type = mesh_type # save pos tuple if mesh_type != "unstructured": - if self.model.dim == 1: - self.pos = (np.array(pos, dtype=np.double).reshape(-1),) - else: - if len(pos) != self.model.dim: - raise ValueError( - "Field: position tuple doesn't match dimension." - ) - self.pos = tuple( - [ - np.array(pos_i, dtype=np.double).reshape(-1) - for pos_i in pos - ] - ) - pos = np.meshgrid(*self.pos, indexing="ij") - shape = tuple([len(pos_i) for pos_i in self.pos]) + pos, shape = format_struct_pos_dim(pos, self.model.dim) + setattr(self, name, pos) + pos = gen_mesh(pos) else: pos = np.array(pos, dtype=np.double).reshape(self.model.dim, -1) - self.pos = pos - shape = np.shape(self.pos[0]) + setattr(self, name, pos) + shape = np.shape(pos[0]) # return isometrized pos tuple and resulting field shape return self.model.isometrize(pos), shape diff --git a/gstools/field/srf.py b/gstools/field/srf.py index caa902cd..c1a189cc 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -13,7 +13,6 @@ import numpy as np from gstools.field.generator import RandMeth, IncomprRandMeth -from gstools.field.tools import reshape_field_from_unstruct_to_struct from gstools.field.base import Field from gstools.field.upscaling import var_coarse_graining, var_no_scaling from gstools.field.condition import ordinary, simple From 898c8b14b0705237231824332365b958d3bd0aac Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 15:40:19 +0100 Subject: [PATCH 07/28] CovModel.plot: update spatial plotting routines to use new Field routines --- gstools/covmodel/plot.py | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/gstools/covmodel/plot.py b/gstools/covmodel/plot.py index 66084aa5..4cd4aca7 100644 --- a/gstools/covmodel/plot.py +++ b/gstools/covmodel/plot.py @@ -24,7 +24,7 @@ import numpy as np import gstools -from gstools.field.tools import reshape_axis_from_struct_to_unstruct + __all__ = [ "plot_variogram", @@ -70,13 +70,11 @@ def plot_vario_spatial( field._value_type = "scalar" if x_max is None: x_max = 3 * model.len_scale - field.mesh_type = "structured" x_s = np.linspace(-x_max, x_max) + x_min - pos = [x_s] * model.dim - x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos) - vario = model.vario_spatial([x, y, z][: model.dim]).reshape(shape) - field.pos = pos - field.field = vario + iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured") + field.field = model.vario_spatial(model.anisometrize(iso_pos)).reshape( + shape + ) return field.plot(fig=fig, ax=ax) @@ -88,13 +86,9 @@ def plot_cov_spatial( field._value_type = "scalar" if x_max is None: x_max = 3 * model.len_scale - field.mesh_type = "structured" x_s = np.linspace(-x_max, x_max) + x_min - pos = [x_s] * model.dim - x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos) - vario = model.cov_spatial([x, y, z][: model.dim]).reshape(shape) - field.pos = pos - field.field = vario + iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured") + field.field = model.cov_spatial(model.anisometrize(iso_pos)).reshape(shape) return field.plot(fig=fig, ax=ax) @@ -106,13 +100,9 @@ def plot_cor_spatial( field._value_type = "scalar" if x_max is None: x_max = 3 * model.len_scale - field.mesh_type = "structured" x_s = np.linspace(-x_max, x_max) + x_min - pos = [x_s] * model.dim - x, y, z, shape = reshape_axis_from_struct_to_unstruct(model.dim, *pos) - vario = model.cor_spatial([x, y, z][: model.dim]).reshape(shape) - field.pos = pos - field.field = vario + iso_pos, shape = field.pre_pos([x_s] * model.dim, "structured") + field.field = model.cor_spatial(model.anisometrize(iso_pos)).reshape(shape) return field.plot(fig=fig, ax=ax) From f9fb673d8ab660fd55d0c62fff20a7e07cfeff4a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:52:34 +0100 Subject: [PATCH 08/28] VTK: refactor VTK export to drop usage of pos2xyz --- gstools/tools/__init__.py | 6 ------ gstools/tools/export.py | 27 ++++++++++++++++----------- 2 files changed, 16 insertions(+), 17 deletions(-) diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index bd80462c..4d8963d9 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -43,8 +43,6 @@ matrix_anisotropify matrix_isometrize matrix_anisometrize - pos2xyz - xyz2pos ang2dir ---- @@ -82,8 +80,6 @@ matrix_isometrize, matrix_anisometrize, rotated_main_axes, - pos2xyz, - xyz2pos, ang2dir, ) @@ -113,7 +109,5 @@ "matrix_isometrize", "matrix_anisometrize", "rotated_main_axes", - "pos2xyz", - "xyz2pos", "ang2dir", ] diff --git a/gstools/tools/export.py b/gstools/tools/export.py index 7085ac05..92d8dd03 100644 --- a/gstools/tools/export.py +++ b/gstools/tools/export.py @@ -19,7 +19,6 @@ import numpy as np from pyevtk.hl import gridToVTK, pointsToVTK -from gstools.tools.geometric import pos2xyz try: import pyvista as pv @@ -43,11 +42,14 @@ def _vtk_structured_helper(pos, fields): """Extract field info for vtk rectilinear grid.""" if not isinstance(fields, dict): fields = {"field": fields} - x, y, z = pos2xyz(pos) - if y is None: - y = np.array([0]) - if z is None: - z = np.array([0]) + if len(pos) > 3: + raise ValueError( + "gstools.vtk_export_structured: " + "vtk export only possible for dim=1,2,3" + ) + x = pos[0] + y = pos[1] if len(pos) > 1 else np.array([0]) + z = pos[2] if len(pos) > 2 else np.array([0]) # need fortran order in VTK for field in fields: fields[field] = fields[field].reshape(-1, order="F") @@ -110,11 +112,14 @@ def vtk_export_structured(filename, pos, fields): # pragma: no cover def _vtk_unstructured_helper(pos, fields): if not isinstance(fields, dict): fields = {"field": fields} - x, y, z = pos2xyz(pos) - if y is None: - y = np.zeros_like(x) - if z is None: - z = np.zeros_like(x) + if len(pos) > 3: + raise ValueError( + "gstools.vtk_export_structured: " + "vtk export only possible for dim=1,2,3" + ) + x = pos[0] + y = pos[1] if len(pos) > 1 else np.zeros_like(x) + z = pos[2] if len(pos) > 2 else np.zeros_like(x) for field in fields: fields[field] = fields[field].reshape(-1) if ( From 21d3fb5bd0b6b49b7b0f3e846dfbb48661a55bc1 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:56:07 +0100 Subject: [PATCH 09/28] SRF: update conditioned SRF to drop pos2xyz usage --- gstools/field/condition.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/gstools/field/condition.py b/gstools/field/condition.py index 9d24ff2d..c107c978 100644 --- a/gstools/field/condition.py +++ b/gstools/field/condition.py @@ -11,8 +11,6 @@ simple """ # pylint: disable=C0103 -from gstools.field.tools import make_isotropic, unrotate_mesh -from gstools.tools.geometric import pos2xyz from gstools.krige import Ordinary, Simple @@ -43,11 +41,7 @@ def ordinary(srf): krige_field, krige_var = krige_ok(srf.pos, srf.mesh_type) # evaluate the field at the conditional points - x, y, z = pos2xyz(srf.cond_pos, max_dim=srf.model.dim) - if srf.model.do_rotation: - x, y, z = unrotate_mesh(srf.model.dim, srf.model.angles, x, y, z) - y, z = make_isotropic(srf.model.dim, srf.model.anis, y, z) - err_data = srf.generator.__call__(x, y, z, "unstructured") + err_data = srf.generator(srf.model.isometrize(srf.cond_pos)) err_ok = Ordinary( model=srf.model, cond_pos=srf.cond_pos, cond_val=err_data @@ -88,11 +82,7 @@ def simple(srf): krige_field, krige_var = krige_sk(srf.pos, srf.mesh_type) # evaluate the field at the conditional points - x, y, z = pos2xyz(srf.cond_pos, max_dim=srf.model.dim) - if srf.model.do_rotation: - x, y, z = unrotate_mesh(srf.model.dim, srf.model.angles, x, y, z) - y, z = make_isotropic(srf.model.dim, srf.model.anis, y, z) - err_data = srf.generator.__call__(x, y, z, "unstructured") + srf.mean + err_data = srf.generator(srf.model.isometrize(srf.cond_pos)) + srf.mean err_sk = Simple( model=srf.model, From cffc8300cc1b3bcf8e6e2926730f393623427d60 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:57:22 +0100 Subject: [PATCH 10/28] Field.pre_pos: generate correct shape for vector fields --- gstools/field/base.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 6bf318eb..a21da04b 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -148,7 +148,7 @@ def mesh( mesh.point_data[name] = field return out - def pre_pos(self, pos, mesh_type="unstructured", name="pos"): + def pre_pos(self, pos, mesh_type="unstructured"): """ Preprocessing positions and mesh_type. @@ -160,9 +160,6 @@ def pre_pos(self, pos, mesh_type="unstructured", name="pos"): mesh_type : :class:`str`, optional 'structured' / 'unstructured' Default: `"unstructured"` - name : :class:`str`, optional - Name of the attribute to be set in the class. - Default: `"pos"` Returns ------- @@ -176,12 +173,15 @@ def pre_pos(self, pos, mesh_type="unstructured", name="pos"): # save pos tuple if mesh_type != "unstructured": pos, shape = format_struct_pos_dim(pos, self.model.dim) - setattr(self, name, pos) + self.pos = pos pos = gen_mesh(pos) else: pos = np.array(pos, dtype=np.double).reshape(self.model.dim, -1) - setattr(self, name, pos) + self.pos = pos shape = np.shape(pos[0]) + # prepend dimension if we have a vector field + if self.value_type == "vector": + shape = (self.model.dim,) + shape # return isometrized pos tuple and resulting field shape return self.model.isometrize(pos), shape From ecdf5272d435f3a5712274e4bfaeb70f4c705db4 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:58:25 +0100 Subject: [PATCH 11/28] SRF: correct shape handling for upscaling variance --- gstools/field/srf.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index c1a189cc..f63857a9 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -158,6 +158,8 @@ def __call__( # upscaled variance if not np.isscalar(point_volumes) or not np.isclose(point_volumes, 0): scaled_var = self.upscaling_func(self.model, point_volumes) + if np.size(scaled_var) > 1: + scaled_var = np.reshape(scaled_var, shape) self.field -= self.mean self.field *= np.sqrt(scaled_var / self.model.sill) self.field += self.mean From 8eeeaa561a0d9eb1da095ac5fa274023605d877e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:58:51 +0100 Subject: [PATCH 12/28] Field.plot: drop usage of pos2xyz --- gstools/field/plot.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 73f5b7e8..97520f79 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -15,7 +15,6 @@ from scipy import interpolate as inter import matplotlib.pyplot as plt from matplotlib.widgets import Slider, RadioButtons -from gstools.tools import pos2xyz from gstools.covmodel.plot import _get_fig_ax __all__ = ["plot_field", "plot_vec_field"] @@ -56,7 +55,7 @@ def _plot_1d(pos, field, fig=None, ax=None): # pragma: no cover """Plot a 1d field.""" fig, ax = _get_fig_ax(fig, ax) title = "Field 1D: " + str(field.shape) - x, __, __ = pos2xyz(pos, max_dim=1) + x = pos[0] x = x.flatten() arg = np.argsort(x) ax.plot(x[arg], field.ravel()[arg]) @@ -71,7 +70,8 @@ def _plot_2d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover """Plot a 2d field.""" fig, ax = _get_fig_ax(fig, ax) title = "Field 2D " + mesh_type + ": " + str(field.shape) - x, y, __ = pos2xyz(pos, max_dim=2) + x = pos[0] + y = pos[1] if mesh_type == "unstructured": cont = ax.tricontourf(x, y, field.ravel(), levels=256) else: @@ -241,7 +241,8 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover fig, ax = _get_fig_ax(fig, ax) title = "Field 2D " + fld.mesh_type + ": " + str(plot_field.shape) - x, y, __ = pos2xyz(fld.pos, max_dim=2) + x = fld.pos[0] + y = fld.pos[1] sp = plt.streamplot( x, From 4c929c47bb963602d3c8fb3a1580fa4794796395 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 17:59:17 +0100 Subject: [PATCH 13/28] Variogram: cleanup --- gstools/variogram/variogram.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index db005ce9..c7bdb4f0 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -14,7 +14,6 @@ import numpy as np -from gstools.field.tools import reshape_axis_from_struct_to_unstruct from gstools.tools.geometric import ( gen_mesh, format_struct_pos_shape, From 97ef9f104623cea6e4db1cb2e40233c8ca6fee5d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 18:01:38 +0100 Subject: [PATCH 14/28] Krige: refactor Kriging routines to drop usage of pos2xyz --- gstools/krige/base.py | 74 ++++++++++++++---------------------------- gstools/krige/tools.py | 47 ++++++++++++--------------- 2 files changed, 44 insertions(+), 77 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 57cd4121..c468bcb6 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -15,12 +15,6 @@ from scipy.spatial.distance import cdist import scipy.linalg as spl -from gstools.field.tools import ( - reshape_field_from_unstruct_to_struct, - make_anisotropic, - rotate_mesh, -) -from gstools.tools.geometric import pos2xyz, xyz2pos from gstools.field.base import Field from gstools.krige.krigesum import krigesum from gstools.krige.tools import ( @@ -175,48 +169,40 @@ def __call__( krige_var : :class:`numpy.ndarray` the kriging error variance """ - self.mesh_type = mesh_type - # internal conversation - x, y, z, self.pos, __, mt_changed, axis_lens = self._pre_pos( - pos, mesh_type, make_unstruct=True - ) - field = np.empty_like(x) - krige_var = np.empty_like(x) + iso_pos, shape = self.pre_pos(pos, mesh_type) + pnt_cnt = len(iso_pos[0]) + + field = np.empty(pnt_cnt, dtype=np.double) + krige_var = np.empty(pnt_cnt, dtype=np.double) # set constant mean if present and wanted if only_mean and self.has_const_mean: field[...] = self.get_mean() - self.mean # mean is added later krige_var[...] = self.model.sill # will result in 0 var. later # execute the kriging routine else: - point_no = len(x) # set chunk size - chunk_size = point_no if chunk_size is None else int(chunk_size) - chunk_no = int(np.ceil(point_no / chunk_size)) - ext_drift = self._pre_ext_drift(point_no, ext_drift) + chunk_size = pnt_cnt if chunk_size is None else int(chunk_size) + chunk_no = int(np.ceil(pnt_cnt / chunk_size)) + ext_drift = self._pre_ext_drift(pnt_cnt, ext_drift) # iterate of chunks for i in range(chunk_no): # get chunk slice for actual chunk chunk_slice = ( i * chunk_size, - min(point_no, (i + 1) * chunk_size), + min(pnt_cnt, (i + 1) * chunk_size), ) c_slice = slice(*chunk_slice) # get RHS of the kriging system k_vec = self._get_krige_vecs( - (x, y, z), chunk_slice, ext_drift, only_mean + iso_pos, chunk_slice, ext_drift, only_mean ) # generate the raw kriging field and error variance field[c_slice], krige_var[c_slice] = krigesum( self._krige_mat, k_vec, self._krige_cond ) # reshape field if we got a structured mesh - if mt_changed: - field = reshape_field_from_unstruct_to_struct( - self.model.dim, field, axis_lens - ) - krige_var = reshape_field_from_unstruct_to_struct( - self.model.dim, krige_var, axis_lens - ) + field = np.reshape(field, shape) + krige_var = np.reshape(krige_var, shape) self._post_field(field, krige_var) return self.field, self.krige_var @@ -280,17 +266,9 @@ def _get_krige_vecs( if self.unbiased: res[self.cond_no, :] = 1 # drift function need the anisotropic and rotated positions - if self.int_drift_no > 0 and not self.model.is_isotropic: - x, y, z = pos2xyz(pos, max_dim=self.model.dim) - y, z = make_anisotropic(self.model.dim, self.model.anis, y, z) - if self.model.do_rotation: - x, y, z = rotate_mesh( - self.model.dim, self.model.angles, x, y, z - ) - pos = xyz2pos(x, y, z, max_dim=self.model.dim) - chunk_pos = list(pos[: self.model.dim]) - for i in range(self.model.dim): - chunk_pos[i] = chunk_pos[i][slice(*chunk_slice)] + if self.int_drift_no > 0: + pos = self.model.anisometrize(pos) + chunk_pos = pos[:, slice(*chunk_slice)] # apply functional drift for i, f in enumerate(self.drift_functions): res[-self.drift_no + i, :] = f(*chunk_pos) @@ -300,13 +278,13 @@ def _get_krige_vecs( res[ext_size:, :] = ext_drift[:, slice(*chunk_slice)] return res - def _pre_ext_drift(self, point_no, ext_drift=None, set_cond=False): + def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): """ Preprocessor for external drifts. Parameters ---------- - point_no : :class:`numpy.ndarray` + pnt_cnt : :class:`numpy.ndarray` Number of points of the mesh. ext_drift : :class:`numpy.ndarray` or :any:`None`, optional the external drift values at the given positions (only for EDK) @@ -324,12 +302,12 @@ def _pre_ext_drift(self, point_no, ext_drift=None, set_cond=False): if ext_drift is not None: if set_cond: ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) - if len(ext_drift.shape) > 2 or ext_drift.shape[1] != point_no: - raise ValueError("Krige: wrong number of cond. drifts.") + if len(ext_drift.shape) > 2 or ext_drift.shape[1] != pnt_cnt: + raise ValueError("Krige: wrong number of ext. drifts.") return ext_drift ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) ext_shape = np.shape(ext_drift) - shape = (self.drift_no, point_no) + shape = (self.drift_no, pnt_cnt) if self.drift_no > 1 and ext_shape[0] != self.drift_no: raise ValueError("Krige: wrong number of external drifts.") if np.prod(ext_shape) != np.prod(shape): @@ -354,7 +332,7 @@ def _post_field(self, field, krige_var): self.field = field else: self.field = field + eval_func( - self.trend_function, self.pos, self.mesh_type + self.trend_function, self.pos, self.model.dim, self.mesh_type ) self.field += self.mean self.krige_var = self.model.sill - krige_var @@ -377,12 +355,9 @@ def _get_dists(self, pos1, pos2=None, pos2_slice=(0, None)): :class:`numpy.ndarray` Matrix containing the pairwise distances. """ - pos1_stack = np.column_stack(pos1[: self.model.dim]) if pos2 is None: - return cdist(pos1_stack, pos1_stack) - p2s = slice(*pos2_slice) - pos2_stack = np.column_stack(pos2[: self.model.dim])[p2s, ...] - return cdist(pos1_stack, pos2_stack) + return cdist(pos1.T, pos1.T) + return cdist(pos1.T, pos2.T[slice(*pos2_slice), ...]) def get_mean(self): """Calculate the estimated mean.""" @@ -478,9 +453,8 @@ def set_drift_functions(self, drift_functions=None): def update(self): """Update the kriging settings.""" - x, y, z, __, __, __, __ = self._pre_pos(self.cond_pos) + self._krige_pos = self.model.isometrize(self.cond_pos) # krige pos are the unrotated and isotropic condition positions - self._krige_pos = (x, y, z)[: self.model.dim] self._krige_mat = self._get_krige_mat() if self.trend_function is no_trend: self._cond_trend = 0.0 diff --git a/gstools/krige/tools.py b/gstools/krige/tools.py index 121e312b..64448f54 100644 --- a/gstools/krige/tools.py +++ b/gstools/krige/tools.py @@ -15,11 +15,8 @@ # pylint: disable=C0103 from itertools import combinations_with_replacement import numpy as np -from gstools.tools.geometric import pos2xyz, xyz2pos -from gstools.field.tools import ( - reshape_axis_from_struct_to_unstruct, - reshape_field_from_unstruct_to_struct, -) +from gstools.tools.geometric import format_struct_pos_dim, gen_mesh + __all__ = ["no_trend", "eval_func", "set_condition", "get_drift_functions"] @@ -44,17 +41,19 @@ def no_trend(*args, **kwargs): return 0.0 -def eval_func(func, pos, mesh_type="structured"): +def eval_func(func, pos, dim, mesh_type="structured"): """ Evaluate a function on a mesh. Parameters ---------- func : :any:`callable` - The function to be called. Should have the signiture f(x, [y, z]) + The function to be called. Should have the signiture f(x, [y, z, ...]). pos : :class:`list` - the position tuple, containing main direction and transversal - directions (x, [y, z]) + The position tuple, containing main direction and transversal + directions (x, [y, z, ...]). + dim : :class:`int` + The spatial dimension. mesh_type : :class:`str`, optional 'structured' / 'unstructured' @@ -63,16 +62,16 @@ def eval_func(func, pos, mesh_type="structured"): :class:`numpy.ndarray` Function values at the given points. """ - x, y, z, dim = pos2xyz(pos, calc_dim=True) - if mesh_type == "structured": - x, y, z, axis_lens = reshape_axis_from_struct_to_unstruct(dim, x, y, z) - res = func(*[x, y, z][:dim]) - if mesh_type == "structured": - res = reshape_field_from_unstruct_to_struct(dim, res, axis_lens) - return res + if mesh_type != "unstructured": + pos, shape = format_struct_pos_dim(pos, dim) + pos = gen_mesh(pos) + else: + pos = np.array(pos, dtype=np.double).reshape(dim, -1) + shape = np.shape(pos[0]) + return np.reshape(func(*pos), shape) -def set_condition(cond_pos, cond_val, max_dim=3): +def set_condition(cond_pos, cond_val, dim): """ Set the conditions for kriging. @@ -82,8 +81,8 @@ def set_condition(cond_pos, cond_val, max_dim=3): the position tuple of the conditions (x, [y, z]) cond_val : :class:`numpy.ndarray` the values of the conditions - max_dim : :class:`int`, optional - Cut of information above the given dimension. Default: 3 + dim : :class:`int`, optional + Spatial dimension Raises ------ @@ -98,15 +97,9 @@ def set_condition(cond_pos, cond_val, max_dim=3): the error checked cond_val """ # convert the input for right shapes and dimension checks - c_x, c_y, c_z = pos2xyz(cond_pos, dtype=np.double, max_dim=max_dim) - cond_pos = xyz2pos(c_x, c_y, c_z) - if len(cond_pos) != max_dim: - raise ValueError( - "Please check your 'cond_pos' parameters. " - + "The dimension does not match with the given one." - ) cond_val = np.array(cond_val, dtype=np.double).reshape(-1) - if not all([len(cond_pos[i]) == len(cond_val) for i in range(max_dim)]): + cond_pos = np.array(cond_pos, dtype=np.double).reshape(dim, -1) + if len(cond_pos[0]) != len(cond_val): raise ValueError( "Please check your 'cond_pos' and 'cond_val' parameters. " + "The shapes do not match." From 0380f06374541380198b324c86ada87b7804104e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 18:02:25 +0100 Subject: [PATCH 15/28] Tests: skip some structured testing (dropped); fix some input missmatches --- tests/test_incomprrandmeth.py | 60 +++++++++++++++++------------------ tests/test_krige.py | 2 +- tests/test_randmeth.py | 58 ++++++++++++++++----------------- tests/test_srf.py | 15 +++++---- 4 files changed, 68 insertions(+), 67 deletions(-) diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index 2ef39d5c..72167f79 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -31,47 +31,47 @@ def setUp(self): ) def test_unstruct_2d(self): - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0, 0], 0.50751115) self.assertAlmostEqual(modes[0, 1], 1.03291018) self.assertAlmostEqual(modes[1, 1], -0.22003005) def test_unstruct_3d(self): - modes = self.rm_3d(self.x_tuple, self.y_tuple, self.z_tuple) + modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) self.assertAlmostEqual(modes[0, 0], 1.49469700) self.assertAlmostEqual(modes[0, 1], 1.38687858) self.assertAlmostEqual(modes[1, 0], -0.27245271) - def test_struct_2d(self): - modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") - self.assertAlmostEqual(modes[0, 0, 0], 0.50751115) - self.assertAlmostEqual(modes[0, 1, 0], 0.69751927) - self.assertAlmostEqual(modes[1, 1, 1], -0.19747468) + # def test_struct_2d(self): + # modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") + # self.assertAlmostEqual(modes[0, 0, 0], 0.50751115) + # self.assertAlmostEqual(modes[0, 1, 0], 0.69751927) + # self.assertAlmostEqual(modes[1, 1, 1], -0.19747468) - def test_struct_3d(self): - modes = self.rm_3d( - self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" - ) - self.assertAlmostEqual(modes[0, 0, 0, 0], 1.49469700) - self.assertAlmostEqual(modes[1, 0, 1, 1], 0.12813365) - self.assertAlmostEqual(modes[1, 1, 0, 1], 0.01443056) - self.assertAlmostEqual(modes[1, 1, 1, 1], -0.12304040) + # def test_struct_3d(self): + # modes = self.rm_3d( + # self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" + # ) + # self.assertAlmostEqual(modes[0, 0, 0, 0], 1.49469700) + # self.assertAlmostEqual(modes[1, 0, 1, 1], 0.12813365) + # self.assertAlmostEqual(modes[1, 1, 0, 1], 0.01443056) + # self.assertAlmostEqual(modes[1, 1, 1, 1], -0.12304040) - def test_struct_unstruct(self): - x_grid = np.arange(0.0, 2.0, 1.0) - y_grid = np.arange(0.0, 2.0, 1.0) - x_tuple = np.array((0.0, 0.0, 1.0, 1.0)) - y_tuple = np.array((0.0, 1.0, 0.0, 1.0)) - unstr_modes = self.rm_2d(x_tuple, y_tuple, mesh_type="unstructured") - str_modes = self.rm_2d(x_grid, y_grid, mesh_type="structured") - for d in range(2): - k = 0 - for i in range(len(x_grid)): - for j in range(len(y_grid)): - self.assertAlmostEqual( - str_modes[d, i, j], unstr_modes[d, k] - ) - k += 1 + # def test_struct_unstruct(self): + # x_grid = np.arange(0.0, 2.0, 1.0) + # y_grid = np.arange(0.0, 2.0, 1.0) + # x_tuple = np.array((0.0, 0.0, 1.0, 1.0)) + # y_tuple = np.array((0.0, 1.0, 0.0, 1.0)) + # unstr_modes = self.rm_2d(x_tuple, y_tuple, mesh_type="unstructured") + # str_modes = self.rm_2d(x_grid, y_grid, mesh_type="structured") + # for d in range(2): + # k = 0 + # for i in range(len(x_grid)): + # for j in range(len(y_grid)): + # self.assertAlmostEqual( + # str_modes[d, i, j], unstr_modes[d, k] + # ) + # k += 1 def test_assertions(self): cov_model_1d = Gaussian(dim=1, var=1.5, len_scale=2.5) diff --git a/tests/test_krige.py b/tests/test_krige.py index 7b8db5aa..e0e8a762 100644 --- a/tests/test_krige.py +++ b/tests/test_krige.py @@ -213,7 +213,7 @@ def test_pseudo(self): krig = krige.Krige( model, self.p_data[:dim], self.p_vals, unbiased=False ) - field, __ = krig([0, 0, 0]) + field, __ = krig([0, 0, 0][:dim]) # with the pseudo-inverse, the estimated value # should be the mean of the 3 redundant input values self.assertAlmostEqual( diff --git a/tests/test_randmeth.py b/tests/test_randmeth.py index 85b2e691..376aaff4 100644 --- a/tests/test_randmeth.py +++ b/tests/test_randmeth.py @@ -30,64 +30,64 @@ def setUp(self): self.rm_3d = RandMeth(self.cov_model_3d, 100, self.seed) def test_unstruct_1d(self): - modes = self.rm_1d(self.x_tuple) + modes = self.rm_1d((self.x_tuple,)) self.assertAlmostEqual(modes[0], 3.19799030) self.assertAlmostEqual(modes[1], 2.44848295) def test_unstruct_2d(self): - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], 1.67318010) self.assertAlmostEqual(modes[1], 2.12310269) def test_unstruct_3d(self): - modes = self.rm_3d(self.x_tuple, self.y_tuple, self.z_tuple) + modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) self.assertAlmostEqual(modes[0], 0.55488481) self.assertAlmostEqual(modes[1], 1.18506639) - def test_struct_1d(self): - modes = self.rm_1d(self.x_grid, mesh_type="structured") - self.assertAlmostEqual(modes[0], 3.19799030) - self.assertAlmostEqual(modes[1], 2.34788923) - - def test_struct_2d(self): - modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") - self.assertAlmostEqual(modes[0, 0], 1.67318010) - self.assertAlmostEqual(modes[1, 0], 1.54740003) - self.assertAlmostEqual(modes[0, 1], 2.02106551) - self.assertAlmostEqual(modes[1, 1], 1.86883255) - - def test_struct_3d(self): - modes = self.rm_3d( - self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" - ) - self.assertAlmostEqual(modes[0, 0, 0], 0.55488481) - self.assertAlmostEqual(modes[0, 1, 0], 0.41858766) - self.assertAlmostEqual(modes[1, 1, 0], 0.95133855) - self.assertAlmostEqual(modes[0, 1, 1], 0.65475042) - self.assertAlmostEqual(modes[1, 1, 1], 1.40915120) + # def test_struct_1d(self): + # modes = self.rm_1d(self.x_grid, mesh_type="structured") + # self.assertAlmostEqual(modes[0], 3.19799030) + # self.assertAlmostEqual(modes[1], 2.34788923) + + # def test_struct_2d(self): + # modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") + # self.assertAlmostEqual(modes[0, 0], 1.67318010) + # self.assertAlmostEqual(modes[1, 0], 1.54740003) + # self.assertAlmostEqual(modes[0, 1], 2.02106551) + # self.assertAlmostEqual(modes[1, 1], 1.86883255) + + # def test_struct_3d(self): + # modes = self.rm_3d( + # self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" + # ) + # self.assertAlmostEqual(modes[0, 0, 0], 0.55488481) + # self.assertAlmostEqual(modes[0, 1, 0], 0.41858766) + # self.assertAlmostEqual(modes[1, 1, 0], 0.95133855) + # self.assertAlmostEqual(modes[0, 1, 1], 0.65475042) + # self.assertAlmostEqual(modes[1, 1, 1], 1.40915120) def test_reset(self): - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], 1.67318010) self.assertAlmostEqual(modes[1], 2.12310269) self.rm_2d.seed = self.rm_2d.seed - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], 1.67318010) self.assertAlmostEqual(modes[1], 2.12310269) self.rm_2d.seed = 74893621 - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], -1.94278053) self.assertAlmostEqual(modes[1], -1.12401651) self.rm_1d.model = self.cov_model_3d - modes = self.rm_1d(self.x_tuple, self.y_tuple, self.z_tuple) + modes = self.rm_1d((self.x_tuple, self.y_tuple, self.z_tuple)) self.assertAlmostEqual(modes[0], 0.55488481) self.assertAlmostEqual(modes[1], 1.18506639) self.rm_2d.mode_no = 800 - modes = self.rm_2d(self.x_tuple, self.y_tuple) + modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], -3.20809251) self.assertAlmostEqual(modes[1], -2.62032778) diff --git a/tests/test_srf.py b/tests/test_srf.py index 082a7ec3..e5efa107 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -304,13 +304,14 @@ def test_assertions(self): self.assertRaises( ValueError, srf, [self.x_grid, self.y_grid, self.z_grid] ) - self.assertRaises( - ValueError, - srf, - [self.x_tuple, self.y_tuple, self.z_tuple], - self.seed, - mesh_type="hyper_mesh", - ) + # everything not "unstructured" is treated as "structured" + # self.assertRaises( + # ValueError, + # srf, + # [self.x_tuple, self.y_tuple, self.z_tuple], + # self.seed, + # mesh_type="hyper_mesh", + # ) if __name__ == "__main__": From 1081a14be9c2023b9cd098d28eca2db60f08a797 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 18:52:16 +0100 Subject: [PATCH 16/28] Field.plot: only plot for dim=1,2,3 --- gstools/field/plot.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 97520f79..29a77f81 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -46,8 +46,10 @@ def plot_field(fld, field="field", fig=None, ax=None): # pragma: no cover ax = _plot_1d(fld.pos, plot_field, fig, ax) elif fld.model.dim == 2: ax = _plot_2d(fld.pos, plot_field, fld.mesh_type, fig, ax) - else: + elif fld.model.dim == 3: ax = _plot_3d(fld.pos, plot_field, fld.mesh_type, fig, ax) + else: + raise ValueError("Field.plot: only possible for dim=1,2,3!") return ax From 7523ae0b10cff67d54c3d6ea0de04738f60a9921 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 18:52:53 +0100 Subject: [PATCH 17/28] COMMIT #1000 HOORAY! --- gstools/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/__init__.py b/gstools/__init__.py index bbf0ff14..c29faee1 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -104,7 +104,7 @@ vario_estimate vario_estimate_axis """ - +# Hooray! from gstools import field, variogram, random, covmodel, tools, krige, transform from gstools.field import SRF from gstools.tools import ( From 18c0c02d6cdb263b33cb38a441849fc49c72ea5a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 22:17:45 +0100 Subject: [PATCH 18/28] Examples: update variogram fitting plots with scatter plots; add models to comparison --- examples/03_variogram/00_fit_variogram.py | 2 +- examples/03_variogram/02_find_best_model.py | 20 +++++++++++--------- examples/03_variogram/04_directional_2d.py | 4 ++-- examples/03_variogram/05_directional_3d.py | 12 ++++++------ 4 files changed, 20 insertions(+), 18 deletions(-) diff --git a/examples/03_variogram/00_fit_variogram.py b/examples/03_variogram/00_fit_variogram.py index c1ddae62..fbaf5a8e 100644 --- a/examples/03_variogram/00_fit_variogram.py +++ b/examples/03_variogram/00_fit_variogram.py @@ -30,5 +30,5 @@ # Plot the fitting result. ax = fit_model.plot(x_max=40) -ax.plot(bin_center, gamma) +ax.scatter(bin_center, gamma) print(fit_model) diff --git a/examples/03_variogram/02_find_best_model.py b/examples/03_variogram/02_find_best_model.py index 1e7b9f87..dc94b0ef 100755 --- a/examples/03_variogram/02_find_best_model.py +++ b/examples/03_variogram/02_find_best_model.py @@ -25,14 +25,16 @@ # Define a set of models to test. models = { - "gaussian": gs.Gaussian, - "exponential": gs.Exponential, - "matern": gs.Matern, - "stable": gs.Stable, - "rational": gs.Rational, - "linear": gs.Linear, - "circular": gs.Circular, - "spherical": gs.Spherical, + "Gaussian": gs.Gaussian, + "Exponential": gs.Exponential, + "Matern": gs.Matern, + "Stable": gs.Stable, + "Rational": gs.Rational, + "Linear": gs.Linear, + "Circular": gs.Circular, + "Spherical": gs.Spherical, + "SuperSpherical": gs.SuperSpherical, + "JBessel": gs.JBessel, } scores = {} @@ -40,7 +42,7 @@ # Iterate over all models, fit their variogram and calculate the r2 score. # plot the estimated variogram -plt.scatter(bin_center, gamma, label="data") +plt.scatter(bin_center, gamma, color="k", label="data") ax = plt.gca() # fit all models to the estimated variogram diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index 7d5dd55d..c74dc7c0 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -48,8 +48,8 @@ fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5]) -ax1.plot(bin_center, dir_vario[0], label="emp. vario: pi/8") -ax1.plot(bin_center, dir_vario[1], label="emp. vario: pi*5/8") +ax1.scatter(bin_center, dir_vario[0], label="emp. vario: pi/8") +ax1.scatter(bin_center, dir_vario[1], label="emp. vario: pi*5/8") ax1.legend(loc="lower right") model.plot("vario_axis", axis=0, ax=ax1, x_max=40, label="fit on axis 0") diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 9b981a05..023954f0 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -86,12 +86,12 @@ ax2.set_title("Tait-Bryan main axis") ax2.legend(loc="lower left") -ax3.plot(bin_center, dir_vario[0], label="0. axis") -ax3.plot(bin_center, dir_vario[1], label="1. axis") -ax3.plot(bin_center, dir_vario[2], label="2. axis") -model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0", linestyle=":") -model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1", linestyle=":") -model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2", linestyle=":") +ax3.scatter(bin_center, dir_vario[0], label="0. axis") +ax3.scatter(bin_center, dir_vario[1], label="1. axis") +ax3.scatter(bin_center, dir_vario[2], label="2. axis") +model.plot("vario_axis", axis=0, ax=ax3, label="fit on axis 0") +model.plot("vario_axis", axis=1, ax=ax3, label="fit on axis 1") +model.plot("vario_axis", axis=2, ax=ax3, label="fit on axis 2") ax3.set_title("Fitting an anisotropic model") ax3.legend() From f78f8c0292c6a2a9f30e48c7234f689d39fc53ee Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 22:18:23 +0100 Subject: [PATCH 19/28] Field: add meshio as dependency; update Field.mesh routine to meshio 4 --- gstools/field/base.py | 15 +++++++-------- requirements.txt | 1 + tests/test_srf.py | 19 +++++++++++++++++++ 3 files changed, 27 insertions(+), 8 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index a21da04b..7b3cfe7a 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -101,7 +101,7 @@ def mesh( See: :any:`Field.__call__` """ - select = _get_select(direction) + select = _get_select(direction)[: self.model.dim] if len(select) < self.model.dim: raise ValueError( "Field.mesh: need at least {} direction(s), got '{}'".format( @@ -117,12 +117,11 @@ def mesh( else: if points == "centroids": # define unique order of cells - cells = list(mesh.cells) offset = [] length = [] pnts = np.empty((0, 3), dtype=np.double) - for cell in cells: - pnt = np.mean(mesh.points[mesh.cells[cell]], axis=1) + for cell in mesh.cells: + pnt = np.mean(mesh.points[cell[1]], axis=1) offset.append(pnts.shape[0]) length.append(pnt.shape[0]) pnts = np.vstack((pnts, pnt)) @@ -134,10 +133,10 @@ def mesh( else: # if multiple values are returned, take the first one field = out[0] - field_dict = {} - for i, cell in enumerate(cells): - field_dict[cell] = field[offset[i] : offset[i] + length[i]] - mesh.cell_data[name] = field_dict + field_list = [] + for i in range(len(offset)): + field_list.append(field[offset[i] : offset[i] + length[i]]) + mesh.cell_data[name] = field_list else: out = self.unstructured(pos=mesh.points.T[select], **kwargs) if isinstance(out, np.ndarray): diff --git a/requirements.txt b/requirements.txt index 56c77216..0f47ca29 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ scipy>=1.1.0 hankel>=1.0.2 emcee>=3.0.0 pyevtk>=1.1.1 +meshio>=4.0.3, <5.0 diff --git a/tests/test_srf.py b/tests/test_srf.py index e5efa107..4f9a6f26 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -8,6 +8,7 @@ import numpy as np from gstools import SRF, Gaussian from gstools import transform as tf +import meshio class TestSRF(unittest.TestCase): @@ -313,6 +314,24 @@ def test_assertions(self): # mesh_type="hyper_mesh", # ) + def test_meshio(self): + points = np.array( + [ + [0.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0], + [1.0, 0.0, 0.0], + ] + ) + cells = [("tetra", np.array([[0, 1, 2, 3]]))] + mesh = meshio.Mesh(points, cells) + model = Gaussian(dim=3, len_scale=0.1) + srf = SRF(model) + srf.mesh(mesh, points="points") + self.assertEqual(len(srf.field), 4) + srf.mesh(mesh, points="centroids") + self.assertEqual(len(srf.field), 1) + if __name__ == "__main__": unittest.main() From 67fd5ce789b9f9ad490f81e49d3200bdc1893c39 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 19 Nov 2020 22:31:11 +0100 Subject: [PATCH 20/28] Field.mesh: better handling of directions --- gstools/field/base.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 7b3cfe7a..d04d067b 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -67,8 +67,8 @@ def unstructured(self, *args, **kwargs): return call(*args, **kwargs) def mesh( - self, mesh, points="centroids", direction="xyz", name="field", **kwargs - ): # pragma: no cover + self, mesh, points="centroids", direction="all", name="field", **kwargs + ): """Generate a field on a given meshio or ogs5py mesh. Parameters @@ -81,11 +81,12 @@ def mesh( (calculated as mean of the cell vertices) or the "points" of the given mesh. Default: "centroids" - direction : :class:`str`, optional + direction : :class:`str` or :class:`list`, optional Here you can state which direction should be choosen for lower dimension. For example, if you got a 2D mesh in xz direction, - you have to pass "xz" - Default: "xyz" + you have to pass "xz". By default, all directions are used. + One can also pass a list of indices. + Default: "all" name : :class:`str`, optional Name to store the field in the given mesh as point_data or cell_data. Default: "field" @@ -101,7 +102,12 @@ def mesh( See: :any:`Field.__call__` """ - select = _get_select(direction)[: self.model.dim] + if isinstance(direction, str) and direction == "all": + select = list(range(self.model.dim)) + elif isinstance(direction, str): + select = _get_select(direction)[: self.model.dim] + else: + select = direction[: self.model.dim] if len(select) < self.model.dim: raise ValueError( "Field.mesh: need at least {} direction(s), got '{}'".format( From 4a601cb107389241c37f77eab12aea307c5ea3c9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 20 Nov 2020 14:02:37 +0100 Subject: [PATCH 21/28] Field.mesh: check mesh dimension --- gstools/field/base.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index d04d067b..6ea836a3 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -125,7 +125,10 @@ def mesh( # define unique order of cells offset = [] length = [] - pnts = np.empty((0, 3), dtype=np.double) + mesh_dim = mesh.points.shape[1] + if mesh_dim < self.model.dim: + raise ValueError("Field.mesh: mesh dimension to low!") + pnts = np.empty((0, mesh_dim), dtype=np.double) for cell in mesh.cells: pnt = np.mean(mesh.points[cell[1]], axis=1) offset.append(pnts.shape[0]) From 05f7da482ae41ca93b4664d4f8da94fb017c21ef Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 20 Nov 2020 15:49:22 +0100 Subject: [PATCH 22/28] Field.mesh: check for meshio mesh --- gstools/field/base.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 6ea836a3..96c72501 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -14,6 +14,7 @@ from functools import partial import numpy as np +import meshio from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, gen_mesh @@ -120,7 +121,7 @@ def mesh( else: pnts = mesh.NODES.T[select] out = self.unstructured(pos=pnts, **kwargs) - else: + elif isinstance(mesh, meshio.Mesh): if points == "centroids": # define unique order of cells offset = [] @@ -154,6 +155,8 @@ def mesh( # if multiple values are returned, take the first one field = out[0] mesh.point_data[name] = field + else: + raise ValueError("Field.mesh: Unknown mesh format!") return out def pre_pos(self, pos, mesh_type="unstructured"): @@ -251,11 +254,7 @@ def _to_vtk_helper( filename, self.pos, {fieldname: field}, self.mesh_type ) else: - print( - "Field.to_vtk: No " - + field_select - + " stored in the class." - ) + print("Field.to_vtk: '{}' not available.".format(field_select)) else: raise ValueError( "Unknown field value type: {}".format(self.value_type) From 9b97f873877c7d759ab7f5f3b6a1d440a1350e90 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 20 Nov 2020 15:49:54 +0100 Subject: [PATCH 23/28] Examples: add an example with meshes --- examples/01_random_field/05_mesh_ensemble.py | 94 ++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100755 examples/01_random_field/05_mesh_ensemble.py diff --git a/examples/01_random_field/05_mesh_ensemble.py b/examples/01_random_field/05_mesh_ensemble.py new file mode 100755 index 00000000..76381729 --- /dev/null +++ b/examples/01_random_field/05_mesh_ensemble.py @@ -0,0 +1,94 @@ +""" +Generating fields on meshes +--------------------------- + +GSTools provides an interface for meshes, to support +`meshio `_ and +`ogs5py `_ meshes. + +When using `meshio`, the generated fields will be stored imediatly in the mesh +container. + +One has two options to generate a field on a given mesh: + +- `points="points"` will generate a field on the mesh points +- `points="centroids"` will generate a field on the cell centroids + +In this example, we will generate a simple mesh with the aid of +`meshzoo `_. +""" + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.tri as tri +import meshzoo +import meshio +import gstools as gs + +# generate a triangulated hexagon with meshzoo +points, cells = meshzoo.ngon(6, 4) +mesh = meshio.Mesh(points, {"triangle": cells}) + +############################################################################### +# Now we prepare the SRF class as always. We will generate an ensemble of +# fields on the generated mesh. + +# number of fields +fields_no = 12 +# model setup +model = gs.Gaussian(dim=2, len_scale=0.5) +srf = gs.SRF(model, mean=1) + +############################################################################### +# To generate fields on a mesh, we provide a separate method: :any:`SRF.mesh`. +# First we generate fields on the mesh-centroids controlled by a seed. +# You can specify the field name by the keyword `name`. + +for i in range(fields_no): + srf.mesh(mesh, points="centroids", name="c-field-{}".format(i), seed=i) + +############################################################################### +# Now we generate fields on the mesh-points again controlled by a seed. + +for i in range(fields_no): + srf.mesh(mesh, points="points", name="p-field-{}".format(i), seed=i) + +############################################################################### +# To get an impression we now want to plot the generated fields. +# Luckily, matplotlib supports triangular meshes. + + +triangulation = tri.Triangulation(points[:, 0], points[:, 1], cells) +# figure setup +cols = 4 +rows = int(np.ceil(fields_no / cols)) + +############################################################################### +# Cell data can be easily visualized with matplotlibs `tripcolor`. +# To highlight the cell structure, we use `triplot`. + +fig = plt.figure(figsize=[2 * cols, 2 * rows]) +for i, field in enumerate(mesh.cell_data, 1): + ax = fig.add_subplot(rows, cols, i) + ax.tripcolor(triangulation, mesh.cell_data[field][0]) + ax.triplot(triangulation, linewidth=0.5, color="k") + ax.set_aspect("equal") +fig.tight_layout() + +############################################################################### +# Point data is plotted via `tricontourf`. + +fig = plt.figure(figsize=[2 * cols, 2 * rows]) +for i, field in enumerate(mesh.point_data, 1): + ax = fig.add_subplot(rows, cols, i) + ax.tricontourf(triangulation, mesh.point_data[field]) + ax.triplot(triangulation, linewidth=0.5, color="k") + ax.set_aspect("equal") +fig.tight_layout() +plt.show() + +############################################################################### +# Last but not least, `meshio` can be used for what is does best: Exporting. +# Tada! + +mesh.write("mesh_ensemble.vtk") From e1dd8d422053d52802eef2f14dabec391304bc47 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 20 Nov 2020 15:50:11 +0100 Subject: [PATCH 24/28] Examples: add a 4D example --- .../01_random_field/06_higher_dimensions.py | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100755 examples/01_random_field/06_higher_dimensions.py diff --git a/examples/01_random_field/06_higher_dimensions.py b/examples/01_random_field/06_higher_dimensions.py new file mode 100755 index 00000000..b9c8a445 --- /dev/null +++ b/examples/01_random_field/06_higher_dimensions.py @@ -0,0 +1,75 @@ +""" +Higher Dimensions +----------------- + +GSTools provides experimental support for higher dimensions. + +Anisotropy is the same as in lower dimensions: + +- in `n` dimensions we need `(n-1)` anisotropy ratios + +Rotation on the other hand is a bit more complex. +With increasing dimensions more and more rotation angles are added in order +to properply describe the rotated axes of anisotropy. + +By design the first rotation angles coincide with the lower ones: + +- 2D (rotation in x-y plane) -> 3D: first angle describes xy-plane rotation +- 3D (tait-bryan angles) -> 4D: first 3 angles coincide with tait-bryan angles + +By increasing the dimension from `n` to `(n+1)`, `n` angles are added: + +- 2D (1 angle) -> 3D: 3 angles (2 added) +- 3D (3 angles) -> 4D: 6 angles (3 added) + +the following list of rotation-planes are described by the list of +angles in the model: + +1. x-y plane +2. x-z plane +3. y-z plane +4. x-v plane +5. y-v plane +6. z-v plane +7. ... + +The rotation direction in these planes have alternating signs +in order to match tait-bryan in 3D. + +Let's have a look at a 4D example, where we naively add a 4th dimension. +""" + +import matplotlib.pyplot as plt +import gstools as gs + +dim = 4 +size = 20 +pos = [range(size)] * dim +model = gs.Exponential(dim=dim, len_scale=5) +srf = gs.SRF(model, seed=20170519) +field = srf.structured(pos) + +############################################################################### +# In order to "prove" correctnes, we can calculate an empirical variogram +# of the generated field and fit our model to it. + +bin_edges = range(size) +bin_center, vario = gs.vario_estimate( + pos, field, bin_edges, sampling_size=2000, mesh_type="structured" +) +model.fit_variogram(bin_center, vario) +print(model) + +############################################################################### +# As you can see, the estimated variance and length scale match our input +# quite good. +# +# Let's have a look at the fit and a x-y cross-section of the 4D field: + +f, a = plt.subplots(1, 2, gridspec_kw={"width_ratios": [2, 1]}, figsize=[9, 3]) +model.plot(x_max=size + 1, ax=a[0]) +a[0].scatter(bin_center, vario) +a[1].imshow(field[:, :, 0, 0].T, origin="lower") +a[0].set_title("isotropic empirical variogram with fitted model") +a[1].set_title("x-y cross-section") +f.show() From 24dfc1df9900f09070515ad16519e988617aa2cd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 20 Nov 2020 15:53:07 +0100 Subject: [PATCH 25/28] Docs: add meshzoo as doc dependency --- docs/requirements_doc.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/requirements_doc.txt b/docs/requirements_doc.txt index c9d3ee24..cdaea2ce 100755 --- a/docs/requirements_doc.txt +++ b/docs/requirements_doc.txt @@ -3,3 +3,4 @@ sphinx-gallery matplotlib pyvista pykrige +meshzoo From de4c923bf7087dbc537d8f273206c3ddf6c758fb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 23 Nov 2020 16:06:39 +0100 Subject: [PATCH 26/28] Geometric: bugfix in 'format_unstruct_pos_shape': wrong dimension determined under some circumstances --- gstools/tools/geometric.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 8d062d2f..a16c7ee4 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -507,24 +507,25 @@ def format_unstruct_pos_shape(pos, shape, check_stacked_shape=False): stacked_shape_size = np.prod(shape[1:]) wrong_shape = False # now we try to be smart + pre_len = len(np.atleast_1d(pos)) + # care about 1D: pos can be given as 1D array here -> convert to 2D array pos = np.array(pos, dtype=np.double, ndmin=2) + post_len = len(pos) + # first array dimension should be spatial dimension (1D is special case) + dim = post_len if pre_len == post_len else 1 pnt_cnt = pos[0].size # case: 1D unstacked - if pos.size == shape_size: - dim = 1 + if dim == 1 and pos.size == shape_size: shape = (1, pos.size) if check_stacked_shape else (pos.size,) # case: 1D and stacked - elif pos.size == stacked_shape_size: - dim = 1 + elif dim == 1 and pos.size == stacked_shape_size: shape = (shape[0], pos.size) wrong_shape = not check_stacked_shape # case: nD unstacked elif pnt_cnt == shape_size: - dim = len(pos) shape = (1, pnt_cnt) if check_stacked_shape else pnt_cnt # case: nD and stacked elif pnt_cnt == stacked_shape_size: - dim = len(pos) shape = (shape[0], pnt_cnt) wrong_shape = not check_stacked_shape # if nothing works, we raise an error From 04c4c6c44bd3c2588d8cbb5537bff1fe71befb82 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Mon, 23 Nov 2020 18:15:02 +0100 Subject: [PATCH 27/28] Fix some typos --- gstools/field/base.py | 6 +++--- gstools/tools/geometric.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 96c72501..cf74c73f 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -25,7 +25,7 @@ class Field: - """A field base class for random and kriging fields ect. + """A base class for random fields, kriging fields, etc. Parameters ---------- @@ -128,7 +128,7 @@ def mesh( length = [] mesh_dim = mesh.points.shape[1] if mesh_dim < self.model.dim: - raise ValueError("Field.mesh: mesh dimension to low!") + raise ValueError("Field.mesh: mesh dimension too low!") pnts = np.empty((0, mesh_dim), dtype=np.double) for cell in mesh.cells: pnt = np.mean(mesh.points[cell[1]], axis=1) @@ -175,7 +175,7 @@ def pre_pos(self, pos, mesh_type="unstructured"): Returns ------- iso_pos : (d, n), :class:`numpy.ndarray` - the isometrize position tuple + the isometrized position tuple shape : :class:`tuple` Shape of the resulting field. """ diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index a16c7ee4..a02e13b0 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -386,7 +386,7 @@ def format_struct_pos_shape(pos, shape, check_stacked_shape=False): """ Format a structured position tuple with given shape. - Shape could be stacked, when multiple fields were given. + Shape could be stacked, when multiple fields are given. Parameters ---------- From 713843e1bf291e2f37b92589508fae2fe4491a98 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 23 Nov 2020 23:29:59 +0100 Subject: [PATCH 28/28] Randmeth: rename summators; remove superfluous tests for structured summation --- gstools/field/generator.py | 8 +++----- gstools/field/summator.pyx | 4 ++-- tests/test_incomprrandmeth.py | 31 ------------------------------- tests/test_randmeth.py | 22 ---------------------- 4 files changed, 5 insertions(+), 60 deletions(-) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 77858f68..213f3841 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -16,7 +16,7 @@ import numpy as np from gstools.covmodel.base import CovModel from gstools.random.rng import RNG -from gstools.field.summator import summate_unstruct, summate_incompr_unstruct +from gstools.field.summator import summate, summate_incompr __all__ = ["RandMeth", "IncomprRandMeth"] @@ -98,9 +98,7 @@ def __call__(self, pos): the random modes """ pos = np.array(pos, dtype=np.double) - summed_modes = summate_unstruct( - self._cov_sample, self._z_1, self._z_2, pos - ) + summed_modes = summate(self._cov_sample, self._z_1, self._z_2, pos) nugget = self._set_nugget(summed_modes.shape) return np.sqrt(self.model.var / self._mode_no) * summed_modes + nugget @@ -367,7 +365,7 @@ def __call__(self, pos): the random modes """ pos = np.array(pos, dtype=np.double) - summed_modes = summate_incompr_unstruct( + summed_modes = summate_incompr( self._cov_sample, self._z_1, self._z_2, pos ) nugget = self._set_nugget(summed_modes.shape) diff --git a/gstools/field/summator.pyx b/gstools/field/summator.pyx index c3da4459..80b5cf4c 100644 --- a/gstools/field/summator.pyx +++ b/gstools/field/summator.pyx @@ -16,7 +16,7 @@ DTYPE = np.double ctypedef np.double_t DTYPE_t -def summate_unstruct( +def summate( const double[:,:] cov_samples, const double[:] z_1, const double[:] z_2, @@ -52,7 +52,7 @@ cdef (double) abs_square(const double[:] vec) nogil: return r -def summate_incompr_unstruct( +def summate_incompr( const double[:,:] cov_samples, const double[:] z_1, const double[:] z_2, diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index 72167f79..1fd6f7ee 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -42,37 +42,6 @@ def test_unstruct_3d(self): self.assertAlmostEqual(modes[0, 1], 1.38687858) self.assertAlmostEqual(modes[1, 0], -0.27245271) - # def test_struct_2d(self): - # modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") - # self.assertAlmostEqual(modes[0, 0, 0], 0.50751115) - # self.assertAlmostEqual(modes[0, 1, 0], 0.69751927) - # self.assertAlmostEqual(modes[1, 1, 1], -0.19747468) - - # def test_struct_3d(self): - # modes = self.rm_3d( - # self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" - # ) - # self.assertAlmostEqual(modes[0, 0, 0, 0], 1.49469700) - # self.assertAlmostEqual(modes[1, 0, 1, 1], 0.12813365) - # self.assertAlmostEqual(modes[1, 1, 0, 1], 0.01443056) - # self.assertAlmostEqual(modes[1, 1, 1, 1], -0.12304040) - - # def test_struct_unstruct(self): - # x_grid = np.arange(0.0, 2.0, 1.0) - # y_grid = np.arange(0.0, 2.0, 1.0) - # x_tuple = np.array((0.0, 0.0, 1.0, 1.0)) - # y_tuple = np.array((0.0, 1.0, 0.0, 1.0)) - # unstr_modes = self.rm_2d(x_tuple, y_tuple, mesh_type="unstructured") - # str_modes = self.rm_2d(x_grid, y_grid, mesh_type="structured") - # for d in range(2): - # k = 0 - # for i in range(len(x_grid)): - # for j in range(len(y_grid)): - # self.assertAlmostEqual( - # str_modes[d, i, j], unstr_modes[d, k] - # ) - # k += 1 - def test_assertions(self): cov_model_1d = Gaussian(dim=1, var=1.5, len_scale=2.5) self.assertRaises(ValueError, IncomprRandMeth, cov_model_1d) diff --git a/tests/test_randmeth.py b/tests/test_randmeth.py index 376aaff4..422567b5 100644 --- a/tests/test_randmeth.py +++ b/tests/test_randmeth.py @@ -44,28 +44,6 @@ def test_unstruct_3d(self): self.assertAlmostEqual(modes[0], 0.55488481) self.assertAlmostEqual(modes[1], 1.18506639) - # def test_struct_1d(self): - # modes = self.rm_1d(self.x_grid, mesh_type="structured") - # self.assertAlmostEqual(modes[0], 3.19799030) - # self.assertAlmostEqual(modes[1], 2.34788923) - - # def test_struct_2d(self): - # modes = self.rm_2d(self.x_grid, self.y_grid, mesh_type="structured") - # self.assertAlmostEqual(modes[0, 0], 1.67318010) - # self.assertAlmostEqual(modes[1, 0], 1.54740003) - # self.assertAlmostEqual(modes[0, 1], 2.02106551) - # self.assertAlmostEqual(modes[1, 1], 1.86883255) - - # def test_struct_3d(self): - # modes = self.rm_3d( - # self.x_grid, self.y_grid, self.z_grid, mesh_type="structured" - # ) - # self.assertAlmostEqual(modes[0, 0, 0], 0.55488481) - # self.assertAlmostEqual(modes[0, 1, 0], 0.41858766) - # self.assertAlmostEqual(modes[1, 1, 0], 0.95133855) - # self.assertAlmostEqual(modes[0, 1, 1], 0.65475042) - # self.assertAlmostEqual(modes[1, 1, 1], 1.40915120) - def test_reset(self): modes = self.rm_2d((self.x_tuple, self.y_tuple)) self.assertAlmostEqual(modes[0], 1.67318010)