diff --git a/README.md b/README.md index 0309bbf..485f0bb 100644 --- a/README.md +++ b/README.md @@ -75,7 +75,7 @@ plot2D([np.array([x.ravel(), y.ravel()]).T], [rf.random_field], title="Random Fi The result is the following random field: -![Random field 2D](./docs/static/random_field.png) +![Random field 2D](./docs/_static/random_field.png) ## Example in 3D @@ -127,4 +127,4 @@ plot3D([np.array([x.ravel(), y.ravel(), z.ravel()]).T], [rf.random_field], title ``` The result is the following random field: -![Random field 3D](./docs/static/random_field_3D.png) +![Random field 3D](./docs/_static/random_field_3D.png) diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/_static/conditioned_random_field_3D.png b/docs/_static/conditioned_random_field_3D.png new file mode 100644 index 0000000..72ce4eb Binary files /dev/null and b/docs/_static/conditioned_random_field_3D.png differ diff --git a/docs/_static/kriging_mean_2D.png b/docs/_static/kriging_mean_2D.png new file mode 100644 index 0000000..8025dc1 Binary files /dev/null and b/docs/_static/kriging_mean_2D.png differ diff --git a/docs/_static/kriging_mean_3D.png b/docs/_static/kriging_mean_3D.png new file mode 100644 index 0000000..cd2fb49 Binary files /dev/null and b/docs/_static/kriging_mean_3D.png differ diff --git a/docs/_static/kriging_std_2D.png b/docs/_static/kriging_std_2D.png new file mode 100644 index 0000000..3a8774b Binary files /dev/null and b/docs/_static/kriging_std_2D.png differ diff --git a/docs/_static/kriging_std_3D.png b/docs/_static/kriging_std_3D.png new file mode 100644 index 0000000..4f95cc7 Binary files /dev/null and b/docs/_static/kriging_std_3D.png differ diff --git a/docs/static/random_field.png b/docs/_static/random_field.png similarity index 100% rename from docs/static/random_field.png rename to docs/_static/random_field.png diff --git a/docs/_static/random_field_1.png b/docs/_static/random_field_1.png new file mode 100644 index 0000000..0538fc6 Binary files /dev/null and b/docs/_static/random_field_1.png differ diff --git a/docs/_static/random_field_2.png b/docs/_static/random_field_2.png new file mode 100644 index 0000000..760b9e9 Binary files /dev/null and b/docs/_static/random_field_2.png differ diff --git a/docs/_static/random_field_3.png b/docs/_static/random_field_3.png new file mode 100644 index 0000000..1914f98 Binary files /dev/null and b/docs/_static/random_field_3.png differ diff --git a/docs/static/random_field_3D.png b/docs/_static/random_field_3D.png similarity index 100% rename from docs/static/random_field_3D.png rename to docs/_static/random_field_3D.png diff --git a/docs/_static/random_field_4.png b/docs/_static/random_field_4.png new file mode 100644 index 0000000..f1fc064 Binary files /dev/null and b/docs/_static/random_field_4.png differ diff --git a/docs/_static/random_field_5.png b/docs/_static/random_field_5.png new file mode 100644 index 0000000..0da8be0 Binary files /dev/null and b/docs/_static/random_field_5.png differ diff --git a/docs/authors.rst b/docs/authors.rst new file mode 100644 index 0000000..cf6cb99 --- /dev/null +++ b/docs/authors.rst @@ -0,0 +1,19 @@ +.. _authors: + +STEM team +========= + +Maintainers +----------- +- Bruno Zuada Coelho `@brunozc `_ +- Aron Noordam `@aronnoordam `_ + +Developers +---------- +- Marjan Fathian `@indigocoral `_ +- Jonathan Nuttall `@mcgicjn2 `_ +- Davide Moretti `@morettid `_ +- Jelme Pennings `@jkpennings `_ +- Bram van den Eijnden `@bramvandeneijnden `_ + + diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..60f76a3 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,62 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +import random_fields + + +project = 'RandomFields' +copyright = '2024, STEM team' +author = 'STEM team' +version = random_fields.__version__ + + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = ['sphinx.ext.autodoc', 'sphinx_rtd_theme', 'sphinx.ext.intersphinx'] + +templates_path = ['_templates'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + +# to add __init__ documentation to the build +autoclass_content = 'both' + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = 'sphinx_rtd_theme' +html_static_path = ['_static'] + +# # Configure the RTD theme options +html_theme_options = { + 'logo_only': False, + 'display_version': True, + 'prev_next_buttons_location': 'bottom', + 'style_external_links': True, # Toc options + 'collapse_navigation': False, + 'sticky_navigation': True, + 'navigation_depth': 2, + 'includehidden': True, + 'titles_only': True, +} + +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. +html_show_sphinx = False +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. +html_show_copyright = True +# If true shows the source link to the rst code +html_show_sourcelink = False + +# Link to the github +html_context = { + "display_github": True, + "github_user": "stemVibrations", + "github_repo": "RanndomFields", + "github_version": "main", + "conf_py_path": "/docs/", +} diff --git a/docs/contributions.rst b/docs/contributions.rst new file mode 100644 index 0000000..90b69c2 --- /dev/null +++ b/docs/contributions.rst @@ -0,0 +1,111 @@ +Code contribution +================= + +Steps for submitting your code +------------------------------ + +When contributing code follow this checklist: + + #. Fork the repository on GitHub. + #. Create an issue with the desired feature or bug fix. + #. Make your modifications or additions in a feature branch. + #. Make changes and commit your changes using a descriptive commit message. + #. Provide tests for your changes, and ensure they all pass. + #. Provide documentation for your changes, in accordance with the style of the rest of the project (see :ref:`style_guide`). + #. Create a pull request to STEM main branch. The STEM team will review and discuss your Pull Request with you. + +For any questions, please get in contact with one of the members of :doc:`authors`. + + +.. _style_guide: + +Code style guide +---------------- +The additional features should follow the style of the STEM project. + +We follow the PEP 8 style guide for Python code, with our custom modifications as defined in the +`Yapf file <../../.style.yapf>`_ and the `flake8 file <../../.flake8>`_. These files can be ran manually by using the +following command from the root directory of the project: + +.. code-block:: + + pre-commit run --all-files + + +The class or function name should be clear and descriptive of the functionality it provides. + +There should be a docstring at the beginning of the class or function describing its purpose and usage. +The docstring should be in the form of a triple-quoted string. + +The class or function must have a type annotation. +The class should specify the attributes and inheritance. +The function should specify the arguments, exceptions and returns (in this order). +The return type annotation should be in the form of a comment after the closing parenthesis of the arguments. + +Please, avoid inheritance, and favour composition when writing your code. + +An example of a class: + +.. code-block:: + + class ClassName(object): + """ + Work in progress: docstring style to be specified + + Inheritance: + -:class:`object` + + Attributes: + - (): ,shape (n,m) + - (): ,shape (m,) + + + """ + def __init__(self,,...): + """ + Constructor of the RandomFields class + + Args: + - : () , shape (n,m,), default + """ + self.attribute_name_1: = + self.attribute_name_2: = + + +An example of a function: + +.. code-block:: + + def generate(self, nodes: npt.NDArray[np.float64]) -> None: + """ + Generate random field + + Args: + - nodes (ndarray): The nodes of the random field. shape (:,`self.n_dim`) + + Raises: + - ValueError: if dimensions of `nodes` do not match dimensions of the model + + Returns: + + """ + + # check dimensions of nodes agrees with dimensions of model + if nodes.shape[1] != self.n_dim: + raise ValueError(f'Dimensions of nodes: {nodes.shape[1]} do not match dimensions of model: {self.n_dim}') + + # scale of fluctuation + scale_fluctuation = np.ones(self.n_dim) * self.vertical_scale_fluctuation + + # apply the anisotropy to the other dimensions + mask = np.arange(len(scale_fluctuation)) != self.v_dim + scale_fluctuation[mask] = scale_fluctuation[mask] * self.anisotropy + + model = self.random_field_model(dim = self.n_dim, + var = self.variance, + len_scale = scale_fluctuation, + angles = self.angle) + self.random_field = gs.SRF(model, + mean = self.mean, + seed = self.seed) + self.random_field(nodes.T) diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..142a41a --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,54 @@ +RandomFields +==== +Python package to generate random fields in 2D or 3D on unstructured grids. + + +User guide +========== + +Background information +...................... +The random fields are generated using [gstools](https://geostat-framework.readthedocs.io/en/latest/). + +Currently, the following covariance models are available: + +* Gaussian +* Exponential +* Matern +* Linear (not for conditioned random fields) + +Additional functoinalities are used to generate conditioned random fields. See :doc:`theory` for details. + + +Tutorials +......... + +* :ref:`tutorial1` + 1D unconditioned and conditioned random field generation. + +* :ref:`tutorial2` + 2D and 3D unconditioned and conditioned random field simulation + +* :ref:`tutorial3` + Non-normal Gaussian fields and conditioning based on CPT data + +How to contribute +================= +If you want to contribute to RandomFields please follow the steps defined in :doc:`contributions`. + + +STEM team +========= +STEM is a research programme that results from a collaboration between the following partners: + +* `ProRail `_ +* `Deltares `_ +* `TNO `_ +* `TU Delft `_ + +The :doc:`authors`. + +Package documentation +===================== + +The :doc:`random_fields` documentation. diff --git a/docs/installation.rst b/docs/installation.rst new file mode 100644 index 0000000..c9dcb89 --- /dev/null +++ b/docs/installation.rst @@ -0,0 +1,10 @@ +RandomFields installation +================= + +.. _python_random_fields: + +RandomFields Python package +------------------- + + +To be copied from STEM \ No newline at end of file diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..954237b --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..2540854 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,2 @@ +Sphinx==6.2.1 +sphinx-rtd-theme==1.2.2 diff --git a/docs/theory.rst b/docs/theory.rst new file mode 100644 index 0000000..7dc7389 --- /dev/null +++ b/docs/theory.rst @@ -0,0 +1,6 @@ + +.. _theory: + +Theory +====== + diff --git a/docs/tutorials.rst b/docs/tutorials.rst new file mode 100644 index 0000000..23faeae --- /dev/null +++ b/docs/tutorials.rst @@ -0,0 +1,252 @@ +RandomFields tutorials +====================== + +.. _tutorial1: + +Example in 1D +------------- +This tutorial shows step by step guide on how to set up the random field generator and visualise the resulting random fields. + +First the necessary packages are imported and paths are defined. + +.. code-block:: python + + import numpy as np + from random_fields.generate_field import RandomFields, ModelName + from random_fields.utils import plot2D + +Then, you need to define the 1D grid: + +.. code-block:: python + + x = np.linspace(0, 100, 501).reshape(-1,1) + +In this example, we create a random field over a line of length 100, divided into steps of 0.2 + +Then we define the random field properties: + +.. code-block:: python + + nb_dimensions = 1 + mean = 10 + variance = 2 + vertical_scale_fluctuation = 10 + anisotropy = [] + angle = [] + model_rf = ModelName.Gaussian + +In this example, we use a a Gaussian covariance model. + +Then we create the random field, in which the index of the 'vertical' dimension needs to be set to `0`: + +.. code-block:: python + + rf = RandomFields(model_rf, nb_dimensions, mean, variance, + vertical_scale_fluctuation, anisotropy, angle, seed=1) + rf.v_dim = 0 + rf.generate(x.reshape(-1,1)) + +To visualise the results, you can run: + +.. code-block:: python + + plot1D([x],[rf.random_field], + title = 'Random Field', + output_folder="./", + output_name="random_field_1.png") + +.. image:: _static/random_field_1.png + + +Conditioning data can be added to the random field. The noise level represents the uncertainty in the individual conditioning points and is formulated as the normalised variance: the variance of the uncertainty in a single data point divided by the variance of the random field. ``noise_level`` can be a scalar, or an array of the same dimensions as ``conditioning _values``. + +.. code-block:: python + + conditioning_points = np.array([[30,50,70]]).T + conditioning_values = np.array([-2.5,0,2.5]) + + rf.set_conditioning_points(conditioning_points,conditioning_values, noise_level = 0.04) + + +Ten conditional random fields can be sampled and visualised as follows: + +.. code-block:: python + + for i in range(10): + rf.seed = i + rf.generate_conditioned(x) + rf_conditioned[i] = rf.conditioned_random_field + + plot1D([x]*10,rf_conditioned, + title = 'Conditioned Random Field', + output_folder="./", + output_name="random_field_2.png") + + +.. image:: _static/random_field_2.png + + +As a by-product of the conditioned random field generator, the kriging mean and variance are provided: + +.. code-block:: python + + mean = rf.kriging_mean + std = rf.kriging_std + + plt.scatter(conditioning_points, conditioning_values, label = 'conditioning points') + plt.plot(x, mean, label = 'kriged field') + plt.fill_between(x.ravel(), mean - 1.65*std, mean + 1.65*std, label = '90% uncertainty bound') + + +.. image:: _static/random_field_3.png + + +Note that only simple kriging is applied for conditioning the random fields. See :ref:`theory` for more details. + + + +.. _tutorial2: + +Example in 2D +------------- + +A (conditioned) 2D random field is constructed as follows: + +.. code-block:: python + + import numpy as np + from random_fields.generate_field import RandomFields, ModelName + from random_fields.utils import plot2D + + + x = np.linspace(0, 100, 51) + y = np.linspace(0, 50, 51) + x, y = np.meshgrid(x, y) + + nb_dimensions = 2 + mean = 10 + variance = 2 + vertical_scale_fluctuation = 10 + anisotropy = [1] + angle = [0] + model_rf = ModelName.Gaussian + +An unconditioned random field is constructed as: + +.. code-block:: python + + rf = RandomFields(model_rf, nb_dimensions, mean, variance, vertical_scale_fluctuation, anisotropy, angle, seed=14) + rf.generate(np.array([x.ravel(), y.ravel()]).T) + + plot2D([np.array([x.ravel(), y.ravel()]).T], [rf.random_field], title="Random Field", output_folder="./", output_name="random_field_4.png") + + +.. image:: _static/random_field_4.png + +Conditioning is consistent with the 1D example: + +.. code-block:: python + + xc = np.array([50.]*51) + yc = np.linspace(0,50,51) + + vc = np.array([15]*51) + + rf.set_conditioning_points(np.array([xc,yc]).T,vc,noise_level = 0.0001) + rf.generate_conditioned(np.array([x.ravel(), y.ravel()]).T) + + plot2D([np.array([x.ravel(), y.ravel()]).T], [rf.conditioned_random_field], + title="Conditioned Random Field", + output_folder="./", + output_name="random_field_5.png") + + +.. image:: _static/random_field_5.png + +The kriging mean and variance can be visualised as: + +.. code-block:: python + + plot2D([np.array([x.ravel(), y.ravel()]).T], [rf.kriging_mean], + title="kriging mean", + output_folder="./", + output_name="kriging_mean_2D.png") + + plot2D([np.array([x.ravel(), y.ravel()]).T], [rf.kriging_std], + title="kriging std", + output_folder="./", + output_name="kriging_std_2D.png") + +.. image:: _static/kriging_mean_2D.png + :width: 49% +.. image:: _static/kriging_std_2D.png + :width: 49% + +In the same way, a 3D conditional random field can be created: + + +.. code-block:: python + + # mesh coordinates + x = np.linspace(0, 100, 51) + y = np.linspace(0, 50, 51) + z = np.linspace(0, 25, 51) + x, y, z = [i.ravel() for i in np.meshgrid(x, y, z)] + + # random field properties + nb_dimensions = 3 + mean = 10 + variance = 2 + vertical_scale_fluctuation = 10 + anisotropy = [2.5, 2.5] + angle = [0, 0] + model_rf = ModelName.Gaussian + + # generate and plot random field + rf = RandomFields(model_rf, nb_dimensions, mean, variance, vertical_scale_fluctuation, + anisotropy, angle, seed=14) + rf.generate(np.array([x, y, z]).T) + plot3D([np.array([x, y, z]).T], [rf.random_field], + title="Random Field", + output_folder="./", + output_name="random_field_3D.png") + + # declae conditioning points + xc = np.array([50.]*51) + yc = np.linspace(0,50,51) + zc = np.array([25]*51) + vc = np.array([15]*51) + rf.set_conditioning_points(np.array([xc,yc,zc]).T,vc,noise_level = 0.0001) + + # generate and plot conditioned random field model + rf.generate_conditioned(np.array([x, y,z]).T) + plot3D([np.array([x, y, z]).T], [rf.conditioned_random_field], + title="Conditioned Random Field", + output_folder="./", + output_name="conditioned_random_field_3D.png") + plot3D([np.array([x, y, z]).T], [rf.kriging_mean], + title="Kriging mean 3D", + output_folder="./", + output_name="kriging_mean_3D.png") + plot3D([np.array([x, y, z]).T], [rf.kriging_std], + title="Kriging std 3D", + output_folder="./", + output_name="kriging_std_3D.png") + + +.. image:: _static/conditioned_random_field_3D.png + :width: 49% +.. image:: _static/kriging_mean_3D.png + :width: 49% +.. image:: _static/kriging_std_3D.png + :width: 49% + + + + +.. _tutorial3: + +Tutorial RandomFields in STEM +----------------------------- + +TODO \ No newline at end of file diff --git a/random_fields/generate_field.py b/random_fields/generate_field.py index 3ec494e..df0a657 100644 --- a/random_fields/generate_field.py +++ b/random_fields/generate_field.py @@ -1,9 +1,15 @@ from enum import Enum from typing import List import numpy as np +import numpy.typing as npt import gstools as gs +from sklearn.gaussian_process import GaussianProcessRegressor as GPR +from sklearn.gaussian_process.kernels import RBF, ConstantKernel , \ + WhiteKernel , Matern + + class ModelName(Enum): """ Name of the model to be used. Options are: "Gaussian", "Exponential", "Matern", "Linear" @@ -14,38 +20,71 @@ class ModelName(Enum): Linear = gs.Linear -class RandomFields: +class RandomFields(): """ - Generate random fields + Generates (conditioned) random fields + + Inheritance: + - None + + Attributes: + - max_conditioning_points (int): maximum number of + - random_field_model_name (str): + - random_field_model (int): + - n_dim (int): number of physical dimesions (1,2 or 3) + - seed (int): seed for the random number generator + - mean (float): unconditioned mean of the random field + - variance (float): unconditioned variance of the random field + - vertical_scale_fluctuation (float): + - anisotropy (list): ratio between horizontal scales and vertical scale of fluctuation, shape (`self.n_dims` - 1) + - angle (list): rotation angles of the principal directions of the scales of fluctuation relative to the vertical, shape (`self.n_dims` - 1) + - v_dim (int): dimesion number corresponding to the vertical scale of fluctuation + - random_field (:class:gstools.SRF) random field generator + - z_kriged_field (array): standard-normal kriged random field + - gaussian_process (:class:sklearn.gaussian_process.GaussianProcessRegressor) + - conditioning_points (array): coordinates of the conditioning points, shape (:,`self.n_dims`) + - conditioning_values (array): values of the conditioning points, shape (:,) + - kriging_mean (array): mean of the conditioned random field, shape (:,) + - kriging_std = (array): standard deviation of the conditioned random field, shape (:,) + - conditioned_random_field (array): conditioned random field, shape (:,) + """ - def __init__(self, model_name: ModelName, n_dim: int, - mean: float, variance: float, - v_scale_fluctuation: float, anisotropy: List[float], angle: List[float], - seed: int = 14, v_dim: int = 1) -> None: + def __init__(self, + model_name: ModelName, + n_dim: int, + mean: float, + variance: float, + v_scale_fluctuation: float, + anisotropy: List[float], + angle: List[float], + seed: int = 14, + v_dim: int = 1, + max_conditioning_points:int = 2000) -> None: """ Initialise generation of random fields - Parameters: - ----------- - model_name: str - Name of the model to be used. Options are: "Gaussian", "Exponential", "Matern", "Linear" - n_dim: int - The dimension of the random field - mean: float - The mean of the random field - variance: float - The variance of the random field - v_scale_fluctuation: float - The vertical scale of fluctuation of the random field - anisotropy: list - The anisotropy of the random field (per dimension) - angle: list - The angle of the random field (per dimension) - seed: int - The seed number for the random number generator - v_dim: int - The dimension of the vertical scale of fluctuation + Args: + - model_name (str): Name of the model to be used. Options are: "Gaussian", "Exponential", "Matern", "Linear" + - n_dim (int): The dimension of the random field + - mean (float): The mean of the random field + - variance (float): The variance of the random field + - v_scale_fluctuation (float): The vertical scale of fluctuation of the random field + - anisotropy (list): The anisotropy of the random field (per dimension) + - angle (list): The angle of the random field (per dimension) + - seed (int): The seed number for the random number generator + - v_dim (int): The dimension of the vertical scale of fluctuation + - max_conditioning_points (int): Maximum number of points to be used as conditioning points: increase leads to + longer computation times and high memory use due to factorisation of dense matrices + + Raises: + - ValueError: if model_name.name if not a member of class ModelName + - ValueError: if the number of dimensions is not in [1,2,3] + - ValueError: if the number of dimensions is nt consistent with the number of random field properties + + Returns: + - self + """ # initialise model if model_name.name not in ModelName.__members__: @@ -57,9 +96,12 @@ def __init__(self, model_name: ModelName, n_dim: int, aux_keys = [anisotropy, angle] for key in aux_keys: if len(key) != n_dim - 1: - raise ValueError(f'Number of dimensions: {n_dim} does not match number of\ - random field properties: {len(variance)}. It should be {n_dim - 1}.') + raise ValueError(f'Number of dimensions: {n_dim} does not match number of ' + + f'random field properties: {len(key)}. It should be {n_dim - 1}.') + + self.max_conditioning_points = max_conditioning_points + self.random_field_model_name = model_name.name self.random_field_model = model_name.value self.n_dim = n_dim self.seed = seed @@ -70,15 +112,27 @@ def __init__(self, model_name: ModelName, n_dim: int, self.angle = angle self.v_dim = v_dim self.random_field = None + self.z_kriged_field = None + self.gaussian_process = GPR() + self.conditioning_points = None + self.conditioning_values = None + self.kriging_mean = None + self.kriging_std = None + self.conditioned_random_field = None + - def generate(self, nodes: np.ndarray) -> None: + def generate(self, nodes: npt.NDArray[np.float64]) -> None: """ Generate random field - Parameters: - ------------ - nodes: list - The nodes of the random field + Args: + - nodes (ndarray): The nodes of the random field. shape (:,`self.n_dim`) + + Raises: + - ValueError: if dimensions of `nodes` do not match dimensions of the model + + Returns: + """ # check dimensions of nodes agrees with dimensions of model if nodes.shape[1] != self.n_dim: @@ -91,6 +145,134 @@ def generate(self, nodes: np.ndarray) -> None: mask = np.arange(len(scale_fluctuation)) != self.v_dim scale_fluctuation[mask] = scale_fluctuation[mask] * self.anisotropy - model = self.random_field_model(dim=self.n_dim, var=self.variance, len_scale=scale_fluctuation, angles=self.angle) - self.random_field = gs.SRF(model, mean=self.mean, seed=self.seed) + model = self.random_field_model(dim = self.n_dim, + var = self.variance, + len_scale = scale_fluctuation, + angles = self.angle) + self.random_field = gs.SRF(model, + mean = self.mean, + seed = self.seed) self.random_field(nodes.T) + + + def set_conditioning_points(self, + points: npt.NDArray[np.float64], + values: npt.NDArray[np.float64], + noise_level = 0.0001 ) -> None: + """ + Initiates the conditioning points and inverts the covariance matrix + + Args: + - points (array): The contitioning point coordinates. Dimensions `(:,self.n_dims)` + - values: array-like: The conditioning point values. + - kernel: (:class: sklearn.gaussian_process.kernels.Kernel): The (calibrated) correlation + kernel from sklearn (default = None). + - noise_level (float or array): normalised variance of the observed data, as a single float + or as an array of the same dimensions as `values` (default = 0.0001). + + Raises: + - Exception: if more than `self.max_conditioning_points` conditioning points are specified, + to prevent excessive computation and memory usage. + - ValueError: if dimensions of `points` do not match dimensions of the model + + Returns: + + """ + + self.noise_level = noise_level + self.conditioning_points = points + self.conditioning_values = values + + # check the maximum number of conditioning points to keep the computation and memory cost reasonable + if points.shape[0] > self.max_conditioning_points: + raise Exception(f'Too many conditioning points! There are {points.shape[0]} ' + + f'points, while the maximum allowed amount is {self.max_conditioning_points}.' + + f'\n Consider increasing `max_conditioning_points` or use fewer conditioning points.' ) + + # check dimensions of conditioning points agrees with dimensions of model + if points.shape[1] != self.n_dim: + raise ValueError(f'Dimensions of conditioning points: {points.shape[1]} do not match dimensions of model: {self.n_dim}') + + + # scale of fluctuation + scale_fluctuation = np.ones(self.n_dim) * self.vertical_scale_fluctuation + + # apply the anisotropy to the other dimensions + mask = np.arange(len(scale_fluctuation)) != self.v_dim + scale_fluctuation[mask] = scale_fluctuation[mask] * self.anisotropy + + # + # correct the length scales between libraries + if self.random_field_model_name == 'Gaussian': + ls_sklearn = scale_fluctuation * 2/np.pi + self.kriging_kernel = WhiteKernel(noise_level) + RBF(length_scale = ls_sklearn) + elif self.random_field_model_name == 'Exponential': + ls_sklearn = scale_fluctuation / np.sqrt(1/2) + self.kriging_kernel = WhiteKernel(noise_level) + Matern(length_scale = ls_sklearn, nu = 0.5) + elif self.random_field_model_name == 'Matern': + ls_sklearn = scale_fluctuation / np.sqrt(1/2) + self.kriging_kernel = WhiteKernel(noise_level) + Matern(length_scale = ls_sklearn) + + # fit a GP against conditioning points standardised by simulation field statistics: NOT by conditioning point statistics + self.gaussian_process = GPR(kernel = self.kriging_kernel) + self.gaussian_process.optimizer = None + self.gaussian_process.fit(points,(values - self.mean) / np.sqrt(self.variance) ) + + # reset variables + self.kriging_mean = np.array([]) + self.kriging_std = np.array([]) + self.conditioned_random_field = np.array([]) + + def generate_conditioned(self, nodes: npt.NDArray[np.float64]) -> None: + """ + Generate conditioned random field + + Args: + - nodes (array): The nodes of the random field, shape (:,`self.n_dims`) + + Raises: + - ValueError: if dimensions of nodes do not match the dimensions of the model + + Returns: + + + """ + # check dimensions of nodes agrees with dimensions of model + if nodes.shape[1] != self.n_dim: + raise ValueError(f'Dimensions of nodes: {nodes.shape[1]} do not match dimensions of model: {self.n_dim}') + # + # create kriged mean field + z_kriged_field,std_kriged_field = self.gaussian_process.predict(nodes,return_std = True) + + std_kriged_field = np.sqrt(std_kriged_field**2 - self.noise_level) + + std_kriged_field *= np.sqrt(self.variance) + self.kriging_mean = self.mean + np.sqrt(self.variance) * z_kriged_field + # + # create single random field at nodes and conditioning points + self.generate(np.vstack([nodes,self.conditioning_points])) + + # Split the generated field into the nodal coordinates and the conditioning points + # Standardize the distribution to marginal ~N(0,1) + z_cond_rf_nodes = (self.random_field.field[:nodes.shape[0]] - self.mean) / np.sqrt(self.variance) + z_rf_cond_points = (self.random_field.field[nodes.shape[0]:] - self.mean) / np.sqrt(self.variance) + + # add the noise to the conditioning points + np.random.seed(self.seed) + z_rf_cond_points += np.sqrt(self.noise_level) * np.random.normal(size = z_rf_cond_points.shape) + + # + # create kriged mean field of random field + gp_rf = GPR(kernel = self.gaussian_process.kernel_) + gp_rf.optimizer = None + gp_rf.fit(self.conditioning_points,z_rf_cond_points) + gp_rf.L_ = self.gaussian_process.L_ + # + # replace kriged mean field (the kriged parts) to create conditioned random field + z_cond_rf_nodes -= gp_rf.predict(nodes) + z_cond_rf_nodes += z_kriged_field + + # + # scale conditioned random field + self.kriging_std = np.sqrt(self.variance) * std_kriged_field + self.conditioned_random_field = self.mean + np.sqrt(self.variance) * z_cond_rf_nodes \ No newline at end of file diff --git a/random_fields/utils.py b/random_fields/utils.py index 3e804ae..3caca1a 100644 --- a/random_fields/utils.py +++ b/random_fields/utils.py @@ -7,20 +7,14 @@ def plot3D(coordinates: list, random_field: list, title: str = "Random Field", output_folder = "./", output_name: str = "random_field.png"): """ - Plot 3D random field - - Parameters - ---------- - coordinates : list - List of coordinates of the random field - random_field : list - List of random field values - title : str - Title of the plot - output_folder : str - Output folder - output_name : str - Output fine name + Plots and saves 3D random field + + Args: + - coordinates (list): List of coordinates of the random field + - random_field (list): List of random field values + - title (str): Title of the plot + - output_folder (str): Output folder + - output_name (str): Output fine name """ # create output folder if not os.path.exists(output_folder): @@ -53,22 +47,21 @@ def plot3D(coordinates: list, random_field: list, title: str = "Random Field", def plot2D(coordinates: list, random_field: list, title: str = "Random Field", - output_folder = "./", output_name: str = "random_field.png"): + output_folder = "./", output_name: str = "random_field.png", + colorbar_label: str = '',conditioning_coordinates: list = [], + conditioning_values:list = [],figsize: tuple = (6,5)): """ - Plot 2D random field - - Parameters - ---------- - coordinates : list - List of coordinates of the random field - random_field : list - List of random field values - title : str - Title of the plot - output_folder : str - Output folder - output_name : str - Output fine name + Plots and saves 2D random field + + Args: + - coordinates (list): List of coordinates of the random field + - random_field (list): List of random field values + - title (str): Title of the plot + - output_folder (str): Output folder + - output_name (str): Output fine name + - colorbar_label (str): Label for the colorbar (optional, deafult = '') + - conditioning_coordinates (list): List of coordinates if the conditioning points (optional, default = []) + - conditioning_values (list): List of the values at the conditioning points (optional, default = []) """ # create output folder @@ -76,7 +69,7 @@ def plot2D(coordinates: list, random_field: list, title: str = "Random Field", os.makedirs(output_folder) # make plot - fig, ax = plt.subplots(figsize=(6, 5)) + fig, ax = plt.subplots(figsize=figsize) ax.set_position([0.1, 0.1, 0.7, 0.8]) vmin = np.min([np.min(aux) for aux in random_field]) @@ -89,9 +82,58 @@ def plot2D(coordinates: list, random_field: list, title: str = "Random Field", ax.set_xlabel('x coordinate') ax.set_ylabel('y coordinate') + for i,coord in enumerate(conditioning_coordinates): + x, y = coord[:, 0], coord[:, 1] + ax.scatter(x, y, c=conditioning_values[i], vmin=vmin, vmax=vmax, cmap="viridis", edgecolors='red', linewidth = 0.25 ) + cax = ax.inset_axes([1.1, 0., 0.05, 1]) norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax) - fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap="viridis"), ax=ax, cax=cax) + fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap="viridis"), ax=ax, cax=cax,label = colorbar_label) + fig.suptitle(title) + plt.savefig(os.path.join(output_folder, output_name)) + + plt.close() + + +def plot1D(coordinates: list, random_field: list, title: str = "Random Field", + output_folder = "./", output_name: str = "random_field.png"): + """ + Plots and saves 1D random field + + Args: + - coordinates (list): List of coordinates of the random field + - random_field (list): List of random field values + - title (str): Title of the plot + - output_folder (str): Output folder + - output_name (str): Output fine name + + """ + + # create output folder + if not os.path.exists(output_folder): + os.makedirs(output_folder) + + # make plot + fig, ax = plt.subplots(figsize=(8, 4)) + ax.set_position([0.1, 0.1, 0.7, 0.8]) + + vmin = np.min([np.min(aux) for aux in random_field]) + vmax = np.max([np.max(aux) for aux in random_field]) + + color = 'black' + if len(coordinates) > 1: + color = 'gray' + + for i, coord in enumerate(coordinates): + x = coord[:, 0] + y = np.array(random_field[i]).ravel() + ax.plot(x, y, color = color) + + ax.set_xlabel('x coordinate') + ax.set_ylabel('random field value') + + plt.tight_layout() + fig.suptitle(title) plt.savefig(os.path.join(output_folder, output_name)) - plt.close() \ No newline at end of file + plt.close() \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 28fa2ce..87c6dbd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ numpy==1.24.2 gstools==1.5.0 matplotlib>=3.7.0 +scikit-learn==1.5.1 diff --git a/requirements_dev.txt b/requirements_dev.txt index e1ee5ed..877675e 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,5 +1,6 @@ numpy==1.24.2 gstools==1.5.0 matplotlib>=3.7.0 +scikit-learn==1.5.1 pytest==7.4.2 -pytest-cov==4.1.0 \ No newline at end of file +pytest-cov==4.1.0 diff --git a/setup.cfg b/setup.cfg index 4970b85..a0756ad 100644 --- a/setup.cfg +++ b/setup.cfg @@ -17,6 +17,7 @@ install_requires = numpy>=1.24 gstools>=1.5.0 matplotlib>=3.7.1 + scikit-learn python_requires = >=3.8 [options.extras_require] diff --git a/tests/data/conditioned_rf_3D.txt b/tests/data/conditioned_rf_3D.txt new file mode 100644 index 0000000..d3e0088 --- /dev/null +++ b/tests/data/conditioned_rf_3D.txt @@ -0,0 +1,3 @@ +1.001037129583365726e+01 1.001779265915069317e+01 1.002616224620184759e+01 1.003297125891605290e+01 1.003561413196588425e+01 1.041995334622643554e+01 1.072045835625976373e+01 1.105935873525189805e+01 1.133506851344082911e+01 1.144208343218670976e+01 1.144208343218670976e+01 1.247399161949403634e+01 1.363774617960853064e+01 1.458450968748578447e+01 1.495198965331282714e+01 1.041995334622643554e+01 1.072045835625976373e+01 1.105935873525189805e+01 1.133506851344082911e+01 1.144208343218670976e+01 1.001037129583365726e+01 1.001779265915069317e+01 1.002616224620184759e+01 1.003297125891605290e+01 1.003561413196588425e+01 1.000936434881164061e+01 1.001606517346010072e+01 1.002362215899146136e+01 1.002977008603369136e+01 1.003215636307182201e+01 1.037918016048852188e+01 1.065050919966851772e+01 1.095650580917420847e+01 1.120544698057176625e+01 1.130207184242689777e+01 1.130207184242689777e+01 1.223379227182344664e+01 1.328455813626671578e+01 1.413940056599700767e+01 1.447120197601240399e+01 1.037918016048852188e+01 1.065050919966851772e+01 1.095650580917420847e+01 1.120544698057176625e+01 1.130207184242689777e+01 1.000936434881164061e+01 1.001606517346010072e+01 1.002362215899146136e+01 1.002977008603369136e+01 1.003215636307182201e+01 1.000827212509241626e+01 1.001419138983034429e+01 1.002086695594758226e+01 1.002629781105297546e+01 1.002840576138263806e+01 1.033495396030360958e+01 1.057463616335353329e+01 1.084494243692374660e+01 1.106484801198021373e+01 1.115020290001145398e+01 1.115020290001145398e+01 1.197325083405968194e+01 1.290145917490132277e+01 1.365659588064252183e+01 1.394969717627938266e+01 1.033495396030360958e+01 1.057463616335353329e+01 1.084494243692374660e+01 1.106484801198021373e+01 1.115020290001145398e+01 1.000827212509241626e+01 1.001419138983034429e+01 1.002086695594758226e+01 1.002629781105297546e+01 1.002840576138263806e+01 1.000936434881164061e+01 1.001606517346010072e+01 1.002362215899146136e+01 1.002977008603369136e+01 1.003215636307182201e+01 1.037918016048852188e+01 1.065050919966851772e+01 1.095650580917420847e+01 1.120544698057176802e+01 1.130207184242689955e+01 1.130207184242689955e+01 1.223379227182344842e+01 1.328455813626671755e+01 1.413940056599701123e+01 1.447120197601240754e+01 1.037918016048852188e+01 1.065050919966851772e+01 1.095650580917420847e+01 1.120544698057176802e+01 1.130207184242689955e+01 1.000936434881164061e+01 1.001606517346010072e+01 1.002362215899146136e+01 1.002977008603369136e+01 1.003215636307182201e+01 1.001037129583365726e+01 1.001779265915069317e+01 1.002616224620184759e+01 1.003297125891605290e+01 1.003561413196588425e+01 1.041995334622643554e+01 1.072045835625976373e+01 1.105935873525189805e+01 1.133506851344083088e+01 1.144208343218671153e+01 1.144208343218671153e+01 1.247399161949403634e+01 1.363774617960853242e+01 1.458450968748578447e+01 1.495198965331282892e+01 1.041995334622643554e+01 1.072045835625976373e+01 1.105935873525189805e+01 1.133506851344083088e+01 1.144208343218671153e+01 1.001037129583365726e+01 1.001779265915069317e+01 1.002616224620184759e+01 1.003297125891605290e+01 1.003561413196588425e+01 +1.999995657040991581e+00 1.999987217896785463e+00 1.999972364185357732e+00 1.999956107022822005e+00 1.999948788263586463e+00 1.992866601328432852e+00 1.978931651574280304e+00 1.954163581335734312e+00 1.926690606096852720e+00 1.914194004823318807e+00 1.914194004823318807e+00 1.735368065190297182e+00 1.364846429243236736e+00 7.781931346967962959e-01 1.990064074795734050e-01 1.992866601328432852e+00 1.978931651574280304e+00 1.954163581335734312e+00 1.926690606096852720e+00 1.914194004823318807e+00 1.999995657040991581e+00 1.999987217896785463e+00 1.999972364185357732e+00 1.999956107022822005e+00 1.999948788263586463e+00 1.999997106204023156e+00 1.999991483048136676e+00 1.999981585768210390e+00 1.999970753373113919e+00 1.999965876781916618e+00 1.995249714589760526e+00 1.985986601500497439e+00 1.969576902793994444e+00 1.951458722886711294e+00 1.943247064898389942e+00 1.943247064898389942e+00 1.827934628395943450e+00 1.604975522519639597e+00 1.318421698327810976e+00 1.166666943669037559e+00 1.995249714589760526e+00 1.985986601500497439e+00 1.969576902793994444e+00 1.951458722886711294e+00 1.943247064898389942e+00 1.999997106204023156e+00 1.999991483048136676e+00 1.999981585768210390e+00 1.999970753373113919e+00 1.999965876781916618e+00 1.999998481827183427e+00 1.999995531752665645e+00 1.999990339355364410e+00 1.999984656383788861e+00 1.999982097995069097e+00 1.997509264526410533e+00 1.992660433959572153e+00 1.984097289174038492e+00 1.974682568141426842e+00 1.970429526712960788e+00 1.970429526712960788e+00 1.911661467756611721e+00 1.803579276266119802e+00 1.677322190641976185e+00 1.617268996638918876e+00 1.997509264526410533e+00 1.992660433959572153e+00 1.984097289174038492e+00 1.974682568141426842e+00 1.970429526712960788e+00 1.999998481827183427e+00 1.999995531752665645e+00 1.999990339355364410e+00 1.999984656383788861e+00 1.999982097995069097e+00 1.999997106204023156e+00 1.999991483048136676e+00 1.999981585768210390e+00 1.999970753373113919e+00 1.999965876781916618e+00 1.995249714589760526e+00 1.985986601500497439e+00 1.969576902793994444e+00 1.951458722886711294e+00 1.943247064898389942e+00 1.943247064898389942e+00 1.827934628395943228e+00 1.604975522519638931e+00 1.318421698327809866e+00 1.166666943669036449e+00 1.995249714589760526e+00 1.985986601500497439e+00 1.969576902793994444e+00 1.951458722886711294e+00 1.943247064898389942e+00 1.999997106204023156e+00 1.999991483048136676e+00 1.999981585768210390e+00 1.999970753373113919e+00 1.999965876781916618e+00 1.999995657040991581e+00 1.999987217896785463e+00 1.999972364185357732e+00 1.999956107022822005e+00 1.999948788263586463e+00 1.992866601328432852e+00 1.978931651574280304e+00 1.954163581335734312e+00 1.926690606096852720e+00 1.914194004823318807e+00 1.914194004823318807e+00 1.735368065190297182e+00 1.364846429243236292e+00 7.781931346967956298e-01 1.990064074795688809e-01 1.992866601328432852e+00 1.978931651574280304e+00 1.954163581335734312e+00 1.926690606096852720e+00 1.914194004823318807e+00 1.999995657040991581e+00 1.999987217896785463e+00 1.999972364185357732e+00 1.999956107022822005e+00 1.999948788263586463e+00 +1.012121647287858828e+01 1.071642554112682433e+01 1.108706117485346532e+01 1.120483867349165230e+01 1.114718970010125254e+01 1.055255019288573948e+01 1.130096185115665541e+01 1.208263360116385599e+01 1.269063991410385661e+01 1.296262105045423851e+01 1.249046303079041920e+01 1.291513149953474482e+01 1.364166258405261445e+01 1.438327350547505468e+01 1.472565541955533774e+01 1.176553517165141116e+01 1.153782641372943552e+01 1.141882964402100775e+01 1.146858862202424589e+01 1.160306592896929345e+01 1.001821773354341083e+01 1.020704159418873935e+01 1.034185632331666760e+01 1.047834728536685134e+01 1.068560737496155078e+01 1.073332562848649019e+01 1.091155365997554050e+01 1.109075129059374376e+01 1.120468639540295008e+01 1.125967414385540621e+01 1.054848864482636017e+01 1.046444561337704116e+01 1.028533875920313889e+01 1.010886811447745792e+01 1.004008457220642470e+01 1.258417437649817927e+01 1.304503419660922425e+01 1.362352677425436553e+01 1.409601635471898362e+01 1.418229174957941652e+01 1.211633594204680264e+01 1.230357056106962332e+01 1.254281838105818281e+01 1.268759867617738024e+01 1.263618346678125270e+01 1.058019197852246762e+01 1.038784009647453033e+01 1.029028550993657731e+01 1.026937331533787301e+01 1.030195309399787540e+01 8.883457927525110875e+00 8.656185766708823337e+00 8.425990024817773261e+00 8.241401247266487928e+00 8.159267121028348058e+00 1.009317641877840721e+01 1.014821670668675857e+01 1.015124644792975417e+01 1.013697310877260094e+01 1.010591861841817440e+01 1.011313890100671742e+01 1.118771109316880796e+01 1.229989740825442013e+01 1.321485100044520422e+01 1.364153170972529949e+01 7.331775718706962763e+00 8.283737075605337452e+00 9.134573268233307175e+00 9.814223526220782645e+00 1.027626074267490708e+01 9.353470713313768670e+00 9.083221003139598793e+00 8.821853352654798641e+00 8.742218087186623521e+00 8.934407406326720036e+00 1.122724383391061131e+01 1.151133619880937609e+01 1.155624267179873677e+01 1.138088448336473135e+01 1.107052906417530735e+01 8.480559650336589428e+00 8.972690491930276124e+00 9.601907045596810519e+00 1.025873694367369460e+01 1.079059351667883604e+01 9.037161152833929734e+00 1.012492132175658988e+01 1.146858042524979382e+01 1.271393564030815959e+01 1.347147530704844343e+01 1.011312904807348723e+01 1.053637517682219205e+01 1.106544555215963399e+01 1.156124807893860407e+01 1.186194429368405601e+01 9.601100668704498631e+00 9.431434050046169304e+00 9.406497778774722818e+00 9.623231490733354931e+00 1.008690053298261091e+01 9.075495998473334325e+00 9.409838886026964389e+00 9.780151585061991781e+00 1.019000517041948584e+01 1.067142913870968357e+01 8.589795774759629410e+00 9.184796607457128559e+00 9.867636524483291538e+00 1.047226059447149638e+01 1.085065198604466197e+01 9.707465276404960264e+00 1.131875262321031173e+01 1.302971604180225818e+01 1.436569418399955644e+01 1.497682066499356246e+01 8.788740439515176917e+00 9.441502269066109676e+00 1.018624732643457342e+01 1.080419838736380278e+01 1.117003755917914631e+01 9.242131500766898711e+00 9.322385502862839957e+00 9.519349270709991373e+00 9.693853564647792354e+00 9.779321257325094052e+00 diff --git a/tests/test_random_field.py b/tests/test_random_field.py index d766764..404a8c2 100644 --- a/tests/test_random_field.py +++ b/tests/test_random_field.py @@ -86,14 +86,62 @@ def test_distribution_RF_unstruc(): x = np.array(nodes_fields)[:, 0] y = np.array(nodes_fields)[:, 1] - rf_emsamble = [] + rf_ensamble = [] param = [] for i in range(nb_runs): rf = RandomFields(ModelName.Gaussian, 2, 10, 2, 0.05, [1], [1], seed=i) rf.generate(np.array([x.ravel(), y.ravel()]).T) - rf_emsamble.append(rf.random_field.field) + rf_ensamble.append(rf.random_field.field) mu_un, std_un = norm.fit(rf.random_field.field) param.append([mu_un, std_un]) np.testing.assert_array_almost_equal(np.mean(np.array(param)[:, 0]), 10, decimal=2) np.testing.assert_array_almost_equal(np.mean(np.array(param)[:, 1])**2, 2, decimal=2) + + +def test_conditioned_RF_3D_mean_variance(): + """test the conditioned random field for correct mean and variance""" + + # mesh coordinates + x = np.linspace(0, 100, 5) + y = np.linspace(0, 50, 5) + z = np.linspace(0, 25, 5) + x, y, z = [i.ravel() for i in np.meshgrid(x, y, z)] + + # random field properties + nb_dimensions = 3 + mean = 10 + variance = 2 + vertical_scale_fluctuation = 10 + anisotropy = [2.5, 2.5] + angle = [0, 0] + model_rf = ModelName.Gaussian + + # generate and plot random field + rf = RandomFields(model_rf, nb_dimensions, mean, variance, vertical_scale_fluctuation, + anisotropy, angle, seed=14) + rf.generate(np.array([x, y, z]).T) + + # declare conditioning points + xc = np.array([50.]*4) + yc = np.linspace(0,50,4) + zc = np.array([25]*4) + vc = np.array([15]*4) + + rf.set_conditioning_points(np.array([xc,yc,zc]).T,vc,noise_level = 0.01) + + # generate conditioned random field model + rf.generate_conditioned(np.array([x, y,z]).T) + + # load reference solution + data_ref = np.loadtxt('./tests/data/conditioned_rf_3D.txt') + kriging_mean_ref = data_ref[0] + kriging_std_ref = data_ref[1] + kriging_field_ref = data_ref[2] + + # test difference in mean, std and single realisation + np.testing.assert_array_almost_equal(kriging_mean_ref, rf.kriging_mean, decimal=4) + np.testing.assert_array_almost_equal(kriging_std_ref, rf.kriging_std, decimal=4) + np.testing.assert_array_almost_equal(kriging_field_ref, rf.conditioned_random_field, decimal=4) + +