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