diff --git a/README.md b/README.md index d086b17..6b4625e 100644 --- a/README.md +++ b/README.md @@ -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. @@ -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) diff --git a/dass/analysis.py b/dass/analysis.py index 3c424e7..2272d03 100644 --- a/dass/analysis.py +++ b/dass/analysis.py @@ -1,6 +1,3 @@ -"""Collection of analysis steps like ES, IES etc. -""" - import numpy as np import numpy.typing as npt @@ -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 diff --git a/notebooks/ES_2D_Heat_Equation.py b/notebooks/ES_2D_Heat_Equation.py index 4ce75bd..0a2aaf5 100644 --- a/notebooks/ES_2D_Heat_Equation.py +++ b/notebooks/ES_2D_Heat_Equation.py @@ -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 @@ -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() - # %% diff --git a/tests/test_analysis.py b/tests/test_analysis.py index 77038eb..edc75db 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -2,7 +2,7 @@ rng = np.random.default_rng() -from dass.analysis import ES, IES +from dass.analysis import ES """ Bayes' theorem states that @@ -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()