Skip to content

Commit

Permalink
Remove analysis.IES
Browse files Browse the repository at this point in the history
IES keeps state and is not as easy to understand as ES.
It is also not suitable for adaptive localization.
  • Loading branch information
dafeda committed Dec 25, 2023
1 parent 6968589 commit 8bbfbb7
Show file tree
Hide file tree
Showing 4 changed files with 4 additions and 127 deletions.
21 changes: 3 additions & 18 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,8 @@
`dass` is tool for learning about data assimilation / history matching created by the developers of [ERT](https://github.com/equinor/ert).
It is inspired by [DAPPER](https://github.com/nansencenter/DAPPER) and [HistoryMatching](https://github.com/patnr/HistoryMatching).

It includes implementations of Ensemble Smoother (ES) as given in [1] and Iterative Ensemble Smoother (IES) as given in [2],
see [dass/analysis.py](dass/analysis.py).
The implementation of ES can easily be extended to the Ensemble Smoother with Multiple Data Assimilation (ES-MDA) as described in [3].
It includes implementations of Ensemble Smoother (ES) as given in [1], see [dass/analysis.py](dass/analysis.py).
The implementation of ES can easily be extended to the Ensemble Smoother with Multiple Data Assimilation (ES-MDA) as described in [2].

For notebooks with examples and tutorials see the `notebooks/` folder.

Expand All @@ -30,23 +29,9 @@ jupyter notebook
# To make sure everything works, run on the of the notebooks in the notebooks/ folder.
```

## On notation

The implementation of ES is based on section 9.5 of [1], while the implementation of IES is based on [2].
The notation used in the two papers differ slightly, so we have made a few tweaks to make them more similar.

- $A$ is used for the prior ensemble. (It's $X$ in [2])
- $E$ is not divided by $\sqrt{N-1}$ as is done in [2], which means that we do not multiply $E$ by $\sqrt{N-1}$ in the definition of $E$.
- We do not use $EE^T / (N-1)$ to estimate the parameter covariance matrix, because we assume a diagonal observation error covariance matrix $C_{dd}$.
We instead scale matrices used in the analysis step such that $C_{dd}$ becomes the identity matrix.
This is what is known as exact inversion.
- $Y$ is used to hold measured responses, which are predictions made by the dynamical model at points in time and space for which we have observations.

## References

[1] - [Data Assimilation
The Ensemble Kalman Filter](https://link.springer.com/book/10.1007/978-3-642-03711-5)

[2] - [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)

[3] - [Ensemble smoother with multiple data assimilation](https://www.sciencedirect.com/science/article/pii/S0098300412000994)
[2] - [Ensemble smoother with multiple data assimilation](https://www.sciencedirect.com/science/article/pii/S0098300412000994)
56 changes: 0 additions & 56 deletions dass/analysis.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,3 @@
"""Collection of analysis steps like ES, IES etc.
"""

import numpy as np
import numpy.typing as npt

Expand Down Expand Up @@ -42,56 +39,3 @@ def ES(
X = np.identity(N) + S.T @ np.linalg.inv(C) @ Dprime

return X


def IES(
Y: npt.NDArray[np.float_],
D: npt.NDArray[np.float_],
Cdd: npt.NDArray[np.float_],
W: npt.NDArray[np.float_],
gamma: float,
) -> npt.NDArray[np.float_]:
"""Iterative Ensemble Smoother based on `evensen2019`.
Parameters
----------
Y : npt.NDArray[np.float_]
Measured responses
D : npt.NDArray[np.float_]
Perturbed observations
Cdd : npt.NDArray[np.float_]
Diagonal observation error covariance
W : npt.NDArray[np.float_]
Coefficient matrix
gamma : float
Step length
Returns
-------
npt.NDArray[np.float_]
Array such that prior multiplied by (I + W) yields posterior
"""
N = W.shape[0]
# Line 4 of `Algorithm 1`.
Y_centered = Y - Y.mean(axis=1, keepdims=True)

Omega = np.identity(N) + (W - W.mean(axis=1, keepdims=True))

S = np.linalg.solve(Omega.T, Y_centered.T).T

H = S @ W + D - Y

# Eq. 42 - Inversion in measurement space
# W = W - gamma * (W - S.T @ np.linalg.inv(S @ S.T + (N - 1) * Cdd) @ H)

# Eq. 50 - Exact inversion without assuming a diagonal error covariance
# notice that the inversion is computed in the ensemble subspace.
W = W - gamma * (
W
- np.linalg.inv(S.T @ np.linalg.inv(Cdd) @ S + (N - 1) * np.identity(N))
@ S.T
@ np.linalg.inv(Cdd)
@ H
)

return W
25 changes: 0 additions & 25 deletions notebooks/ES_2D_Heat_Equation.py
Original file line number Diff line number Diff line change
Expand Up @@ -506,25 +506,6 @@ def plot_responses(
err_prior = alpha_t.ravel() - A.mean(axis=1)
np.sqrt(np.mean(err_prior * err_prior))

# %% [markdown]
# ## IES

# %%
# Step length in Gauss Newton
gamma = 1.0

# Coefficient matrix as defined in Eq. 16 and Eq. 17.
W = np.zeros(shape=(N, N))

W = analysis.IES(Y, D, Cdd, W, gamma)
X_IES = np.identity(N) + W
A_IES = A @ X_IES
A_IES = A_IES.clip(min=1e-8)

# %%
err_posterior_ies = alpha_t.ravel() - A_IES.mean(axis=1)
np.sqrt(np.mean(err_posterior_ies * err_posterior_ies))

# %% [markdown]
# ## Graphical comparison of prior and posterior mean-fields

Expand Down Expand Up @@ -613,10 +594,4 @@ def plot_responses(

plot_responses(np.unique(k_levels), d, Y_df, u_t, Y_df_post)

# %% [markdown]
# ## Check that single iteration of IES with step length 1.0 is the same as ES.

# %%
assert np.isclose(X_IES, X, atol=1e-5).all()

# %%
29 changes: 1 addition & 28 deletions tests/test_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

rng = np.random.default_rng()

from dass.analysis import ES, IES
from dass.analysis import ES

"""
Bayes' theorem states that
Expand Down Expand Up @@ -56,30 +56,3 @@ def test_ES_with_localisation():
A_ES = A @ X

assert np.isclose(A_ES_local, A_ES).all()


def test_IES():
# Check that single iteration of IES with step length 1.0 is the same as ES.
gamma = 1.0
W = np.zeros(shape=(N, N))
W = IES(Y, D, Cdd, W, gamma)
X_IES = np.identity(N) + W
A_IES = A @ X_IES

# Potentially flaky, but still useful.
assert np.isclose(A_IES[0, :].mean(), observation / 2, rtol=0.1)
assert np.isclose(A_IES[1, :].mean(), observation / 2, rtol=0.1)
assert (np.abs(np.cov(A_IES) - np.identity(nparam)) < 0.15).all()


def test_ES_vs_IES():
X = ES(Y, D, Cdd)
A_ES = A @ X

gamma = 1.0
W = np.zeros(shape=(N, N))
W = IES(Y, D, Cdd, W, gamma)
X_IES = np.identity(N) + W
A_IES = A @ X_IES

assert np.isclose(A_ES, A_IES).all()

0 comments on commit 8bbfbb7

Please sign in to comment.