From e019aada49d7987140d43cbb790f0b73534b1a83 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. --- .gitignore | 3 + 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 | 16 +- docs/source/bibliography.rst | 2 + docs/source/conf.py | 252 +++++++++++++++++- 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 | 25 +- .../_iterative_ensemble_smoother.py | 123 +++++++-- src/iterative_ensemble_smoother/esmda.py | 3 +- .../esmda_inversion.py | 23 +- 18 files changed, 556 insertions(+), 84 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/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 d9cc0f67..c2896973 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -1,13 +1,11 @@ -API reference -============= +.. _api_ref: -.. currentmodule:: iterative_ensemble_smoother +============= +API Reference +============= -.. autoclass:: iterative_ensemble_smoother.SIES - :members: +.. automodule:: iterative_ensemble_smoother -.. autoclass:: iterative_ensemble_smoother.ES - :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..ac704b6e 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -1,26 +1,260 @@ +#!/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 --------------------- + +# first, we define new methods for any new sections and add them to the 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 acd7d338..88a960e5 100644 --- a/src/iterative_ensemble_smoother/__init__.py +++ b/src/iterative_ensemble_smoother/__init__.py @@ -1,7 +1,24 @@ -""" 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 + + ES + SIES + ESMDA + """ try: from ._version import version as __version__ # type: ignore diff --git a/src/iterative_ensemble_smoother/_iterative_ensemble_smoother.py b/src/iterative_ensemble_smoother/_iterative_ensemble_smoother.py index 4a4771cb..c66e622e 100644 --- a/src/iterative_ensemble_smoother/_iterative_ensemble_smoother.py +++ b/src/iterative_ensemble_smoother/_iterative_ensemble_smoother.py @@ -18,27 +18,28 @@ class SIES: """ - Initialize a Subspace Iterative Ensemble Smoother (SIES) instance. + Subspace Iterative Ensemble Smoother (SIES). 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 + and Reservoir History Matching + written by :cite:t:`evensen2019efficient`. The default step length is described in equation (49) in the paper Formulating the history matching problem with consistent error statistics - written by Geir Evensen (2021), - URL: https://link.springer.com/article/10.1007/s10596-021-10032-7 + written by :cite:t:`evensen2021formulating`. - Parameters + Attributes ---------- + iteration: int + Current iteration number (starts at 1). steplength_schedule : Optional[Callable[[int], float]], optional A function that takes as input the iteration number (starting at 1) and returns steplength (a float in the range (0, 1]). The default is None, which defaults to using an exponential decay. See the references or the function `steplength_exponential`. - seed : Optional[int], 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. Examples -------- @@ -53,7 +54,20 @@ def __init__( *, steplength_schedule: Optional[Callable[[int], float]] = None, seed: Optional[int] = None, - ): + ) -> None: + """ + Initialize the instance. + + Parameters + ---------- + steplength_schedule : Optional[Callable[[int], float]], optional + A function that takes as input the iteration number (starting at 1) and + returns steplength (a float in the range (0, 1]). + The default is None, which defaults to using an exponential decay. + See the references or the function `steplength_exponential`. + seed : Optional[int], optional + Integer used to seed the random number generator. The default is None. + """ self.iteration = 1 self.steplength_schedule = steplength_schedule self.rng = np.random.default_rng(seed) @@ -106,7 +120,9 @@ def fit( inversion: str = "exact", param_ensemble: Optional[npt.NDArray[np.double]] = None, ) -> None: - """Perform one Gauss-Newton step and update the coefficient matrix W. + r""" + Perform one Gauss-Newton step and update the coefficient matrix W. + To apply the coefficient matrix W to the ensemble, call update() after fit(). Parameters @@ -114,16 +130,16 @@ def fit( response_ensemble : npt.NDArray[np.double] A 2D array of responses from the model g(X) of shape (observations, ensemble_size). - This matrix is Y in Evensen (2019). + This matrix is Y in :cite:t:`evensen2019efficient`. observation_errors : npt.NDArray[np.double] Either a 1D array of standard deviations, or a 2D covariance matrix. - 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 e is multivariate normal with covariance - or standard devaiations given by observation_errors. + This is C_dd in :cite:t:`evensen2019efficient`, 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 e is multivariate normal + with covariance or standard devaiations given by observation_errors. observation_values : npt.NDArray[np.double] A 1D array of observations, with shape (observations,). - This is d in Evensen (2019). + This is d in :cite:t:`evensen2019efficient`. truncation : float, optional A value in the range [0, 1], used to determine the number of significant singular values. The default is 0.98. @@ -135,9 +151,12 @@ def fit( that are are active. The default is None, which means all realizations are active. Inactive realizations are ignored (not updated). Must be of shape (ensemble_size,). - inversion : InversionType, optional - The type of subspace inversion used in the algorithm. - The default is InversionType.EXACT. + inversion : str, optional + The type of subspace inversion used in the algorithm. To choose between + \"naive\", "\exact\", \"exact_r\" and \"subspace_re\". For large-scale + problems, the two last options are preferred as they rely on svd for matrix + inversion and uses rescaling to avoid loss of information during the + truncation of small singular values. The default is \"exact\". param_ensemble : Optional[npt.NDArray[np.double]], optional All parameters input to dynamical model used to generate responses. Must be passed if total number of parameters is less than @@ -267,7 +286,8 @@ def fit( self.coefficient_matrix[np.ix_(ensemble_mask, ensemble_mask)] = W def update(self, param_ensemble: npt.NDArray[np.double]) -> npt.NDArray[np.double]: - """Update the parameter ensemble (X in Evensen (2019)). + """ + Update the parameter ensemble (X in :cite:t:`evensen2019efficient`). Parameters ---------- @@ -298,10 +318,20 @@ def __repr__(self) -> str: class ES: + """Implement an Ensemble smoother.""" + def __init__( self, seed: Optional[int] = None, ) -> None: + """ + Initialize the instance. + + Parameters + ---------- + seed : Optional[int], optional + _description_, by default None + """ self.smoother: Optional[SIES] = None self.seed = seed @@ -315,6 +345,45 @@ def fit( inversion: str = "exact", param_ensemble: Optional[npt.NDArray[np.double]] = None, ) -> None: + r""" + Perform one Gauss-Newton step and update the coefficient matrix W. + + To apply the coefficient matrix W to the ensemble, call update() after fit(). + + Parameters + ---------- + response_ensemble : npt.NDArray[np.double] + A 2D array of responses from the model g(X) of shape + (observations, ensemble_size). + This matrix is Y in :cite:t:`evensen2019efficient`. + observation_errors : npt.NDArray[np.double] + Either a 1D array of standard deviations, or a 2D covariance matrix. + This is C_dd in :cite:t:`evensen2019efficient`, 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 e is multivariate normal + with covariance or standard devaiations given by observation_errors. + observation_values : npt.NDArray[np.double] + A 1D array of observations, with shape (observations,). + This is d in :cite:t:`evensen2019efficient`. + truncation : float, optional + A value in the range [0, 1], used to determine the number of + significant singular values. The default is 0.98. + step_length : Optional[float], optional + If given, this value (in the range [0, 1]) overrides the step length + schedule that was provided at initialization. The default is None. + inversion : str, optional + The type of subspace inversion used in the algorithm. To choose between + \"naive\", "\exact\", \"exact_r\" and \"subspace_re\". For large-scale + problems, the two last options are preferred as they rely on svd for matrix + inversion and uses rescaling to avoid loss of information during the + truncation of small singular values. The default is \"exact\". + param_ensemble : Optional[npt.NDArray[np.double]], optional + All parameters input to dynamical model used to generate responses. + Must be passed if total number of parameters is less than + ensemble_size - 1 and the dynamical model g(x) is non-linear. + This is X in Evensen (2019), and has shape (parameters, ensemble_size). + The default is None. See section 2.4.3 in Evensen (2019). + """ self.smoother = SIES(seed=self.seed) self.smoother.fit( response_ensemble, @@ -328,6 +397,20 @@ def fit( ) def update(self, param_ensemble: npt.NDArray[np.double]) -> npt.NDArray[np.double]: + """ + Update the parameter ensemble (X in :cite:t:`evensen2019efficient`). + + Parameters + ---------- + param_ensemble : npt.NDArray[np.double] + The (prior) parameter ensemble. The same `param_ensemble` should be + used in each Gauss-Newton step. + + Returns + ------- + np.ndarray + Updated parameter ensemble. + """ assert self.smoother is not None return self.smoother.update(param_ensemble) diff --git a/src/iterative_ensemble_smoother/esmda.py b/src/iterative_ensemble_smoother/esmda.py index 90c46048..9f9b65aa 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 @@ -38,7 +39,7 @@ class ESMDA: """Initialize 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 ---------- 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