diff --git a/notebooks/ES_2D_Heat_Equation.py b/notebooks/ES_2D_Heat_Equation.py index b25ca13..4ce75bd 100644 --- a/notebooks/ES_2D_Heat_Equation.py +++ b/notebooks/ES_2D_Heat_Equation.py @@ -6,7 +6,7 @@ # extension: .py # format_name: percent # format_version: '1.3' -# jupytext_version: 1.14.5 +# jupytext_version: 1.16.0 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -37,8 +37,6 @@ from p_tqdm import p_map -import iterative_ensemble_smoother as ies - # %% # %load_ext autoreload # %autoreload 2 @@ -487,22 +485,6 @@ def plot_responses( # Setting negative values to a small positive value but not zero because we want to be able to divide by them. A_ES = A_ES.clip(min=1e-8) -# %% [markdown] -# ## Testing the new iterative_ensemble_smoother package -# -# As part of ERT development, we wrote an efficient implementation of the iterative ensemble smoother. -# This package is available via pypi and you can easily test it out here if you wish. -# -# ```python -# A_ES_ert = ies.ensemble_smoother_update_step( -# Y, -# A, -# obs_std[~is_outlier], -# obs_value[~is_outlier], -# inversion=ies.InversionType.EXACT, -# ) -# ``` - # %% [markdown] # ## Numerical comparison of prior and posterior using RMSE # diff --git a/notebooks/EnKF_2D_Heat_Equation.py b/notebooks/EnKF_2D_Heat_Equation.py deleted file mode 100644 index 282cc09..0000000 --- a/notebooks/EnKF_2D_Heat_Equation.py +++ /dev/null @@ -1,317 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.14.1 -# kernelspec: -# display_name: Python 3 (ipykernel) -# language: python -# name: python3 -# --- - -# %% [markdown] -# # EnKF using 2D heat equation - -# %% -# %load_ext autoreload - -# %% -import numpy as np - -np.set_printoptions(suppress=True) -rng = np.random.default_rng() - -import matplotlib.pyplot as plt - -plt.rcParams["figure.figsize"] = (6, 6) -plt.rcParams.update({"font.size": 10}) -# Ignore error when drawing many figures -plt.rcParams.update({"figure.max_open_warning": 0}) -from ipywidgets import interact -import ipywidgets as widgets - -from p_tqdm import p_map - -from scipy.ndimage import gaussian_filter - -# %% -# %autoreload 2 -from dass import pde, utils, analysis, taper - -# %% [markdown] -# # EnKF -# -# Constant parameter, unknown field. -# Stop forward model at every point in time where we have observations. - -# %% [markdown] -# ## Define parameters, set true initial conditions and calculate the true temperature field - -# %% -N = 100 - -k_start = 0 -k_end = 300 - -# Number of grid-cells in x and y direction -nx = 10 - -# Set the coefficient of heat transfer for each grid cell. -# Using trick from page 15 of "An Introduction to the Numerics of Flow in Porous Media using Matlab". -# It's a nice way to generate realistic-looking parameter fields. -# In real life we use third-party tools to generate good (whatever that means) prior parameter fields. -alpha_t = np.exp( - 5 - * gaussian_filter(gaussian_filter(rng.random(size=(nx, nx)), sigma=2.0), sigma=1.0) -) - - -dx = 1 -dt = dx**2 / (4 * np.max(alpha_t)) - -# True initial temperature field. -u_top = 100.0 -u_init = np.empty((k_end, nx, nx)) -u_init.fill(0.0) - -# Set the boundary conditions -u_init[:, 1 : (nx - 1), 0] = u_top - -# How much noise to add to heat equation, also called model noise. -# scale = 0.1 -scale = None - -u_t = pde.heat_equation(u_init, alpha_t, dx, dt, k_start, k_end, rng=rng, scale=scale) - - -# %% [markdown] -# ## Interactive plot of true temperature field - - -# %% -def interactive_truth(k): - fig, ax = plt.subplots() - fig.suptitle("True temperature field") - p = ax.pcolormesh(u_t[k].T, cmap=plt.cm.jet) - ax.set_title(f"k = {k}") - ax.invert_yaxis() - utils.colorbar(p) - fig.tight_layout() - - -interact( - interactive_truth, - k=widgets.IntSlider(min=k_start, max=k_end - 1, step=1, value=0), -) - -# %% [markdown] -# ## Define random seeds because multiprocessing -# -# https://numpy.org/doc/stable/reference/random/parallel.html#seedsequence-spawning - -# %% -ss = np.random.SeedSequence(12345) -child_seeds = ss.spawn(N) -streams = [np.random.default_rng(s) for s in child_seeds] - -# %% [markdown] -# ## Define placement of sensors and generate synthetic observations based on the true temperature field - -# %% -# placement of sensors, i.e, where the observations are done -pad = 1 -coords = np.array([(x, y) for x in range(pad, nx - pad) for y in range(pad, nx - pad)]) -ncoords = coords.shape[0] -nmeas = 40 -coords_idx = np.random.choice(np.arange(ncoords), size=nmeas, replace=False) -obs_coordinates = [utils.Coordinate(xc, yc) for xc, yc in coords[coords_idx]] - -# At which times observations are taken -obs_times = np.linspace(5, k_end, 50, endpoint=False, dtype=int) - -d = utils.observations(obs_coordinates, obs_times, u_t, lambda value: abs(0.05 * value)) -# number of measurements -m = d.shape[0] -print("Number of observations: ", m) - -k_levels = d.index.get_level_values("k").to_list() -x_levels = d.index.get_level_values("x").to_list() -y_levels = d.index.get_level_values("y").to_list() - -# Plot temperature field and show placement of sensors. -obs_coordinates_from_index = set(zip(x_levels, y_levels)) -x, y = zip(*obs_coordinates_from_index) - -fig, ax = plt.subplots() -p = ax.pcolormesh(u_t[0].T, cmap=plt.cm.jet) -ax.invert_yaxis() -ax.set_title("True temperature field with sensor placement") -utils.colorbar(p) -ax.plot([i + 0.5 for i in x], [j + 0.5 for j in y], "s", color="white", markersize=5) - - -# %% -def gen_field(): - u = np.empty((k_end, nx, nx)) - - u_top = 100.0 - u_init = np.empty((k_end, nx, nx)) - u_init.fill(0.0) - - # Set the boundary conditions - u_init[:, 1 : (nx - 1), 0] = u_top + rng.normal() - - return u - - -fields = [gen_field() for _ in range(N)] - - -# %% -def matrix_from_fields(fields, k): - nx = fields[0][0].shape[0] - A = np.zeros(shape=(nx * nx, N)) - for f in range(len(fields)): - A[:, f] = fields[f][k].ravel() - return A - - -# %% [markdown] -# ## Plot tapering function used for localisation - -# %% -fig, ax = plt.subplots() -ax.plot(taper.gauss(np.linspace(-nx, nx), 3.0)) - -# %% -A_no_update = {} -localize = True -k_start = 0 -for k_obs in d.index.get_level_values("k").unique().to_list(): - fields = p_map( - pde.heat_equation, - fields, - [alpha_t] * N, - [dx] * N, - [dt] * N, - [k_start] * N, - [k_obs + 1] * N, - streams, - [np.sqrt(dt)] * N, - desc=f"Running forward model from {k_start} to {k_obs}", - ) - - d_k = d.query(f"k == {k_obs}") - m = d_k.shape[0] - - A = matrix_from_fields(fields, k_obs) - A_no_update[k_obs] = A.copy() - - # measure response - Y = np.zeros(shape=(m, N)) - for i in range(N): - Y[:, i] = A[:, i].reshape(nx, nx)[ - d_k.index.get_level_values("x").to_list(), - d_k.index.get_level_values("y").to_list(), - ] - - Cdd = np.diag(d_k.sd.values**2) - - E = rng.multivariate_normal(mean=np.zeros(len(Cdd)), cov=Cdd, size=N).T - E = E - E.mean(axis=1, keepdims=True) - assert E.shape == (m, N) - - D = np.ones((m, N)) * d_k.value.values.reshape(-1, 1) + E - - S = Y - Y.mean(axis=1, keepdims=True) - - A_centered = A - A.mean(axis=0, keepdims=True) - Dprime = D - Y - - if localize: - for i in range(nx**2): - state_idx = np.unravel_index(i, shape=(nx, nx)) - dist = np.sqrt( - (state_idx[0] - obs_coordinates[0].x) ** 2 - + (state_idx[1] - obs_coordinates[0].y) ** 2 - ) - taper_coeff = taper.gauss(dist, 2.0) - - K = A_centered[i, :] @ S.T @ np.linalg.pinv(S @ S.T + (N - 1) * Cdd) - - # K = ( - # A_centered[i, :] - # @ S.T - # @ np.linalg.pinv(S @ S.T + (N - 1) * (1 / taper_coeff**2) * Cdd) - # ) - - A[i, :] = A[i, :] + np.sqrt(taper_coeff) * K @ (Dprime) - else: - K = A_centered @ S.T @ np.linalg.pinv(S @ S.T + (N - 1) * Cdd) - A = A + K @ Dprime - - for i in range(N): - fields[i][k_obs:] = A[:, i].reshape(nx, nx) - - k_start = k_obs - -# %% [markdown] -# ## Plot difference between prior and posterior of a single update - -# %% -k_obs = 5 -prior_mean_field = A_no_update[k_obs].mean(axis=1).reshape(nx, nx) -posterior_mean_field = matrix_from_fields(fields, k_obs).mean(axis=1).reshape(nx, nx) - -fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 6)) -p1 = axes[0].pcolormesh(posterior_mean_field - prior_mean_field) -utils.colorbar(p1) - -axes[1].plot(posterior_mean_field.ravel() - prior_mean_field.ravel()) - -fig.tight_layout() - - -# %% [markdown] -# ## Interactive plotting - - -# %% -def updated_vs_truth(k): - vmin = 0 - vmax = 100 - fig, axes = plt.subplots(nrows=1, ncols=2) - fig.suptitle(f"k = {k}") - - A_with_update = matrix_from_fields(fields, k) - - axes[0].set_title("With update") - axes[0].invert_yaxis() - p0 = axes[0].pcolormesh( - A_with_update.mean(axis=1).reshape(nx, nx).T, - cmap=plt.cm.viridis, - vmin=vmin, - vmax=vmax, - ) - - axes[1].set_title("Truth") - axes[1].invert_yaxis() - p1 = axes[1].pcolormesh(u_t[k].T, cmap=plt.cm.viridis, vmin=vmin, vmax=vmax) - - utils.colorbar(p1) - - fig.tight_layout() - - -interact( - updated_vs_truth, - k=widgets.IntSlider( - min=0, max=d.index.get_level_values("k").max(), step=1, value=0 - ), -) - -# %% diff --git a/notebooks/IES_ESMDA_Burgers.py b/notebooks/IES_ESMDA_Burgers.py deleted file mode 100644 index ffd2094..0000000 --- a/notebooks/IES_ESMDA_Burgers.py +++ /dev/null @@ -1,440 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.14.1 -# kernelspec: -# display_name: Python 3 (ipykernel) -# language: python -# name: python3 -# --- - -# %% [markdown] -# # Iterative Ensemble Smoother (IES) and Ensemble Smoother with Multiple Data Assimilation (ES-MDA) with the Burgers PDE -# -# Using IES and ES-MDA to estimate fluid viscosity. - -# %% -import matplotlib.pyplot as plt - -import numpy as np - -rng = np.random.default_rng(12345) - -import pandas as pd - -from p_tqdm import p_map -from ipywidgets import interact -import ipywidgets as widgets - -from iterative_ensemble_smoother import IterativeEnsembleSmoother - -# %% -# %load_ext autoreload -# %autoreload 2 -from dass import pde, utils, analysis, taper - -# %% [markdown] -# ## Define true parameter and run forward model, i.e., Burgers' equation. - -# %% -nx = 41 -ny = 41 - -# time steps -k_start = 0 -k_end = 500 - -nu_true = 0.1 -u = pde.burgers(nx=nx, ny=ny, nt=k_end, nu=nu_true) - - -# %% -def interactive_truth(k): - fig, ax = plt.subplots() - fig.suptitle("True response") - p = ax.pcolormesh(u[k]) - ax.set_title(f"k = {k}") - utils.colorbar(p) - fig.tight_layout() - - -interact( - interactive_truth, - k=widgets.IntSlider(min=k_start, max=k_end - 1, step=1, value=0), -) - -# %% [markdown] -# ## Define placement of sensors and generate synthetic observations based on the true field values - -# %% -# placement of sensors, i.e, where the observations are done -pad = 1 -coords = np.array([(x, y) for x in range(pad, nx - pad) for y in range(pad, nx - pad)]) -ncoords = coords.shape[0] -nmeas = 200 -coords_idx = np.random.choice(np.arange(ncoords), size=nmeas, replace=False) -obs_coordinates = [utils.Coordinate(xc, yc) for xc, yc in coords[coords_idx]] - -# At which times observations are taken -obs_times = np.linspace(5, k_end, 20, endpoint=False, dtype=int) - -d = utils.observations(obs_coordinates, obs_times, u, lambda value: abs(0.01 * value)) -# number of measurements -m = d.shape[0] -print("Number of observations: ", m) - -# Plot temperature field and show placement of sensors. -obs_coordinates = set(zip(d.index.get_level_values("x"), d.index.get_level_values("y"))) -x, y = zip(*obs_coordinates) - -fig, ax = plt.subplots() -p = ax.pcolormesh(u[0], cmap=plt.cm.viridis) -ax.set_title("True velocity field with sensor placement") -utils.colorbar(p) -ax.plot([i + 0.5 for i in x], [j + 0.5 for j in y], "s", color="white", markersize=5) - -# %% -# Number of realisations -N = 50 - -# %% -ss = np.random.SeedSequence(12345) -child_seeds = ss.spawn(N) -streams = [np.random.default_rng(s) for s in child_seeds] - -# %% -# List of matrices of size (nx, nx) containing priors. -# The reason for having a list is that `p_map` requires it. -# `p_map` runs stuff in parallel. -A_IES = rng.uniform(low=0.0, high=0.6, size=N) - -# %% -fig, ax = plt.subplots() -ax.hist(A_IES, bins=30) - - -# %% -def get_steplength(iteration_nr: int) -> float: - """ - This is an implementation of Eq. (49), which calculates a suitable step length for - the update step, from the book: - - Geir Evensen, Formulating the history matching problem with consistent error statistics, - Computational Geosciences (2021) 25:945 –970 - """ - max_steplength: float = 0.6 - min_steplength: float = 0.3 - dec_steplength: float = 2.5 - - steplength = min_steplength + (max_steplength - min_steplength) * pow( - 2, -(iteration_nr - 1) / (dec_steplength - 1) - ) - return steplength - - -# %% -number_of_ies_iterations = 3 - -# Coefficient matrix as defined in Eq. 16 and Eq. 17. -W = np.zeros(shape=(N, N)) - -# Use same prior for all methods. -A_IES_ert = A_IES.copy() -A_ESMDA = A_IES.copy() - -# %% -df_IES = pd.DataFrame({"Iteration": [0] * len(A_IES), "Value": A_IES}) -df_IES_ert = pd.DataFrame({"Iteration": [0] * len(A_IES_ert), "Value": A_IES_ert}) -df_ESMDA = pd.DataFrame({"Iteration": [0] * len(A_ESMDA), "Value": A_ESMDA}) - -# %% [markdown] -# ## Running ES for reference - -# %% -A_ES = A_IES.copy() - -fwd_runs = p_map( - pde.burgers, - [nx] * N, - [ny] * N, - [k_end] * N, - A_ES, - desc=f"Running forward model.", -) - -# Assume diagonal ensemble covariance matrix for the measurement perturbations. -Cdd = np.diag(d.sd.values**2) - -# 9.4 Ensemble representation for measurements -E = rng.multivariate_normal(mean=np.zeros(len(Cdd)), cov=Cdd, size=N).T -E = E - E.mean(axis=1, keepdims=True) -assert E.shape == (m, N) - -D = np.ones((m, N)) * d.value.values.reshape(-1, 1) + E - -Y = np.array( - [ - fwd_run[ - d.index.get_level_values("k").to_list(), - d.index.get_level_values("y").to_list(), - d.index.get_level_values("x").to_list(), - ] - for fwd_run in fwd_runs - ] -).T - -# We expect that a certain number of points gets deactivated, -# since not all cells are changing values due to the way the model works. -enough_ens_var_idx = Y.var(axis=1) > 1e-6 -print(f"{list(enough_ens_var_idx).count(False)} measurements will be deactivated.") -Y = Y[enough_ens_var_idx, :] -D = D[enough_ens_var_idx, :] -Cdd = Cdd[enough_ens_var_idx, :] -Cdd = Cdd[:, enough_ens_var_idx] - -X = analysis.ES(Y, D, Cdd) -A_ES = A_ES @ X - -# %% [markdown] -# ## Running `dass`'s implementation of IES - -# %% -for i in range(number_of_ies_iterations): - # Step length in Gauss Newton - gamma = get_steplength(i) - - fwd_runs = p_map( - pde.burgers, - [nx] * N, - [ny] * N, - [k_end] * N, - A_IES, - desc=f"Running forward model.", - ) - - # Assume diagonal ensemble covariance matrix for the measurement perturbations. - Cdd = np.diag(d.sd.values**2) - - # 9.4 Ensemble representation for measurements - E = rng.multivariate_normal(mean=np.zeros(len(Cdd)), cov=Cdd, size=N).T - E = E - E.mean(axis=1, keepdims=True) - assert E.shape == (m, N) - - D = np.ones((m, N)) * d.value.values.reshape(-1, 1) + E - - Y = np.array( - [ - fwd_run[ - d.index.get_level_values("k").to_list(), - d.index.get_level_values("y").to_list(), - d.index.get_level_values("x").to_list(), - ] - for fwd_run in fwd_runs - ] - ).T - - assert Y.shape == ( - m, - N, - ), "Measured responses must be a matrix with dimensions (number of observations x number of realisations)" - - # We expect that a certain number of points gets deactivated, - # since not all cells are changing values due to the way the model works. - enough_ens_var_idx = Y.var(axis=1) > 1e-6 - print(f"{list(enough_ens_var_idx).count(False)} measurements will be deactivated.") - Y = Y[enough_ens_var_idx, :] - D = D[enough_ens_var_idx, :] - Cdd = Cdd[enough_ens_var_idx, :] - Cdd = Cdd[:, enough_ens_var_idx] - - W = analysis.IES(Y, D, Cdd, W, gamma) - X_IES = np.identity(N) + W - A_IES = A_IES @ X_IES - - df_iter = pd.DataFrame({"Iteration": [i + 1] * len(A_IES), "Value": A_IES}) - df_IES = pd.concat([df_IES, df_iter]) - - -# %% -df_IES_pivot = df_IES.pivot(columns="Iteration", values="Value") -df_IES_pivot.columns.name = "" -df_IES_pivot.columns = ["Value_" + str(s) for s in df_IES_pivot.columns] - -# %% -fig, ax = plt.subplots(1, number_of_ies_iterations + 1) -fig.set_size_inches(8, 3) -df_IES_pivot.hist(ax=ax) -fig.tight_layout() - -# %% [markdown] -# ## Running ERT's implementation of IES - -# %% -for i in range(number_of_ies_iterations): - # Step length in Gauss Newton - gamma = get_steplength(i) - - fwd_runs = p_map( - pde.burgers, - [nx] * N, - [ny] * N, - [k_end] * N, - A_IES_ert.ravel(), - desc=f"Running forward model.", - ) - - Y = np.array( - [ - fwd_run[ - d.index.get_level_values("k").to_list(), - d.index.get_level_values("y").to_list(), - d.index.get_level_values("x").to_list(), - ] - for fwd_run in fwd_runs - ] - ).T - - assert Y.shape == ( - m, - N, - ), "Measured responses must be a matrix with dimensions (number of observations x number of realisations)" - - # We expect that a certain number of points gets deactivated, - # since not all cells are changing values due to the way the model works. - enough_ens_var_idx = Y.var(axis=1) > 1e-6 - print(f"{list(enough_ens_var_idx).count(False)} measurements will be deactivated.") - Y = Y[enough_ens_var_idx, :] - - A_IES_ert = IterativeEnsembleSmoother(N).update_step( - Y, - A_IES_ert.reshape(1, -1), - d.sd.values[enough_ens_var_idx], - d.value.values[enough_ens_var_idx], - step_length=gamma, - ) - - df_iter = pd.DataFrame( - {"Iteration": [i + 1] * len(A_IES_ert.ravel()), "Value": A_IES_ert.ravel()} - ) - df_IES_ert = pd.concat([df_IES_ert, df_iter]) - -# %% -df_IES_ert_pivot = df_IES_ert.pivot(columns="Iteration", values="Value") -df_IES_ert_pivot.columns.name = "" -df_IES_ert_pivot.columns = ["Value_" + str(s) for s in df_IES_ert_pivot.columns] - -# %% -fig, ax = plt.subplots(1, number_of_ies_iterations + 1) -fig.set_size_inches(8, 3) -df_IES_pivot.hist(ax=ax) -fig.tight_layout() - -# %% [markdown] -# ## Ensemble Smoother with Multiple Data Assimilation (ES-MDA) - -# %% -weights = [4, 2, 1] - -# %% -# Same function is used in ERT. -from typing import List - - -def normalize_weights(weights: List[float]) -> List[float]: - """Scale weights such that their reciprocals sum to 1.0, - i.e., sum(1.0 / x for x in weights) == 1.0. - See for example Equation 38 of evensen2018 - Analysis of iterative - ensemble smoothers for solving inverse problems. - """ - if not weights: - return [] - weights = [weight for weight in weights if abs(weight) != 0.0] - - length = sum(1.0 / x for x in weights) - return [x * length for x in weights] - - -# %% -normalized_weights = normalize_weights(weights) - -# %% -for i, w in enumerate(normalized_weights): - fwd_runs = p_map( - pde.burgers, - [nx] * N, - [ny] * N, - [k_end] * N, - A_ESMDA, - desc=f"Running forward model.", - ) - - # Assume diagonal ensemble covariance matrix for the measurement perturbations. - # Inflating measurement errors by a factor sqrt(normalized_weights) as shown - # in for example evensen2018 - Analysis of iterative ensemble smoothers for solving inverse problems. - Cdd = np.diag((d.sd.values * np.sqrt(w)) ** 2) - - # 9.4 Ensemble representation for measurements - E = rng.multivariate_normal(mean=np.zeros(len(Cdd)), cov=Cdd, size=N).T - E = E - E.mean(axis=1, keepdims=True) - assert E.shape == (m, N) - - D = np.ones((m, N)) * d.value.values.reshape(-1, 1) + E - - Y = np.array( - [ - fwd_run[ - d.index.get_level_values("k").to_list(), - d.index.get_level_values("y").to_list(), - d.index.get_level_values("x").to_list(), - ] - for fwd_run in fwd_runs - ] - ).T - - assert Y.shape == ( - m, - N, - ), "Measured responses must be a matrix with dimensions (number of observations x number of realisations)" - - # We expect that a certain number of points gets deactivated, - # since not all cells are changing values due to the way the model works. - enough_ens_var_idx = Y.var(axis=1) > 1e-6 - print(f"{list(enough_ens_var_idx).count(False)} measurements will be deactivated.") - Y = Y[enough_ens_var_idx, :] - D = D[enough_ens_var_idx, :] - Cdd = Cdd[enough_ens_var_idx, :] - Cdd = Cdd[:, enough_ens_var_idx] - - X_ESMDA = analysis.ES(Y, D, Cdd) - A_ESMDA = A_ESMDA @ X_ESMDA - - df_iter = pd.DataFrame( - {"Iteration": [i + 1] * len(A_ESMDA.ravel()), "Value": A_ESMDA.ravel()} - ) - df_ESMDA = pd.concat([df_ESMDA, df_iter]) - -# %% -df_ESMDA_pivot = df_ESMDA.pivot(columns="Iteration", values="Value") -df_ESMDA_pivot.columns.name = "" -df_ESMDA_pivot.columns = ["Value_" + str(s) for s in df_ESMDA_pivot.columns] - -# %% -fig, ax = plt.subplots(1, len(normalized_weights) + 1) -fig.suptitle("ES-MDA") -fig.set_size_inches(8, 3) -df_ESMDA_pivot.hist(ax=ax) -fig.tight_layout() - -# %% -fig, ax = plt.subplots(1, number_of_ies_iterations + 1) -fig.suptitle("IES") -fig.set_size_inches(8, 3) -df_IES_pivot.hist(ax=ax) -fig.tight_layout() - -# %% diff --git a/notebooks/Poly.py b/notebooks/Poly.py deleted file mode 100644 index 4ad10ff..0000000 --- a/notebooks/Poly.py +++ /dev/null @@ -1,158 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.14.4 -# kernelspec: -# display_name: Python 3 (ipykernel) -# language: python -# name: python3 -# --- - -# %% [markdown] -# # Poly case - -# %% -import numpy as np -import pandas as pd - -np.set_printoptions(suppress=True) -rng = np.random.default_rng(12345) - -import matplotlib.pyplot as plt - -plt.rcParams["figure.figsize"] = (6, 6) -plt.rcParams.update({"font.size": 10}) -from ipywidgets import interact -import ipywidgets as widgets - -from p_tqdm import p_map - -from scipy.ndimage import gaussian_filter - -import iterative_ensemble_smoother as ies - -# %% -# %load_ext autoreload -# %autoreload 2 -from dass import pde, utils, analysis, taper - -# %% -N = 200 - - -# %% -def poly(a, b, c, x): - return a * x**2 + b * x + c - - -# %% -a_t = 0.5 -b_t = 1.0 -c_t = 3.0 - -x_observations = [0, 2, 4, 6, 8] -observations = [ - ( - poly(a_t, b_t, c_t, x) + rng.normal(loc=0, scale=0.2 * poly(a_t, b_t, c_t, x)), - 0.2 * poly(a_t, b_t, c_t, x), - x, - ) - for x in x_observations -] - -d = pd.DataFrame(observations, columns=["value", "sd", "x"]) - -d = d.set_index("x") - -m = d.shape[0] - -# %% -fig, ax = plt.subplots() -x_plot = np.linspace(0, 10, 50) -ax.set_title("Truth and noisy observations") -ax.set_xlabel("Time step") -ax.set_ylabel("Response") -ax.plot(x_plot, poly(a_t, b_t, c_t, x_plot)) -ax.plot(d.index.get_level_values("x"), d.value.values, "o") -ax.grid() - -# %% -# Assume diagonal ensemble covariance matrix for the measurement perturbations. -Cdd = np.diag(d.sd.values**2) - -E = rng.multivariate_normal(mean=np.zeros(len(Cdd)), cov=Cdd, size=N).T -E = E - E.mean(axis=1, keepdims=True) -assert E.shape == (m, N) - -D = np.ones((m, N)) * d.value.values.reshape(-1, 1) + E - -# %% -coeff_a = rng.normal(0, 1, size=N) -coeff_b = rng.normal(0, 1, size=N) -coeff_c = rng.normal(0, 1, size=N) - -# %% -A = np.concatenate( - (coeff_a.reshape(-1, 1), coeff_b.reshape(-1, 1), coeff_c.reshape(-1, 1)), axis=1 -).T - -# %% -fwd_runs = p_map( - poly, - coeff_a, - coeff_b, - coeff_c, - [np.arange(max(x_observations) + 1)] * N, - desc=f"Running forward model.", -) - -# %% -Y = np.array( - [fwd_run[d.index.get_level_values("x").to_list()] for fwd_run in fwd_runs] -).T - -assert Y.shape == ( - m, - N, -), "Measured responses must be a matrix with dimensions (number of observations x number of realisations)" - -# %% -X = analysis.ES(Y, D, Cdd) -A_ES = A @ X - -# %% -fig, ax = plt.subplots(nrows=3, ncols=2) -fig.set_size_inches(10, 10) - -ax[0, 0].set_title("a - prior") -ax[0, 0].hist(A[0, :]) - -ax[0, 1].set_title("a - posterior") -ax[0, 1].hist(A_ES[0, :]) -ax[0, 1].axvline(a_t, color="k", linestyle="--", label="truth") -ax[0, 1].legend() - -ax[1, 0].set_title("b - prior") -ax[1, 0].hist(A[1, :]) - -ax[1, 1].set_title("b - posterior") -ax[1, 1].hist(A_ES[1, :]) -ax[1, 1].axvline(b_t, color="k", linestyle="--", label="truth") -ax[1, 1].legend() - -ax[2, 0].set_title("c - prior") -ax[2, 0].hist(A[2, :]) - -ax[2, 1].set_title("c - posterior") -ax[2, 1].hist(A_ES[2, :]) -ax[2, 1].axvline(c_t, color="k", linestyle="--", label="truth") -ax[2, 1].legend() - -fig.tight_layout() - -# %% diff --git a/setup.py b/setup.py index dd6f45d..0bfefbb 100644 --- a/setup.py +++ b/setup.py @@ -13,6 +13,5 @@ "jupytext", "pandas>=1.4.2,<2.0.0", "scipy>=1.9.0,<2.0.0", - "iterative_ensemble_smoother>=0.0.5,<1.0.0", ], )