Localization II (#7)

docs/_static/style.css

@@ -0,0 +1,33 @@
+/* "emsig"-colours
+ blue         : #447fc1
+ blue-lighter : #61a2d8
+ blue-darker  : #285985
+ red          : #ef413d
+*/
+
+/* Button colours */
+
+.sd-btn-info {
+  font-weight: bold;
+  background: #447fc1 !important;
+  border-color: #447fc1 !important;
+}
+
+.sd-btn-info:hover {
+  background: #ef413d !important;
+  border-color: #ef413d !important;
+}
+
+/* Card headers */
+.sd-card-header {
+  font-size: var(--pst-font-size-h5);
+  font-weight: bold;
+  padding: 1rem 2rem;
+  text-align: center;
+}
+
+/* fa stuff */
+.fa {
+  margin: 0 1rem 0 0;
+  color: #ef413d;
+}

docs/conf.py

@@ -102,6 +102,7 @@
     "logo": {
         "text": "resmda",
     },
+    'navigation_with_keys': True,
     "github_url": "https://github.com/tuda-geo/resmda",
     # "use_edit_page_button": True,
 }
@@ -116,3 +117,8 @@
 html_use_modindex = True
 html_file_suffix = '.html'
 htmlhelp_basename = 'resmda'
+html_css_files = [
+    "style.css",
+    "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/" +
+    "css/font-awesome.min.css"
+]

docs/index.rst

@@ -18,11 +18,63 @@ resmda Documentation
    api/index
 
 
-**TODOs**
+.. grid:: 1 1 3 3
+    :gutter: 2
 
-- Add tests
-- Finish docstrings & code documentation
-- Compare with existing esmda's
-- docs:
+    .. grid-item-card::
 
-  - only deploy for versions
+        :fa:`book;fa-4x`
+
+        User Guide
+        ^^^^^^^^^^
+
+        The manual contains installation instructions, some background theory,
+        important references, license information, and more.
+
+        +++
+
+        .. button-ref:: manual
+            :expand:
+            :color: info
+            :click-parent:
+
+            To the user guide
+
+    .. grid-item-card::
+
+        :fa:`image;fa-4x`
+
+        Gallery
+        ^^^^^^^
+
+        The gallery contains examples on the usage of resmda, and is generally
+        the best way to get started. Download them and modify them to your
+        needs!
+
+        +++
+
+        .. button-ref:: resmda_gallery
+            :expand:
+            :color: info
+            :click-parent:
+
+            To the gallery
+
+    .. grid-item-card::
+
+        :fa:`code;fa-4x`
+
+        API reference
+        ^^^^^^^^^^^^^
+
+        The API reference of resmda includes almost every function and class.
+        Some of the underlying theory is also described in the docstring.
+
+        +++
+
+        .. button-ref:: api
+            :expand:
+            :color: info
+            :click-parent:
+
+            To the API reference

docs/manual/about.rst

@@ -4,9 +4,12 @@ About
 A simple 2D reservoir simulator and a straight-forward implementation of the
 basic *Ensemble Smoother with Multiple Data Assimilation* (ES-MDA) algorithm.
 
+.. _esmda:
+
 ES-MDA
 ------
 
+
 In the following an introduction to the ES-MDA (Ensemble Smoother with Multiple
 Data Assimilation) algorithm following [EmRe13]_:
 

docs/manual/index.rst

@@ -16,4 +16,5 @@ User Guide
    about
    installation
    references
+   license
 

docs/manual/installation.rst

@@ -1,8 +1,16 @@
 Installation
 ============
 
-You can install resmda simply via ``pip``:
+You can install the latest release of resmda simply via ``pip``:
 
 .. code-block:: console
 
    pip install resmda
+
+or clone the repository and run within the command
+
+.. code-block:: console
+
+   make install
+
+to get the latest version.

docs/manual/license.rst

