diff --git a/.gitignore b/.gitignore index 8b9dee82..bc1db546 100644 --- a/.gitignore +++ b/.gitignore @@ -176,4 +176,7 @@ poetry.toml # LSP config files pyrightconfig.json +# For the documentation +_autosummary + # End of https://www.toptal.com/developers/gitignore/api/python diff --git a/README.md b/README.md index e066cec8..554b13c2 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,15 @@ Iterative Ensemble Smoother =========================== +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://github.com/equinor/iterative_ensemble_smoother/blob/main/COPYING) +[![Stars](https://img.shields.io/github/stars/equinor/iterative_ensemble_smoother.svg?style=social&label=Star&maxAge=2592000)](https://github.com/equinor/iterative_ensemble_smoother/stargazers) +[![Python](https://img.shields.io/pypi/pyversions/iterative_ensemble_smoother.svg)](https://pypi.org/pypi/iterative_ensemble_smoother) +[![PyPI](https://img.shields.io/pypi/v/iterative_ensemble_smoother.svg)](https://pypi.org/pypi/iterative_ensemble_smoother) +[![Downloads](https://static.pepy.tech/badge/iterative_ensemble_smoother)](https://pepy.tech/project/iterative_ensemble_smoother) +[![Build Status](https://github.com/equinor/iterative_ensemble_smoother/actions/workflows/upload_to_pypi.yml/badge.svg)](https://github.com/equinor/iterative_ensemble_smoother/actions/workflows/main.yml) [![Precommit: enabled](https://img.shields.io/badge/pre--commit-enabled-brightgreen?logo=pre-commit)](https://github.com/pre-commit/pre-commit) [![Ruff](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/astral-sh/ruff/main/assets/badge/v2.json)](https://github.com/astral-sh/ruff) +[![Mypy](https://www.mypy-lang.org/static/mypy_badge.svg)](https://mypy-lang.org/) [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black) [![docs](https://readthedocs.org/projects/iterative_ensemble_smoother/badge/?version=latest&style=plastic)](https://iterative-ensemble-smoother.readthedocs.io/) @@ -21,31 +28,9 @@ pip install iterative_ensemble_smoother ## Usage -The following illustrates usage but does not actually create the inputs and can not be run. +**iterative_ensemble_smoother** mainly implements `SIES` and `ESMDA` classes. Check out +the examples section to see how to use them. -The example below illustrates the usage of the package. -Note that it does not create the necessary inputs and cannot be run directly. -For complete information and runnable examples, please refer to [the documentation](http://iterative_ensemble_smoother.rtfd.io). - -```python -import iterative_ensemble_smoother as ies - -# `prior` is a matrix of size (num_params, ensemble_size). -# In geostatistical applications, it typically consists of Gaussian random fields -# that model rock properties like porosities and permeabilities. - -# `response_ensemble` is a matrix of size (num_obs, ensemble_size) and is -# generated by running a dynamical model, such as a flow simulator. -# In geostatistical applications, this is typically a porous flow simulator like Eclipse or OPM flow. - -# `obs_errors` and `obs_values` are 1D array of size `num_obs` that hold observations from real-life -# sensors and their uncertainty specifications. -# In geostatistical applications, the observations are typically pressures, temperatures, production rates etc. - -smoother = ies.ES() -smoother.fit(response_ensemble, obs_errors, obs_values) -posterior = smoother.update(prior) -``` ## Building from source @@ -57,7 +42,7 @@ cd iterative_ensemble_smoother pip install . ``` -### Building the documentation +## Building the documentation ```bash apt install pandoc # Pandoc is required to build the documentation. diff --git a/docs/source/_templates/autosummary/base.rst b/docs/source/_templates/autosummary/base.rst new file mode 100644 index 00000000..8d96add1 --- /dev/null +++ b/docs/source/_templates/autosummary/base.rst @@ -0,0 +1,11 @@ +.. + Template for the html base rendering + + Modified from + https://github.com/sphinx-doc/sphinx/tree/master/sphinx/ext/autosummary/templates/autosummary/base.rst + +{{ fullname | escape | underline}} + +.. currentmodule:: {{ module }} + +.. auto{{ objtype }}:: {{ objname }} diff --git a/docs/source/_templates/autosummary/class.rst b/docs/source/_templates/autosummary/class.rst new file mode 100644 index 00000000..318721b0 --- /dev/null +++ b/docs/source/_templates/autosummary/class.rst @@ -0,0 +1,21 @@ +.. + Template for the html class rendering + + Modified from + https://github.com/sphinx-doc/sphinx/tree/master/sphinx/ext/autosummary/templates/autosummary/class.rst + +{{ fullname | escape | underline}} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :members: + :private-members: + :show-inheritance: + :inherited-members: + + {% block methods %} + .. automethod:: __init__ + {% endblock %} + + .. rubric:: {{ _('Methods definition') }} diff --git a/docs/source/_templates/autosummary/module.rst b/docs/source/_templates/autosummary/module.rst new file mode 100644 index 00000000..85ba96c8 --- /dev/null +++ b/docs/source/_templates/autosummary/module.rst @@ -0,0 +1,67 @@ +.. + Template for the html module rendering + + Modified from + https://github.com/sphinx-doc/sphinx/tree/master/sphinx/ext/autosummary/templates/autosummary/module.rst + +{{ fullname | escape | underline}} + +.. automodule:: {{ fullname }} + + {% block attributes %} + {% if attributes %} + .. rubric:: {{ _('Module Attributes') }} + + .. autosummary:: + {% for item in attributes %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block functions %} + {% if functions %} + .. rubric:: {{ _('Functions') }} + + .. autosummary:: + {% for item in functions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block classes %} + {% if classes %} + .. rubric:: {{ _('Classes') }} + + .. autosummary:: + {% for item in classes %} + :template: class.rst + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + + {% block exceptions %} + {% if exceptions %} + .. rubric:: {{ _('Exceptions') }} + + .. autosummary:: + {% for item in exceptions %} + {{ item }} + {%- endfor %} + {% endif %} + {% endblock %} + +{% block modules %} +{% if modules %} +.. rubric:: Modules + +.. autosummary:: + :toctree: + :recursive: +{% for item in modules %} + {{ item }} +{%- endfor %} +{% endif %} +{% endblock %} diff --git a/docs/source/api.rst b/docs/source/api.rst index a00951cb..c2896973 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -1,10 +1,11 @@ -API reference +.. _api_ref: + +============= +API Reference ============= -.. currentmodule:: iterative_ensemble_smoother +.. automodule:: iterative_ensemble_smoother -.. autoclass:: iterative_ensemble_smoother.SIES - :members: +.. raw:: latex -.. autoclass:: iterative_ensemble_smoother.ESMDA - :members: + \clearpage diff --git a/docs/source/bibliography.rst b/docs/source/bibliography.rst index 0d21440d..eb35d1fa 100644 --- a/docs/source/bibliography.rst +++ b/docs/source/bibliography.rst @@ -2,3 +2,5 @@ Bibliography ============ .. bibliography:: + :all: + :labelprefix: Bib- diff --git a/docs/source/conf.py b/docs/source/conf.py index fd4b3d99..53b5801d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,26 +1,168 @@ +""" +iterative_ensemble_smoother documentation build configuration file. +""" +import datetime +import os +import re +import sys from subprocess import check_output +import iterative_ensemble_smoother as ies + +package_path = os.path.abspath("..") +sys.path.insert(0, package_path) + project = "iterative_ensemble_smoother" -copyright = "2022, Equinor" author = "Equinor" -release = "0.1.1" - +copyright = f"2022-{datetime.datetime.today().year}, {author}" +# The default replacements for |version| and |release|, also used in various +# other places throughout the built documents. +version = re.sub(r"\.dev.*$", r".dev", ies.__version__) +release = version +# convert the python file to a notebook check_output(["jupytext", "Polynomial.py", "-o", "Polynomial.ipynb"]) +# Do the same for this file check_output(["jupytext", "Oscillator.py", "-o", "Oscillator.ipynb"]) - +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom ones. extensions = [ - "sphinx.ext.autodoc", - "sphinx.ext.doctest", - "nbsphinx", + "sphinx.ext.napoleon", # autodoc understands numpy docstrings + "sphinx.ext.autodoc", # Core library for html generation from docstrings + "sphinx.ext.autosummary", # Create neat summary tables + "sphinx.ext.doctest", # Test snippets in the documentation + "sphinx.ext.viewcode", # Add links to highlighted source code + "sphinx.ext.intersphinx", # Link to other projects’ documentation + "sphinx.ext.autosectionlabel", # Allow reference sections using its title + # allows BibTeX citations to be inserted into documentation generated by Sphinx "sphinxcontrib.bibtex", + "jupyter_sphinx", "numpydoc", + "m2r2", # compatibility with markdown + "myst_nb", # to include jupyter notebooks ] -bibtex_bibfiles = ["refs.bib"] + +# ----------------------------------------------------------------------------- +# Autosummary +# ----------------------------------------------------------------------------- + +# autosummaries from source-files +autosummary_generate = True +# dont show __init__ docstring +autoclass_content = "class" +# sort class members +autodoc_member_order = "groupwise" +# autodoc_member_order = 'bysource' + +autosummary_generate = True # Turn on sphinx.ext.autosummary +# This is because numpydoc screew up with autosummary +# numpydoc_show_class_members=False + +# Napoleon settings +napoleon_google_docstring = True +napoleon_numpy_docstring = True +napoleon_include_init_with_doc = False +napoleon_include_private_with_doc = True +napoleon_include_special_with_doc = False +# napoleon_use_admonition_for_examples = False +# napoleon_use_admonition_for_notes = False +# napoleon_use_admonition_for_references = False +napoleon_use_ivar = True +# napoleon_use_param = True +# napoleon_use_rtype = True + +# Add any paths that contain templates here, relative to this directory. +templates_path = ["_templates"] + +# Add any external modules you want to refer to in the docs here. +intersphinx_mapping = { + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), +} + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# source_suffix = ['.rst', '.md'] +source_suffix = [".rst", ".md", ".ipynb"] + +# The encoding of source files. +# source_encoding = 'utf-8-sig' + +# The master toctree document. +master_doc = "index" + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. language = "en" -html_theme = "sphinx_rtd_theme" +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +exclude_patterns = ["build", "_templates/*.rst'", "Thumbs.db", ".DS_Store"] + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = "sphinx" + +# ----------------------------------------------------------------------------- +# HTML output +# ----------------------------------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# html_theme = "sphinx_rtd_theme" +html_theme = "pydata_sphinx_theme" + +# TODO: when ies gets a nice logo +# html_logo = '_static/logo.svg' +# html_favicon = '_static/favicon.ico' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +# +html_theme_options = { + # "google_analytics_id": "UA-140243896-1", + "show_prev_next": False, + "github_url": "https://github.com/equinor/iterative_ensemble_smoother", + "icon_links": [ + { + "name": "Support", + "url": "https://github.com/equinor/iterative_ensemble_smoother/issues", + "icon": "fa fa-comment fa-fw", + }, + # { + # "name": "The Paper", + # "url": "https://doi.org/10.21105/joss.01450", + # "icon": "fa fa-file-text fa-fw", + # }, + ], +} + +# Title displayed on the html page +html_title = "Iterative Ensemble Smoother" + +# ----------------------------------------------------------------------------- +# Bibliography +# ----------------------------------------------------------------------------- + +bibtex_bibfiles = ["./refs.bib"] +bibtex_default_style = "unsrt" +bibtex_reference_style = "author_year" +suppress_warnings = ["bibtex.duplicate_citation", "autosectionlabel.*"] numpydoc_class_members_toctree = False + +# ----------------------------------------------------------------------------- +# Myst configuration +# ----------------------------------------------------------------------------- + +# see https://myst-parser.readthedocs.io/en/latest/syntax/optional.html#syntax-mathjax +# for math parsing in jupyter notebooks +myst_enable_extensions = [ + "amsmath", + "dollarmath", +] diff --git a/docs/source/examples.rst b/docs/source/examples.rst index 55edb7b0..d6506b3e 100644 --- a/docs/source/examples.rst +++ b/docs/source/examples.rst @@ -7,3 +7,5 @@ Example Usage Polynomial Oscillator + +* :ref:`genindex` diff --git a/docs/source/glossary.rst b/docs/source/glossary.rst index 4e85de9d..8551ee35 100644 --- a/docs/source/glossary.rst +++ b/docs/source/glossary.rst @@ -17,8 +17,8 @@ Glossary forward model A forward model maps model parameters to predicted measurements - (:term:`prediction`). See, for instance, :cite:`evensen2018analysis` and - :cite:`evensen2019efficient`. + (:term:`prediction`). See, for instance, :cite:t:`evensen2018analysis` and + :cite:t:`evensen2019efficient`. error function The error function, or erf_, is a function that for a diff --git a/docs/source/index.rst b/docs/source/index.rst index 37976b64..fc61b3f1 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,6 +1,6 @@ -=========================== -Iterative Ensemble Smoother -=========================== +**Date**: |today| **Version**: |version| + +.. mdinclude:: ../../README.md .. toctree:: :hidden: @@ -12,26 +12,6 @@ Iterative Ensemble Smoother bibliography -A library for the iterative ensemble smoothers. - -.. currentmodule:: iterative_ensemble_smoother - -Currently, two main algorithms are implemented: - -* :class:`SIES` - Subspace Iterative Ensemble Smoother - based on the method developed in :cite:`evensen2018analysis`. -* :class:`ESMDA` - Ensemble Smoother with Multiple Data Assimilation - based on the method developed in :cite:`EMERICK2013`. - - -Installation -============ - -Install from `PyPI `_: - -.. code-block:: console - - pip install iterative_ensemble_smoother - - Indices and tables ================== diff --git a/docs/source/refs.bib b/docs/source/refs.bib index aa587b04..74daf105 100644 --- a/docs/source/refs.bib +++ b/docs/source/refs.bib @@ -43,3 +43,18 @@ @article{EMERICK2013 url = {https://www.sciencedirect.com/science/article/pii/S0098300412000994}, keywords = {History matching, Ensemble smoother, Ensemble Kalman filter, Multiple data assimilation}, } + +@article{emerickHistoryMatchingTimelapse2012, + title = {History Matching Time-Lapse Seismic Data Using the Ensemble {{Kalman}} Filter with Multiple Data Assimilations}, + author = {Emerick, Alexandre A. and Reynolds, Albert C.}, + year = {2012}, + month = jun, + journal = {Computational Geosciences}, + volume = {16}, + number = {3}, + pages = {639--659}, + issn = {1573-1499}, + doi = {10.1007/s10596-012-9275-5}, + langid = {english}, + keywords = {Ensemble Kalman filter,Multiple data assimilations,Time-lapse seismic}, +} diff --git a/pyproject.toml b/pyproject.toml index 95e99431..8fea3bfe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,19 +39,22 @@ text = "GPL-3.0" [project.optional-dependencies] doc = [ "sphinx", - "sphinx-rtd-theme", - "nbsphinx", + "pydata_sphinx_theme", + "jupyter_sphinx", "sphinxcontrib.bibtex", "matplotlib", + "pygments", "jupytext", "pandas", - "scipy", "p_tqdm", "ipykernel", "ipywidgets", "numpydoc", + "scipy", + "jinja2", + "m2r2", + "myst_nb" ] - dev = [ "pytest", "pytest-snapshot", diff --git a/src/iterative_ensemble_smoother/__init__.py b/src/iterative_ensemble_smoother/__init__.py index 6c9f1215..6464a67c 100644 --- a/src/iterative_ensemble_smoother/__init__.py +++ b/src/iterative_ensemble_smoother/__init__.py @@ -1,7 +1,31 @@ -""" Implementation of the iterative ensemble smoother history matching algorithms -from Evensen et al. "Efficient Implementation of an Iterative Ensemble Smoother -for Data Assimilation and Reservoir History Matching" -https://www.frontiersin.org/articles/10.3389/fams.2019.00047/full +""" +Purpose +======= + +**iterative_ensemble_smoother** is an open-source, pure python, +and object-oriented library that provides +a user friendly implementation of history matching algorithms +from :cite:t:`evensen2019efficient`. + +The following functionalities are directly provided on module-level. + +Classes +======= + +.. autosummary:: + :toctree: _autosummary + + SIES + ESMDA + +Functions +========= + +.. autosummary:: + :toctree: _autosummary + + steplength_exponential + """ try: from ._version import version as __version__ diff --git a/src/iterative_ensemble_smoother/esmda.py b/src/iterative_ensemble_smoother/esmda.py index 8255d5f3..2707bfea 100644 --- a/src/iterative_ensemble_smoother/esmda.py +++ b/src/iterative_ensemble_smoother/esmda.py @@ -21,6 +21,7 @@ https://helper.ipam.ucla.edu/publications/oilws3/oilws3_14147.pdf """ + import numbers from typing import Optional, Union @@ -36,9 +37,10 @@ class ESMDA: - """Initialize Ensemble Smoother with Multiple Data Assimilation (ES-MDA). + """ + Implement an Ensemble Smoother with Multiple Data Assimilation (ES-MDA). - The implementation follows the 2013 paper by Emerick et al. + The implementation follows :cite:t:`EMERICK2013`. Parameters ---------- @@ -89,6 +91,7 @@ def __init__( seed: Union[np.random._generator.Generator, int, None] = None, inversion: str = "exact", ) -> None: + """Initialize the instance.""" # Validate inputs if not (isinstance(covariance, np.ndarray) and covariance.ndim in (1, 2)): raise TypeError( diff --git a/src/iterative_ensemble_smoother/esmda_inversion.py b/src/iterative_ensemble_smoother/esmda_inversion.py index 5d1586d9..549a9c85 100644 --- a/src/iterative_ensemble_smoother/esmda_inversion.py +++ b/src/iterative_ensemble_smoother/esmda_inversion.py @@ -77,7 +77,7 @@ def empirical_cross_covariance( def normalize_alpha(alpha: npt.NDArray[np.double]) -> npt.NDArray[np.double]: """Assure that sum_i (1/alpha_i) = 1. - This is Eqn (22) in the 2013 Emerick paper. + This is Eqn (22) in :cite:t:`EMERICK2013`. Examples -------- @@ -135,10 +135,10 @@ def singular_values_to_keep( # # C_MD @ inv(C_DD + alpha * C_D) @ (D - Y) # -# where C_MD = empirical_cross_covariance(X, Y) = -# center(X) @ center(Y).T / (X.shape[1] - 1) -# C_DD = empirical_cross_covariance(Y, Y) = -# center(Y) @ center(Y).T / (Y.shape[1] - 1) +# where C_MD = empirical_cross_covariance(X, Y) = center(X) @ center(Y).T +# / (X.shape[1] - 1) +# C_DD = empirical_cross_covariance(Y, Y) = center(Y) @ center(Y).T +# / (Y.shape[1] - 1) # # The methods can be classified as # - exact : with truncation=1.0, these methods compute the exact solution @@ -264,7 +264,9 @@ def inversion_exact_rescaled( ) -> npt.NDArray[np.double]: """Compute a rescaled inversion. - See Appendix A.1 in Emerick et al (2012) for details regarding this approach.""" + See Appendix A.1 in :cite:t:`emerickHistoryMatchingTimelapse2012` + for details regarding this approach. + """ C_DD = empirical_cross_covariance(Y, Y) if C_D.ndim == 2: @@ -338,9 +340,6 @@ def inversion_exact_subspace_woodbury( (A + U @ U.T)^-1 = A^-1 - A^-1 @ U @ (1 + U.T @ A^-1 @ U )^-1 @ U.T @ A^-1 to compute inv(C_DD + alpha * C_D). - - - """ # Woodbury: @@ -397,7 +396,7 @@ def inversion_subspace( X: npt.NDArray[np.double], truncation: float = 1.0, ) -> npt.NDArray[np.double]: - """See Appendix A.2 in Emerick et al (2012) + """See Appendix A.2 in :cite:t:`emerickHistoryMatchingTimelapse2012`. This is an approximate solution. The approximation is that when U, w, V.T = svd(D_delta) @@ -484,10 +483,10 @@ def inversion_rescaled_subspace( X: npt.NDArray[np.double], truncation: float = 1.0, ) -> npt.NDArray[np.double]: - """See Appendix A.2 in Emerick et al (2012) + """ + See Appendix A.2 in :cite:t:`emerickHistoryMatchingTimelapse2012`. Subspace inversion with rescaling. - """ # TODO: I don't understand why this approach is not approximate, when # `inversion_subspace` is approximate diff --git a/src/iterative_ensemble_smoother/sies.py b/src/iterative_ensemble_smoother/sies.py index 21bc8f61..2ba5a1e1 100644 --- a/src/iterative_ensemble_smoother/sies.py +++ b/src/iterative_ensemble_smoother/sies.py @@ -17,50 +17,13 @@ class SIES: - """ - Initialize a Subspace Iterative Ensemble Smoother (SIES) instance. - - This is an implementation of the algorithm described in the paper: - Efficient Implementation of an Iterative Ensemble Smoother for - Data Assimilation and Reservoir History Matching - written by Evensen et al (2019), - URL: https://www.frontiersin.org/articles/10.3389/fams.2019.00047/full - - Parameters - ---------- - parameters : npt.NDArray[np.double] - A 2D array of shape (num_parameters, ensemble_size). Each row corresponds - to a parameter in the model, and each column corresponds to an ensemble - member (realization). This is X in Evensen (2019). - covariance : npt.NDArray[np.double] - Either a 1D array of diagonal covariances, or a 2D covariance matrix. - The shape is either (num_observations,) or (num_observations, num_observations). - This is C_dd in Evensen (2019), and represents observation or measurement - errors. We observe d from the real world, y from the model g(x), and - assume that d = y + e, where the error e is multivariate normal with - covariance given by `covariance`. - observations : npt.NDArray[np.double] - A 1D array of observations, with shape (num_observations,). - This is d in Evensen (2019). - inversion : str - The type of inversion used in the algorithm. Every inversion method - scales the variables. The default is `subspace.` - The options are: - - - `direct`: Solve Eqn (42) directly, which involves inverting a - matrix of shape (num_parameters, num_parameters). - - `subspace_exact` : Solve Eqn (42) using Eqn (50), i.e., the Woodbury - lemma to invert a matrix of size (ensemble_size, ensemble_size). - - `subspace_projected` : Solve Eqn (42) using Section 3.3, i.e., - by projecting the covariance onto S. This approach utilizes the - truncation factor `truncation`. - - truncation : float - How much of the total energy (singular values squared) to keep in the - SVD when `inversion` equals `subspace_projected`. Choosing 1.0 - retains all information, while 0.0 removes all information. - seed : Union[None, int, np.random._generator.Generator], optional - Integer used to seed the random number generator. The default is None. + r""" + Implement a Sequential Iterative Ensemble Smoother (SIES) for Data Assimilation. + + This is an implementation of the algorithm described in the paper: "Efficient + Implementation of an Iterative Ensemble Smoother for Data Assimilation + and Reservoir History Matching", written by :cite:t:`evensen2019efficient`. + """ inversion_funcs = { @@ -78,7 +41,50 @@ def __init__( inversion: str = "subspace_exact", truncation: float = 1.0, seed: Union[None, int, np.random._generator.Generator] = None, - ): + ) -> None: + """ + Initialize the instance. + + Parameters + ---------- + parameters : npt.NDArray[np.double] + A 2D array of shape (num_parameters, ensemble_size). Each row corresponds + to a parameter in the model, and each column corresponds to an ensemble + member (realization). This is X in Evensen (2019). + covariance : npt.NDArray[np.double] + Either a 1D array of diagonal covariances, or a 2D covariance matrix. + The shape is either (num_observations,) or (num_observations, + num_observations). + This is C_dd in Evensen (2019), and represents observation or measurement + errors. We observe d from the real world, y from the model g(x), and + assume that d = y + e, where the error e is multivariate normal with + covariance given by `covariance`. + observations : npt.NDArray[np.double] + A 1D array of observations, with shape (num_observations,). + This is d in Evensen (2019). + inversion : str + The type of inversion used in the algorithm. Every inversion method + scales the variables. The default is `subspace.` + The options are: + + * `direct`: + Solve Eqn (42) directly, which involves inverting a + matrix of shape (num_parameters, num_parameters). + * `subspace_exact` : + Solve Eqn (42) using Eqn (50), i.e., the Woodbury + lemma to invert a matrix of size (ensemble_size, ensemble_size). + * `subspace_projected` : + Solve Eqn (42) using Section 3.3, i.e., by projecting the covariance + onto S. This approach utilizes the truncation factor `truncation`. + + truncation : float + How much of the total energy (singular values squared) to keep in the + SVD when `inversion` equals `subspace_projected`. Choosing 1.0 + retains all information, while 0.0 removes all information. + The default is 1.0. + seed : Union[None, int, np.random._generator.Generator], optional + Integer used to seed the random number generator. The default is None. + """ _validate_inputs( parameters=parameters, covariance=covariance, observations=observations ) @@ -203,7 +209,6 @@ def sies_iteration( This method implements lines 4-9 in Algorithm 1. It returns an updated X and updates the internal state W. - Parameters ---------- responses : npt.NDArray[np.double] @@ -240,7 +245,7 @@ def sies_iteration( def propose_W( self, responses: npt.NDArray[np.double], step_length: float = 0.5 ) -> npt.NDArray[np.double]: - """Returns a proposal for W_i, without updating the internal W. + """Return a proposal for W_i, without updating the internal W. This is an implementation of lines 4-8 in Algorithm 1. @@ -322,7 +327,7 @@ def propose_W_masked( ensemble_mask: npt.NDArray[np.bool_], step_length: float = 0.5, ) -> npt.NDArray[np.double]: - """Returns a proposal for W_i, without updating the internal W. + """Return a proposal for W_i, without updating the internal W. This is an implementation of lines 4-8 in Algorithm 1. diff --git a/src/iterative_ensemble_smoother/utils.py b/src/iterative_ensemble_smoother/utils.py index b2e55040..ab2190a3 100644 --- a/src/iterative_ensemble_smoother/utils.py +++ b/src/iterative_ensemble_smoother/utils.py @@ -14,13 +14,12 @@ def steplength_exponential( max_steplength: float = 0.6, halflife: float = 1.5, ) -> float: - """ - This is an implementation of Eq. (49), which calculates a suitable step length for - the update step, from the book: + r""" + Compute a suitable step length for the update step. - Geir Evensen, Formulating the history matching problem with - consistent error statistics, - Computational Geosciences (2021) 25:945 –970 + This is an implementation of Eq. (49), which calculates a suitable step length for + the update step, from the book: \"Formulating the history matching problem with + consistent error statistics", written by :cite:t:`evensen2021formulating`. Examples --------