iterative_ensemble_smoother.ESMDA#

class iterative_ensemble_smoother.ESMDA(covariance: ndarray[Any, dtype[float64]], observations: ndarray[Any, dtype[float64]], alpha: int | ndarray[Any, dtype[float64]] = 5, seed: Generator | int | None = None, inversion: str = 'exact')[source]#

Bases: object

Implement an Ensemble Smoother with Multiple Data Assimilation (ES-MDA).

The implementation follows Emerick and Reynolds [2013].

Parameters:

  • covariance (np.ndarray) – Either a 1D array of diagonal covariances, or a 2D covariance matrix. The shape is either (num_observations,) or (num_observations, num_observations). This is C_D in Emerick (2013), and represents observation or measurement errors. We observe d from the real world, y from the model g(x), and assume that d = y + e, where the error e is multivariate normal with covariance given by covariance.

  • observations (np.ndarray) – 1D array of shape (num_observations,) representing real-world observations. This is d_obs in Emerick (2013).

  • alpha (int or 1D np.ndarray, optional) – Multiplicative factor for the covariance. If an integer alpha is given, an array with length alpha and elements alpha is constructed. If an 1D array is given, it is normalized so sum_i 1/alpha_i = 1 and used. The default is 5, which corresponds to np.array([5, 5, 5, 5, 5]).

  • seed (integer or numpy.random._generator.Generator, optional) – A seed or numpy.random._generator.Generator used for random number generation. The argument is passed to numpy.random.default_rng(). The default is None.

  • inversion (str, optional) – Which inversion method to use. The default is “exact”. See the dictionary ESMDA._inversion_methods for more information.

Examples

>>> covariance = np.diag([1, 1, 1])
>>> observations = np.array([1, 2, 3])
>>> esmda = ESMDA(covariance, observations)

Methods

assimilate(X, Y, *[, overwrite, truncation])

Assimilate data and return an updated ensemble X_posterior.

compute_transition_matrix(Y, *, alpha[, ...])

Return a matrix K such that X_posterior = X_prior + X_prior @ K.

perturb_observations(*, size, alpha)

Create a matrix D with perturbed observations.

num_assimilations
__init__(covariance: ndarray[Any, dtype[float64]], observations: ndarray[Any, dtype[float64]], alpha: int | ndarray[Any, dtype[float64]] = 5, seed: Generator | int | None = None, inversion: str = 'exact') → None[source]#

Initialize the instance.

Methods definition

assimilate(X: ndarray[Any, dtype[float64]], Y: ndarray[Any, dtype[float64]], *, overwrite: bool = False, truncation: float = 1.0) → ndarray[Any, dtype[float64]][source]#

Assimilate data and return an updated ensemble X_posterior.

Parameters:

  • X (np.ndarray) – A 2D array of shape (num_parameters, ensemble_size). Each row corresponds to a parameter in the model, and each column corresponds to an ensemble member (realization).

  • Y (np.ndarray) – 2D array of shape (num_parameters, ensemble_size), containing responses when evaluating the model at X. In other words, Y = g(X), where g is the forward model.

  • overwrite (bool) – If True, then arguments X and Y may be overwritten and mutated. If False, then the method will not mutate inputs in any way. Setting this to True might save memory.

  • truncation (float) – How large a fraction of the singular values to keep in the inversion routine. Must be a float in the range (0, 1]. A lower number means a more approximate answer and a slightly faster computation.

Returns:

X_posterior – 2D array of shape (num_inputs, num_ensemble_members).

Return type:

np.ndarray

compute_transition_matrix(Y: ndarray[Any, dtype[float64]], *, alpha: float, truncation: float = 1.0) → ndarray[Any, dtype[float64]][source]#

Return a matrix K such that X_posterior = X_prior + X_prior @ K.

The purpose of this method is to facilitate row-by-row, or batch-by-batch, updates of X. This is useful if X is too large to fit in memory.

Parameters:

  • Y (np.ndarray) – 2D array of shape (num_parameters, ensemble_size), containing responses when evaluating the model at X. In other words, Y = g(X), where g is the forward model.

  • alpha (float) – The covariance inflation factor. The sequence of alphas should obey the equation sum_i (1/alpha_i) = 1. However, this is NOT enforced in this method call. The user/caller is responsible for this.

  • truncation (float) – How large a fraction of the singular values to keep in the inversion routine. Must be a float in the range (0, 1]. A lower number means a more approximate answer and a slightly faster computation.

Returns:

K – A matrix K such that X_posterior = X_prior + X_prior @ K. It has shape (num_ensemble_members, num_ensemble_members).

Return type:

np.ndarray

perturb_observations(*, size: Tuple[int, int], alpha: float) → ndarray[Any, dtype[float64]][source]#

Create a matrix D with perturbed observations.

In the Emerick (2013) paper, the matrix D is defined in section 6. See section 2(b) of the ES-MDA algorithm in the paper.

Parameters:

  • size (Tuple[int, int]) – The size, a tuple with (num_observations, ensemble_size).

  • alpha (float) – The covariance inflation factor. The sequence of alphas should obey the equation sum_i (1/alpha_i) = 1. However, this is NOT enforced in this method call. The user/caller is responsible for this.

Returns:

D – Each column consists of perturbed observations, scaled by alpha.

Return type:

np.ndarray