@@ -0,0 +1,34 @@
+.. _license:
+
+License
+=======
+
+Everything of resmda is open-source: Code is released under the Apache 2.0
+license, text (e.g., the manual) is released under the CC-BY 4.0 license.
+
+
+License of Text and Website
+---------------------------
+
+This work is licensed under Attribution 4.0 International (CC-BY-4.0).
+
+    https://creativecommons.org/licenses/by/4.0/
+
+
+License of Code
+---------------
+
+Copyright 2024 Dieter Werthmüller, Gabriel Serrao Seabra
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    https://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+

examples/README.rst

@@ -1,4 +1,8 @@
-.. _sphx_glr_gallery:
+.. _resmda_gallery:
 
+=======
 Gallery
 =======
+
+:Release: |version|
+:Date: |today|

examples/basicESMDA.py

@@ -83,7 +83,7 @@ def plot_result(mpost, dpost, dobs, title, ylim):
 
     # Plot Likelihood
     ax2.plot(
-        *pseudopdf(resmda.utils.rng().normal(dobs, size=(ne, dobs.size))),
+        *pseudopdf(resmda.rng.normal(dobs, size=(ne, dobs.size))),
         'C2', lw=2, label='Datum'
     )
 
@@ -133,7 +133,7 @@ def plot_result(mpost, dpost, dobs, title, ylim):
 obs_std = 1.0
 
 # Prior: Let's start with ones as our prior guess
-mprior = resmda.utils.rng().normal(loc=1.0, scale=obs_std, size=(ne, 1))
+mprior = resmda.rng.normal(loc=1.0, scale=obs_std, size=(ne, 1))
 
 
 ###############################################################################

examples/basicreservoir.py

@@ -2,24 +2,16 @@
 2D Reservoir ESMDA example
 ==========================
 
-Data Assimilation using ESMDA in Reservoir Simulation
------------------------------------------------------
+Ensemble Smoother Multiple Data Assimilation (ES-MDA) in Reservoir Simulation.
 
