Clean up and document; tests

README.rst

@@ -10,5 +10,6 @@ presented by Emerick and Reynolds, 2013.
 - **Bug reports:** https://github.com/tuda-geo/resmda/issues
 
 
-Available through pip: ``pip install resmda``.
+Available through pip and conda:
+``pip install resmda`` / ``conda install -c conda-forge resmda``.
 

docs/index.rst

@@ -68,7 +68,7 @@ resmda Documentation
         ^^^^^^^^^^^^^
 
         The API reference of resmda includes almost every function and class.
-        Some of the underlying theory is also described in the docstring.
+        Some of the underlying theory is also described in the docstrings.
 
         +++
 

docs/manual/about.rst

@@ -10,15 +10,15 @@ ES-MDA
 ------
 
 
-In the following an introduction to the ES-MDA (Ensemble Smoother with Multiple
-Data Assimilation) algorithm following [EmRe13]_:
+In the following an introduction to the Ensemble Smoother with Multiple Data
+Assimilation (ES-MDA) algorithm following [EmRe13]_:
 
 In history-matching problems, it is common to consider solely the
 parameter-estimation problem and thereby neglecting model uncertainties. Thus,
-unlike EnKF, the parameters and states are always consistent (Thulin et al.,
-2007). This fact helps to explain the better data matches obtained by ES-MDA
-compared to EnKF. The analyzed vector of model parameters :math:`m^a` is given
-in that case by
+unlike with the ensemble Kalman filter (EnKF), the parameters and states are
+always consistent (Thulin et al., 2007). This fact helps to explain the better
+data matches obtained by ES-MDA compared to EnKF. The analyzed vector of model
+parameters :math:`m^a` is given in that case by
 
 .. math::
     m_j^a = m_j^f + C_\text{MD}^f \left(C_\text{DD}^f + \alpha C_\text{D}

docs/manual/installation.rst

@@ -7,6 +7,12 @@ You can install the latest release of resmda simply via ``pip``
 
    pip install resmda
 
+and ``conda``
+
+.. code-block:: console
+
+   conda install -c conda-forge resmda
+
 or clone the repository and install it manually with
 
 .. code-block:: console

docs/manual/license.rst

@@ -3,8 +3,8 @@
 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.
+Everything of resmda is open-source: The code is released under the Apache 2.0
+license, the text (e.g., the manual) is released under the CC-BY 4.0 license.
 
 
 License of Text and Website

examples/fluvialreservoir.py

@@ -161,7 +161,7 @@ def sim(x):
 nd_positions = np.tile(np.array([ox, oy]), time.size).T
 
 # Create matrix
-loc_mat = resmda.build_localization_matrix(RP.cov, nd_positions, (nx, ny))
+loc_mat = resmda.localization_matrix(RP.cov, nd_positions, (nx, ny))
 
 
 ###############################################################################

examples/localization.py

@@ -96,7 +96,7 @@ def sim(x):
 nd_positions = np.tile(np.array([ox, oy]), time.size).T
 
 # Create matrix
-loc_mat = resmda.build_localization_matrix(RP.cov, nd_positions, (nx, ny))
+loc_mat = resmda.localization_matrix(RP.cov, nd_positions, (nx, ny))
 
 # QC localization matrix
 fig, ax = plt.subplots(1, 1, constrained_layout=True)

requirements-dev.txt

@@ -16,6 +16,7 @@ sphinx_gallery>=0.16
 pydata_sphinx_theme
 sphinx_automodapi
 ipykernel
+ipympl
 pooch
 
 # FOR TESTING

resmda/__init__.py

@@ -15,8 +15,8 @@
 # the License.
 
 from resmda import utils
-from resmda.utils import Report
-from resmda.data_assimilation import esmda, build_localization_matrix
+from resmda.utils import localization_matrix, Report
+from resmda.data_assimilation import esmda
 from resmda.reservoir_simulator import Simulator, RandomPermeability
 
 
@@ -25,8 +25,7 @@
 
 
 __all__ = ['reservoir_simulator', 'data_assimilation', 'utils',
-           'Simulator', 'RandomPermeability',
-           'esmda', 'build_localization_matrix',
-           'rng', 'Report']
+           'esmda', 'Simulator', 'RandomPermeability',
+           'localization_matrix', 'rng', 'Report']
 
 __version__ = utils.__version__

resmda/data_assimilation.py

@@ -16,9 +16,9 @@
 
 import numpy as np
 
-from resmda.utils import rng
+from resmda import utils
 
-__all__ = ['esmda', 'build_localization_matrix']
+__all__ = ['esmda']
 
 
 def __dir__():
@@ -47,23 +47,24 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
         Observed data of shape ``(nd)``.
     sigma : {float, ndarray}
         Standard deviation(s) of the observation noise.
-    alphas : {int, ndarray}, default: 4
+    alphas : {int, array-like}, default: 4
         Inflation factors for ES-MDA.
     data_prior : ndarray, default: None
-        Prior data ensemble, of shape (ne, nd).
+        Prior data ensemble, of shape ``(ne, nd)``.
     callback_post : function, default: None
         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 also ``forward(model_post)``.
     return_steps : bool, default: False
-        If true, returns model (and data) of all ES-MDA steps.
-        Setting ``return_steps`` to True enforces ``return_post_data=True``.
+        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.
+        Seed or random generator for reproducibility; see
+        :func:`resmda.utils.rng`.
     localization_matrix : {ndarray, None}, default: None
         If provided, apply localization to the Kalman gain matrix, of shape
-        (model-shape, nd).
+        ``(model-shape, nd)``.
 
 
     Returns
@@ -78,6 +79,9 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
     ne = model_prior.shape[0]
     nd = data_obs.size
 
+    # Get random number generator
+    rng = utils.rng(random)
+
     # Expand sigma if float
     if np.asarray(sigma).size == 1:
         sigma = np.zeros(nd) + sigma
@@ -106,7 +110,7 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
         # For each ensemble member, perturb the observation vector using
         # d_uc = d_obs + sqrt(α_i) * C_D^0.5 z_d; z_d ~ N(0, I_N_d)
 
-        zd = rng(random).normal(size=(ne, nd))
+        zd = rng.normal(size=(ne, nd))
         data_pert = data_obs + np.sqrt(alpha) * sigma * zd
 
         # == Step (c) of Emerick & Reynolds, 2013 ==
@@ -166,37 +170,3 @@ def esmda(model_prior, forward, data_obs, sigma, alphas=4, data_prior=None,
         return model_post, data_post
     else:
         return model_post
-
-
-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
-    data positions.
-
-    Parameters
-    ----------
-    cov_matrix : ndarray
-        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,
-        of size (nd, 2).
-    shape : tuple
-        Dimensions of the grid (nx, ny).
-
-    Returns
-    -------
-    loc_matrix : ndarray
-        Localization matrix of shape (nx, ny, nd).
-
-    """
-    # Convert 2D positions of data points to 1D indices suitable for accessing
-    # the covariance matrix
-    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, and reshape.
-    return cov_matrix[:, indices.astype(int)].reshape((*shape, -1), order='F')