Skip to content

Commit

Permalink
Merge e7c2b83 into c1b56da
Browse files Browse the repository at this point in the history
  • Loading branch information
MuellerSeb authored Jan 14, 2021
2 parents c1b56da + e7c2b83 commit 44befd8
Show file tree
Hide file tree
Showing 16 changed files with 394 additions and 322 deletions.
28 changes: 17 additions & 11 deletions examples/06_conditioned_fields/00_condition_ensemble.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,9 @@
Conditioning with Ordinary Kriging
----------------------------------
Here we use ordinary kriging in 1D (for plotting reasons) with 5 given observations/conditions,
Here we use ordinary kriging in 1D (for plotting reasons)
with 5 given observations/conditions,
to generate an ensemble of conditioned random fields.
The estimated mean can be accessed by ``srf.mean``.
"""
import numpy as np
import matplotlib.pyplot as plt
Expand All @@ -16,30 +16,35 @@
gridx = np.linspace(0.0, 15.0, 151)

###############################################################################
# The conditioned spatial random field class depends on a Krige class in order
# to handle the conditions.
# This is created as described in the kriging tutorial.
#
# Here we use a gaussian covariance model and ordinary kriging for conditioning
# the spatial random field.

# spatial random field class
model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5)
srf = gs.SRF(model)
srf.set_condition(cond_pos, cond_val, "ordinary")
krige = gs.krige.Ordinary(model, cond_pos, cond_val)
csrf = gs.CondSRF(krige)

###############################################################################

fields = []
for i in range(100):
# print(i) if i % 10 == 0 else None
fields.append(srf(gridx, seed=i))
fields.append(csrf(gridx, seed=i))
label = "Conditioned ensemble" if i == 0 else None
plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label)
plt.plot(gridx, np.full_like(gridx, srf.mean), label="estimated mean")
plt.plot(gridx, csrf.krige(gridx, only_mean=True), label="estimated mean")
plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean")
plt.plot(gridx, srf.krige_field, linestyle="dashed", label="kriged field")
plt.plot(gridx, csrf.krige.field, linestyle="dashed", label="kriged field")
plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions")
# 99 percent confidence interval
conf = gs.tools.confidence_scaling(0.99)
plt.fill_between(
gridx,
srf.krige_field - conf * np.sqrt(srf.krige_var),
srf.krige_field + conf * np.sqrt(srf.krige_var),
csrf.krige.field - conf * np.sqrt(csrf.krige.krige_var),
csrf.krige.field + conf * np.sqrt(csrf.krige.krige_var),
alpha=0.3,
label="99% confidence interval",
)
Expand All @@ -48,4 +53,5 @@

###############################################################################
# As you can see, the kriging field coincides with the ensemble mean of the
# conditioned random fields and the estimated mean is the mean of the far-field.
# conditioned random fields and the estimated mean
# is the mean of the far-field.
6 changes: 4 additions & 2 deletions gstools/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,12 +27,13 @@
Spatial Random Field
^^^^^^^^^^^^^^^^^^^^
Class for random field generation
Classes for (conditioned) random field generation
.. currentmodule:: gstools.field
.. autosummary::
SRF
CondSRF
Covariance Base-Class
^^^^^^^^^^^^^^^^^^^^^
Expand Down Expand Up @@ -113,7 +114,7 @@
"""
# Hooray!
from gstools import field, variogram, random, covmodel, tools, krige, transform
from gstools.field import SRF
from gstools.field import SRF, CondSRF
from gstools.tools import (
rotated_main_axes,
EARTH_RADIUS,
Expand Down Expand Up @@ -188,6 +189,7 @@

__all__ += [
"SRF",
"CondSRF",
"rotated_main_axes",
"EARTH_RADIUS",
"vtk_export",
Expand Down
4 changes: 0 additions & 4 deletions gstools/covmodel/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -1135,10 +1135,6 @@ def __ne__(self, other):
"""Compare CovModels."""
return not self.__eq__(other)

def __str__(self):
"""Return String representation."""
return self.__repr__()

def __repr__(self):
"""Return String representation."""
return model_repr(self)
9 changes: 3 additions & 6 deletions gstools/covmodel/plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -72,8 +72,7 @@ def plot_vario_spatial(
model, x_min=0.0, x_max=None, fig=None, ax=None
): # pragma: no cover
"""Plot spatial variogram of a given CovModel."""
field = gstools.field.base.Field(model)
field._value_type = "scalar"
field = gstools.field.base.Field(model, "scalar")
if x_max is None:
x_max = 3 * model.len_scale
x_s = np.linspace(-x_max, x_max) + x_min
Expand All @@ -88,8 +87,7 @@ def plot_cov_spatial(
model, x_min=0.0, x_max=None, fig=None, ax=None
): # pragma: no cover
"""Plot spatial covariance of a given CovModel."""
field = gstools.field.base.Field(model)
field._value_type = "scalar"
field = gstools.field.base.Field(model, "scalar")
if x_max is None:
x_max = 3 * model.len_scale
x_s = np.linspace(-x_max, x_max) + x_min
Expand All @@ -102,8 +100,7 @@ def plot_cor_spatial(
model, x_min=0.0, x_max=None, fig=None, ax=None
): # pragma: no cover
"""Plot spatial correlation of a given CovModel."""
field = gstools.field.base.Field(model)
field._value_type = "scalar"
field = gstools.field.base.Field(model, "scalar")
if x_max is None:
x_max = 3 * model.len_scale
x_s = np.linspace(-x_max, x_max) + x_min
Expand Down
4 changes: 3 additions & 1 deletion gstools/field/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,12 @@
.. autosummary::
SRF
CondSRF
----
"""

from gstools.field.srf import SRF
from gstools.field.cond_srf import CondSRF

__all__ = ["SRF"]
__all__ = ["SRF", "CondSRF"]
31 changes: 11 additions & 20 deletions gstools/field/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@

__all__ = ["Field"]

VALUE_TYPES = ["scalar", "vector"]


class Field:
"""A base class for random fields, kriging fields, etc.
Expand All @@ -31,21 +33,18 @@ class Field:
----------
model : :any:`CovModel`
Covariance Model related to the field.
mean : :class:`float`, optional
Mean value of the field.
"""

def __init__(self, model, mean=0.0):
def __init__(self, model, value_type="scalar"):
# initialize attributes
self.pos = None
self.mesh_type = None
self.field = None
# initialize private attributes
self._mean = None
self._model = None
self.mean = mean
self.model = model
self._value_type = None
self.value_type = value_type

def __call__(*args, **kwargs):
"""Generate the field."""
Expand Down Expand Up @@ -376,15 +375,6 @@ def plot(self, field="field", fig=None, ax=None): # pragma: no cover

return r

@property
def mean(self):
""":class:`float`: The mean of the field."""
return self._mean

@mean.setter
def mean(self, mean):
self._mean = float(mean)

@property
def model(self):
""":any:`CovModel`: The covariance model of the field."""
Expand All @@ -404,15 +394,16 @@ def value_type(self):
""":class:`str`: Type of the field values (scalar, vector)."""
return self._value_type

def __str__(self):
"""Return String representation."""
return self.__repr__()
@value_type.setter
def value_type(self, value_type):
""":class:`str`: Type of the field values (scalar, vector)."""
if value_type not in VALUE_TYPES:
raise ValueError("Field: value type not in {}".format(VALUE_TYPES))
self._value_type = value_type

def __repr__(self):
"""Return String representation."""
return "Field(model={0}, mean={1:.{p}})".format(
self.model, self.mean, p=self.model._prec
)
return "Field(model={0})".format(self.model)


def _names(name, cnt):
Expand Down
168 changes: 168 additions & 0 deletions gstools/field/cond_srf.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,168 @@
# -*- coding: utf-8 -*-
"""
GStools subpackage providing a class for conditioned spatial random fields.
.. currentmodule:: gstools.field.cond_srf
The following classes are provided
.. autosummary::
CondSRF
"""
# pylint: disable=C0103

import numpy as np
from gstools.field.generator import RandMeth
from gstools.field.base import Field
from gstools.krige import Krige

__all__ = ["CondSRF"]

GENERATOR = {
"RandMeth": RandMeth,
}


class CondSRF(Field):
"""A class to generate conditioned spatial random fields (SRF).
Parameters
----------
model : :any:`CovModel`
Covariance Model of the spatial random field.
generator : :class:`str`, optional
Name of the field generator to be used.
At the moment, the following generators are provided:
* "RandMeth" : The Randomization Method.
See: :any:`RandMeth`
* "IncomprRandMeth" : The incompressible Randomization Method.
This is the original algorithm proposed by Kraichnan 1970
See: :any:`IncomprRandMeth`
* "VectorField" : an alias for "IncomprRandMeth"
* "VelocityField" : an alias for "IncomprRandMeth"
Default: "RandMeth"
**generator_kwargs
Keyword arguments that are forwarded to the generator in use.
Have a look at the provided generators for further information.
"""

def __init__(self, krige, generator="RandMeth", **generator_kwargs):
if not isinstance(krige, Krige):
raise ValueError("CondSRF: krige should be an instance of Krige.")
self.krige = krige
# initialize attributes
self.pos = None
self.mesh_type = None
self.field = None
self.raw_field = None
# initialize private attributes
self._value_type = None
self._generator = None
# initialize attributes
self.set_generator(generator, **generator_kwargs)

def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs):
"""Generate the conditioned spatial random field.
The field is saved as `self.field` and is also returned.
Parameters
----------
pos : :class:`list`
the position tuple, containing main direction and transversal
directions
seed : :class:`int`, optional
seed for RNG for reseting. Default: keep seed from generator
mesh_type : :class:`str`
'structured' / 'unstructured'
**kwargs
keyword arguments that are forwarded to the kriging routine in use.
Returns
-------
field : :class:`numpy.ndarray`
the conditioned SRF
"""
kwargs["mesh_type"] = mesh_type
kwargs["only_mean"] = False # overwrite if given
kwargs["return_var"] = True # overwrite if given
# update the model/seed in the generator if any changes were made
self.generator.update(self.model, seed)
# get isometrized positions and the resulting field-shape
iso_pos, shape = self.pre_pos(pos, mesh_type)
# generate the field
self.raw_field = np.reshape(self.generator(iso_pos), shape)
field, krige_var = self.krige(pos, **kwargs)
var_scale, nugget = self.get_scaling(krige_var, shape)
self.field = field + var_scale * self.raw_field + nugget
return self.field

def get_scaling(self, krige_var, shape):
"""
Get scaling coefficients for the random field.
Parameters
----------
krige_var : :class:`numpy.ndarray`
Kriging variance.
shape : :class:`tuple` of :class:`int`
Field shape.
Returns
-------
var_scale : :class:`numpy.ndarray`
Variance scaling factor for the random field.
nugget : :class:`numpy.ndarray` or `class:`int
Nugget to be added to the field.
"""
if self.model.nugget > 0:
var_scale = np.maximum(krige_var - self.model.nugget, 0)
nug_scale = np.sqrt((krige_var - var_scale) / self.model.nugget)
var_scale = np.sqrt(var_scale / self.model.var)
nugget = nug_scale * self.generator.get_nugget(shape)
else:
var_scale = np.sqrt(krige_var / self.model.var)
nugget = 0
return var_scale, nugget

def set_generator(self, generator, **generator_kwargs):
"""Set the generator for the field.
Parameters
----------
generator : :class:`str`, optional
Name of the generator to use for field generation.
Default: "RandMeth"
**generator_kwargs
keyword arguments that are forwarded to the generator in use.
"""
if generator in GENERATOR:
gen = GENERATOR[generator]
self._generator = gen(self.model, **generator_kwargs)
self.value_type = self._generator.value_type
else:
raise ValueError("gstools.CondSRF: Unknown generator " + generator)
if self.value_type != "scalar":
raise ValueError("CondSRF: only scalar field value type allowed.")

@property
def generator(self):
""":any:`callable`: The generator of the field."""
return self._generator

@property
def model(self):
""":any:`CovModel`: The covariance model of the field."""
return self.krige.model

@model.setter
def model(self, model):
pass

def __repr__(self):
"""Return String representation."""
return "CondSRF(krige={0}, generator={1})".format(
self.krige, self.generator
)
Loading

0 comments on commit 44befd8

Please sign in to comment.