-*Advanced Data Assimilation using Ensemble Smoother Multiple Data Assimilation
-(ESMDA) in Reservoir Simulation.*
-
-.. math::
-    m_j^a = m_j^f + C_{MD}^f (C_{DD}^f + \alpha_i C_D)^{-1} (d_{uc,j} -
-    d_j^f);\quad \text{for} \quad j=1, 2, \dots, N_e
-
-- Prior model: M (Ne, Nx, Ny)
-- Prior data: D (Ne, Nt)
 """
 import numpy as np
 import matplotlib.pyplot as plt
 
 import resmda
 
+# For reproducibility, we instantiate a random number generator with a fixed
+# seed. For production, remove the seed!
 rng = np.random.default_rng(1848)
 
 # sphinx_gallery_thumbnail_number = 4

resmda/__init__.py

@@ -20,9 +20,13 @@
 from resmda.reservoir_simulator import Simulator, RandomPermeability
 
 
+# Initialize a random number generator.
+rng = utils.rng()
+
+
 __all__ = ['reservoir_simulator', 'data_assimilation', 'utils',
            'Simulator', 'RandomPermeability',
            'esmda', 'build_localization_matrix',
-           'Report']
+           'rng', 'Report']
 
 __version__ = utils.__version__

resmda/data_assimilation.py

@@ -18,20 +18,20 @@
 
 from resmda.utils import rng
 
-__all__ = ['esmda']
+__all__ = ['esmda', 'build_localization_matrix']
 
 
 def __dir__():
     return __all__
 
 
 def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
-          callback_post=None, return_post_data=True, return_steps=False,
-          random=None, localization_matrix=None):
-    """ESMDA algorithm (Emerick and Reynolds, 2013) with optional localization.
-
-    ES-MDA as presented by [EmRe13]_.
+          localization_matrix=None, callback_post=None, return_post_data=True,
+          return_steps=False, random=None):
+    """ES-MDA algorithm ([EmRe13]_) with optional localization.
 
+    Consult the section :ref:`esmda` in the manual for the theory and more
+    information about ES-MDA.
 
     Parameters
     ----------
@@ -52,16 +52,18 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
     data_prior : ndarray, default: None
         Prior data ensemble, of shape (ne, nd).
     callback_post : function, default: None
-        Function to be executed after each ESMDA iteration.
+        Function to be executed after each ES-MDA iteration to the posterior
+        model, ``callback_post(model_post)``.
     return_post_data : bool, default: True
-        If true, returns data
+        If true, returns also ``forward(model_post)``.
     return_steps : bool, default: False
-        If true, returns model (and data) of all ESMDA steps.
-        Setting ``return_steps`` to True wil enforce ``return_post_data``.
+        If true, returns model (and data) of all ES-MDA steps.
+        Setting ``return_steps`` to True enforces ``return_post_data=True``.
     random : {None, int,  np.random.Generator}, default: None
         Seed or random generator for reproducibility.
     localization_matrix : {ndarray, None}, default: None
-        If provided, apply localization to the Kalman gain matrix.
+        If provided, apply localization to the Kalman gain matrix, of shape
+        (model-shape, nd).
 
 
     Returns
@@ -91,7 +93,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
 
     # Loop over alphas
     for i, alpha in enumerate(alphas):
-        print(f"ESMDA step {i+1: 3d}; α={alpha}")
+        print(f"ES-MDA step {i+1: 3d}; α={alpha}")
 
         # == Step (a) of Emerick & Reynolds, 2013 ==
         # Run the ensemble from time zero.
@@ -132,7 +134,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
 
         # Apply localization if provided
         if localization_matrix is not None:
-            K *= localization_matrix[..., None]
+            K *= localization_matrix
 
         # Update the ensemble parameters
         model_post += np.moveaxis(K @ (data_pert - data_prior).T, -1, 0)
@@ -166,29 +168,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
         return model_post
 
 
-def convert_positions_to_indices(positions, grid_dimensions):
-    """Convert 2D grid positions to 1D indices assuming zero-indexed positions.
-
-    Parameters
-    ----------
-    positions : ndarray
-        Array of (x, y) positions in the grid.
-    grid_dimensions : tuple
-        Dimensions of the grid (nx, ny).
-
-    Returns
-    -------
-    indices : ndarray
-        Array of indices corresponding to the positions in a flattened array.
-
-    """
-    nx, ny = grid_dimensions
-    # Ensure the positions are zero-indexed and correctly calculated for
-    # row-major order
-    return positions[:, 1] * nx + positions[:, 0]
-
-
-def build_localization_matrix(cov_matrix, data_positions, grid_dimensions):
+def build_localization_matrix(cov_matrix, data_positions, shape):
     """Build a localization matrix
 
     Build a localization matrix from a full covariance matrix based on specific
@@ -197,22 +177,26 @@ def build_localization_matrix(cov_matrix, data_positions, grid_dimensions):
     Parameters
     ----------
     cov_matrix : ndarray
-        The full (nx*ny) x (nx*ny) covariance matrix.
+        The lower triangular covariance matrix (nx*ny, nx*ny).
     data_positions : ndarray
-        Positions in the grid for each data point (e.g., wells), zero-indexed.
-    grid_dimensions : tuple
+        Positions in the grid for each data point (e.g., wells), zero-indexed,
+        of size (nd, 2).
+    shape : tuple
         Dimensions of the grid (nx, ny).
 
     Returns
     -------
     loc_matrix : ndarray
-        Localization matrix of shape (nx*ny, number of data positions).
+        Localization matrix of shape (nx, ny, nd).
 
     """
     # Convert 2D positions of data points to 1D indices suitable for accessing
     # the covariance matrix
-    indices = convert_positions_to_indices(
-            data_positions, grid_dimensions).astype(int)
+    indices = data_positions[:, 1] * shape[0] + data_positions[:, 0]
+
+    # Create full matrix from lower triangular matrix
+    cov_matrix = cov_matrix + np.tril(cov_matrix, -1).T
+
     # Extract the columns from the covariance matrix corresponding to each data
-    # point's position
-    return cov_matrix[:, indices]
+    # point's position, and reshape.
+    return cov_matrix[:, indices.astype(int)].reshape((*shape, -1), order='F')