From 6930cd69e0eabca58ee6486e968cf38deb305e73 Mon Sep 17 00:00:00 2001 From: antoinecollet5 Date: Tue, 3 Oct 2023 23:31:02 +0200 Subject: [PATCH] DOCS: improved documentation using citet and citep for references, makefiles for unix and windows to avoid having to use tox, use pydata theme rather than read the docs and improve docstrings for classes SIES and ES. Add badges for project visibility. --- .gitignore | 3 + README.md | 41 ++- docs/Makefile | 20 ++ docs/make.bat | 36 +++ docs/source/_templates/autosummary/base.rst | 5 + docs/source/_templates/autosummary/class.rst | 15 ++ docs/source/_templates/autosummary/module.rst | 61 +++++ docs/source/api.rst | 13 +- docs/source/bibliography.rst | 2 + docs/source/conf.py | 250 +++++++++++++++++- docs/source/examples.rst | 2 + docs/source/glossary.rst | 4 +- docs/source/index.rst | 24 +- docs/source/refs.bib | 15 ++ pyproject.toml | 11 +- src/iterative_ensemble_smoother/__init__.py | 32 ++- src/iterative_ensemble_smoother/esmda.py | 7 +- .../esmda_inversion.py | 23 +- src/iterative_ensemble_smoother/sies.py | 119 ++++++--- src/iterative_ensemble_smoother/utils.py | 11 +- 20 files changed, 568 insertions(+), 126 deletions(-) create mode 100644 docs/Makefile create mode 100644 docs/make.bat create mode 100644 docs/source/_templates/autosummary/base.rst create mode 100644 docs/source/_templates/autosummary/class.rst create mode 100644 docs/source/_templates/autosummary/module.rst 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..c8c3fb2c 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,21 @@ 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 +34,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 +48,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/Makefile b/docs/Makefile new file mode 100644 index 00000000..1001f0b4 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHINXBUILD = python -msphinx +SPHINXPROJ = iteratrive_ensemble_smoother +SOURCEDIR = source +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/make.bat b/docs/make.bat new file mode 100644 index 00000000..b1289356 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,36 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=python -msphinx +) +set SOURCEDIR=source +set BUILDDIR=build +set SPHINXPROJ=iterative_ensemble_smoother + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The Sphinx module was not found. Make sure you have Sphinx installed, + echo.then set the SPHINXBUILD environment variable to point to the full + echo.path of the 'sphinx-build' executable. Alternatively you may add the + echo.Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/docs/source/_templates/autosummary/base.rst b/docs/source/_templates/autosummary/base.rst new file mode 100644 index 00000000..b7556ebf --- /dev/null +++ b/docs/source/_templates/autosummary/base.rst @@ -0,0 +1,5 @@ +{{ 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..a2951b2b --- /dev/null +++ b/docs/source/_templates/autosummary/class.rst @@ -0,0 +1,15 @@ +{{ fullname | escape | underline}} + +.. currentmodule:: {{ module }} + +.. autoclass:: {{ objname }} + :members: <-- add at least this line + :private-members: + :show-inheritance: <-- plus I want to show inheritance... + :inherited-members: <-- ...and inherited members too + + {% 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..a3397aac --- /dev/null +++ b/docs/source/_templates/autosummary/module.rst @@ -0,0 +1,61 @@ +{{ 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..d5bf677c 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,26 +1,258 @@ +#!/usr/bin/env python +# +# pyesmda documentation build configuration file, created by +# sphinx-quickstart on Fri Jun 9 13:47:02 2017. +# +# This file is execfile()d with the current directory set to its +# containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +# If extensions (or modules to document with autodoc) are in another +# directory, add these directories to sys.path here. If the directory is +# relative to the documentation root, use os.path.abspath to make it +# absolute, like shown here. + +import datetime + +# FIX: https://github.com/mgaitan/sphinxcontrib-mermaid/issues/72 +import errno +import os +import sys from subprocess import check_output +import sphinx.util.osutil +from sphinx.ext.napoleon.docstring import GoogleDocstring + +import iterative_ensemble_smoother as ies + +sphinx.util.osutil.ENOENT = errno.ENOENT + + +package_path = os.path.abspath("..") +sys.path.insert(0, package_path) + + +def skip(app, what, name, obj, skip, options): + if name in ["__call__"]: + return False + return skip + + +def setup(app): + app.connect("autodoc-skip-member", skip) + + project = "iterative_ensemble_smoother" -copyright = "2022, Equinor" author = "Equinor" -release = "0.1.1" - +copyright = f"2022-{datetime.datetime.today().year}, {author}" +release = ies.__version__ +version = ies.__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.todo", # Support for todo items + "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", + "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 + "jupyter_sphinx", "numpydoc", + "m2r2", # compatibility with markdown + "myst_nb", # to include jupyter nteboo ] -bibtex_bibfiles = ["refs.bib"] + +suppress_warnings = [ + "nbsphinx", +] + +# ----------------------------------------------------------------------------- +# 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" + +# There are two options for replacing |today|: either, you set today to some +# non-false value, then it is used: +# today = '' +# Else, today_fmt is used as the format for a strftime call. +# today_fmt = '%B %d, %Y' + +# 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 reST default role (used for this markup: `text`) to use for all +# documents. +# default_role = None + +# If true, '()' will be appended to :func: etc. cross-reference text. +# add_function_parentheses = True + +# If true, the current module name will be prepended to all description +# unit titles (such as .. function::). +# add_module_names = True + +# If true, sectionauthor and moduleauthor directives will be shown in the +# output. They are ignored by default. +# show_authors = False + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = "sphinx" + +# A list of ignored prefixes for module index sorting. +# modindex_common_prefix = [] + +# If true, keep warnings as "system message" paragraphs in the built documents. +# keep_warnings = False + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = True # set False by default, enable for debugging + + +# -- Options for 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" + +# 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", + # }, + ], +} + +# 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 + +# Issue with attributes section, see: +# https://github.com/sphinx-doc/sphinx/issues/2115 +# Solution: +# https://michaelgoerz.net/notes/extending-sphinx-napoleon-docstring-sections.html +# -- Extensions to the Napoleon GoogleDocstring class --------------------- + + +def parse_keys_section(self, section): + return self._format_fields("Keys", self._consume_fields()) + + +GoogleDocstring._parse_keys_section = parse_keys_section + + +def parse_attributes_section(self, section): + return self._format_fields("Attributes", self._consume_fields()) + + +GoogleDocstring._parse_attributes_section = parse_attributes_section + + +def parse_class_attributes_section(self, section): + return self._format_fields("Class Attributes", self._consume_fields()) + + +GoogleDocstring._parse_class_attributes_section = parse_class_attributes_section + +# we now patch the parse method to guarantee that the the above methods are +# assigned to the _section dict + + +def patched_parse(self): + self._sections["keys"] = self._parse_keys_section + self._sections["class attributes"] = self._parse_class_attributes_section + self._unpatched_parse() + + +GoogleDocstring._unpatched_parse = GoogleDocstring._parse +GoogleDocstring._parse = patched_parse 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..f2fd76b0 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -1,6 +1,4 @@ -=========================== -Iterative Ensemble Smoother -=========================== +.. mdinclude:: ../../README.md .. toctree:: :hidden: @@ -12,26 +10,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 90c46048..19242474 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 ---------- @@ -82,6 +84,7 @@ def __init__( seed: Union[np.random._generator.Generator, int, None] = None, inversion: str = "exact", ) -> None: + """Initialize the instance.""" # Validate inputs if not (isinstance(C_D, np.ndarray) and C_D.ndim in (1, 2)): raise TypeError("Argument `C_D` must be a NumPy array of dimension 1 or 2.") 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..95b365ba 100644 --- a/src/iterative_ensemble_smoother/sies.py +++ b/src/iterative_ensemble_smoother/sies.py @@ -17,50 +17,61 @@ class SIES: - """ - Initialize a Subspace Iterative Ensemble Smoother (SIES) instance. + 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`. - 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 + Note + ---- + The class attributes follow the paper nomenclature. - Parameters + Attributes ---------- - parameters : npt.NDArray[np.double] + X : 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] + A : npt.NDArray[np.double] + A 2D anomaly matrix of shape (num_parameters, ensemble_size) defined as + + .. math:: + + \mathbf{A} = \mathbf{X}\left(\mathbf{I_{N_{e}}} + - \dfrac{1}{N_{e}} \mathbf{11}^{T} + \right) / \sqrt{N_{e}-1} + + with :math:`N_{e}` the number of members in the ensemble. + Note that :math:`\mathbf{A}\mathbf{A}^{T} = \mathbf{C}_{xx}`, + the adjusted parameters empirical covariance matrix. + W : npt.NDArray[np.double] + Current weight matrix W, which represents the current X_i as a linear + combination of the prior. See equation (17) and line 9 in Algorithm 1. + C_dd : 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 + It 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`. - + C_dd_cholesky : npt.NDArray[np.double] + :math:`\mathbf{C}_{dd}` cholesky factorization lower part. + If :math:`\mathbf{C}_{dd}` is a 1D array of diagonal covariances, then it is + simply the square roots of these variances (standard deviations of observation + errors). + d : npt.NDArray[np.double] + 1D array of observations. + inversion : Callable + Inversion routine used in the algorithm. 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. + rng: np.random.RandomState + Pseudorandom number generator state used to generate samples. + """ inversion_funcs = { @@ -78,7 +89,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 +257,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 +293,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 +375,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 --------