Four Dimensional Variation (4DVar) DA

This module implements several variants of the four-dimensional variational (4D-Var) data assimilation scheme. 4D-Var refers to simulation/prediction in three-dimensional spatial coordinates + one temporal dimension. The name 4D-Var is popular in the atmospheric weather prediction (NWP) literature. This algorithm is equivalent to standard regularized least-squares inversion. Multiple observations (over a time window) are used to update the model parameter/initial condition. When the time window collapses to one observation instance, the 4D-Var reduces to 3D-Var scheme

class VanillaFourDVarConfigs(*, debug=False, verbose=False, output_dir='./_PYOED_RESULTS_', name='4D-Var: Four Dimensional Variational Data Assimilation', model=None, prior=None, observation_error_model=None, observation_operator=None, observations=None, optimizer=<class 'pyoed.optimization.scipy_optimization.ScipyOptimizer'>, optimizer_configs=<class 'pyoed.optimization.scipy_optimization.ScipyOptimizerConfigs'>, window=None)

Bases: VariationalSmootherConfigs

Configurations class for the VanillaDVar abstract base class.

Parameters:

• verbose (bool) – a boolean flag to control verbosity of the object.

• debug (bool) – a boolean flag that enables adding extra functionlity in a debug mode

• output_dir (str | Path) – the base directory where the output files will be saved.

• name (str | None) – name of the DA (inverse problem solver) approach/method.

• model (None | SimulationModel) – the simulation model.

• prior (None | ErrorModel) – Background/Prior model (e.g., GaussianErrorModel)

• observation_operator (None | ObservationOperator) – operator to map model state to observation space

• observation_error_model (None | ErrorModel) – Observation error model (e.g., GaussianErrorModel)

• observations (None | Any) – Observational data (the data type is very much dependent of the DA method)

• optimizer (None | Optimizer) –

  the optimization routine (optimizer) to be registered and later used for solving the OED problem. This can be one of the following:

  • None: In this case, no optimizer is registered, and the solve() won’t be functional until an optimization routine is registered.

  • An optimizer instance (object that inherits :py:class`Optimizer`). In this case, the optimizer is registered as is and is updated with the passed configurations if available.

  • The class (subclass of Optimizer) to be used to instantiate the optimizer.

• optimizer_configs (None | OptimizerConfigs) –

  the configurations of the optimization routine. This can be one of the following:

  • None, in this case configurations are discarded, and whatever default configurations of the selected/passed optimizer are employed.

  • A dict holding full/partial configurations of the selected optimizer. These are either used to instantiate or update the optimizer configurations based on the type of the passed optimizer.

  • A class providing implementations of the configurations (this must be a subclass of OptimizerConfigs.

  • An instance of a subclass of OptimizerConfigs which is to set/udpate optimizer configurations.

  Note

  Not all DA (inverse problem) objects are optimization-based. For example, particle-based (EnKF, PF, etc.) employ a sample to estimate the flow of the distribution through the model dynamics (prior -> posterior). Thus, the optimizer (and configs) in this case (the default) are set to None. For optimization-based methods such as 3DVar, 4DVar, etc., an optimizer must be registered for the inverse problem to be solved.

• window (None | Iterable) – the assimilation window (t0, tf)

__init__(*, debug=False, verbose=False, output_dir='./_PYOED_RESULTS_', name='4D-Var: Four Dimensional Variational Data Assimilation', model=None, prior=None, observation_error_model=None, observation_operator=None, observations=None, optimizer=<class 'pyoed.optimization.scipy_optimization.ScipyOptimizer'>, optimizer_configs=<class 'pyoed.optimization.scipy_optimization.ScipyOptimizerConfigs'>, window=None)

class VanillaFourDVar(configs=None)

Bases: VariationalSmoother

A class implementing the vanilla 4D-Var DA scheme to invert for model initial state. This version is very elementary and is developed as guidelines to create more advanced versions, e.g., parameter-retrieval, or extended (parameter-state) retrieval for inifinite dimensional formulations; see remrks below for assumptions about the dynamical model.

This scheme assumes Gaussian observational noise. The model unknown state is regularized (to avoid overfitting to noise) by assuming a prior (usually a Gaussian), or by specifying a regularization matrix and following an \(\ell_2\) regularization approach.

__init__(configs=None)

register_optimizer(optimizer='ScipyOptimizer', optimizer_configs=None)

Register (and return) an optimization routine, and make sure the objective function in the optimization routine is set to the objctive function of this objective objective_function_value(). The objective function of the optimizer (and its derivatives) are set to the objective function and derivatives implemented here unless the objective function is passed explicitly in the configurations. This method sets a reference to the optimizer and update the underlying configurations accordingly.

Note

If the optimizer passed is an instance of (derived from) Optimizer, and valid configurations optimizer_configs is passed, the optimizer is updated with the passed configurations by calling optimizer.update_configurations(optimizer_configs).

Warning

If the optimizer passed is an instance of (derived from) Optimizer, the responsibility is on the developer/user to validate the contents of the optimizer.

Parameters:

• optimizer (str | Optimizer | Type[Optimizer]) –

  the optimization routine (optimizer) to register. This can be one of the following:

  • An optimizer instance (object that inherits :py:class`Optimizer`). In this case, the optimizer is registered as is and is updated with the passed configurations if available.

  • The class (subclass of Optimizer) to be used to instantiate the optimizer. To see available built-in optimization routines, see/call show_supported_optimizers().

  • The name of the optimizer (str). This has to match the name of one of the available optimization routine. To see available built-in optimization routines, see/call show_supported_optimizers().

• optimizer_configs (None | dict | OptimizerConfigs | Type[OptimizerConfigs]) –

  the configurations of the optimization routine. This can be one of the following:

  • None, in this case configurations are discarded, and whatever default configurations of the selected/passed optimizer are employed.

  • A dict holding full/partial configurations of the selected optimizer. These are either used to instantiate or update the optimizer configurations based on the type of the passed optimizer.

  • A class providing implementations of the configurations (this must be a subclass of OptimizerConfigs.

  • An instance of a subclass of OptimizerConfigs which is to set/udpate optimizer configurations.

Returns:

the registered optimizer.

Raises:

TypeError – if the type of passed optimizer and/or configurations are/is not supported

Return type:

Type[Optimizer]

validate_configurations(configs, raise_for_invalid=True)

Validate passed configurations.

Parameters:

configs (dict | VanillaFourDVarConfigs) – configurations to validate. If a VanillaFourDVarConfigs object is passed, validation is performed on the entire set of configurations. However, if a dictionary is passed, validation is performed only on the configurations corresponding to the keys in the dictionary.

update_configurations(**kwargs)

Take any set of keyword arguments, and lookup each in the configurations, and update as nessesary/possible/valid

Raises:

TypeError – if any of the passed keys in kwargs is invalid/unrecognized

Remarks:

• Generally, we don’t want actual implementations in abstract classes, however, this one is provided as good guidance. Derived classes can rewrite it and/or provide additional updates.

register_model(model=None)

Register (and return) the simulation model to be registered. This calls InverseProblem.register_model and adds extra assertions/functionality specific for filters.

Raises:

TypeError – if the type of passed model is not supported

register_prior(prior=None, initiate_posterior=True)

Update prior/regularization information. Since, changing the prior implies a different posterior, the argument initiate_posterior is added with defaul value True to make sure the posterio is overwritten/reset

Raises:

TypeError is raised if:

• the passed prior type is not supported

• a valid simulation model has not been registered yet

• the prior is not compatible with the registered simulation model

Note

The prior here is asserted on the model state or model parmeter based on the invert_for flag selected in the configurations dictionary

register_observation_operator(observation_operator)

The observation operator is a function/map that takes as input a model state, and returns an observation instance.

Parameters:

observation_operator (None | ObservationOperator) – the obsevation operator

Raises:

• TypeError – if observation_operator is not derived from ObservationOperator

• TypeError – if the passed operator doesn’t provide a function to apply the tangent linear of the observation operator ‘Jacobian_T_matvec’

• TypeError – if the model state does not match the observation operator domain size

Note

We are assuming the observation operator is time-independent, however, it provides a method apply and/or __call__ which converts a model state into an equivalent observation

solve(init_guess=None, skip_map_estimate=False, update_posterior=False)

Solve the inverse problem, i.e., find the model state given the registered observation and prior information

Parameters:

• init_guess – initial guess of the model state

• skip_map_estimate (bool) – use the prior mean as a map estimte (ignore observations) and upate posterior. This useful for OED when the model is linear

• update_posterior (bool) – if True, the posterior mean and covariance operator are updated (given/around the solution of the inverse problem)

Returns:

the analysis state (MAP estimate of the posterior)

Raises:

• TypeError if any of the 4D-Var elements (model, observation operator, prior, solver, or data) is missing

• ValueError if the underlying optimization routine is not supported yet

• AttributeError is raised if the passed model doesn not provide a method ‘Jacobian_T_matvec’ to apply the tangent linear model (TLM/Jacobian) transpose of the right-hand-side of the discretized PDE to a model state.

show_registered_elements(display=True)

Compose and (optionally) print out a message containing the elements of the inverse problem and show what’s registered and what’s not

Parameters:

display – if True print out the composed message about registered/unregistered elements

Returns:

the composed message

Return type:

None

check_registered_elements()

Check if all inversion elements are registered or not

objective_function_value(state, data_misfit_only=False, return_traject=False)

Evaluate the value of the 4D-Var objective function (data misfit + prior/regularization) given the passed state as an estimate of the true model state/parameter

Parameters:

• state – state giving an estimate of the model initial state/parameter

• return_traject (bool) – if True the checkpoints and checkpointed states (forward traject only at the observation time points) are returned along with the value of the objective function

Returns:

value of the 4D-Var objective function. The checkpoints and checkpointed states are returned only if return_traject is True

Raises:

• TypeError is raised if state is not a valid model state

• ValueError is raised if observations are not yet registered

• AssertionError is raised if less than two valid observations are registered or observation time(s) fall outside the assimilation window

• KeyError or AttributeError (as raised by the model) when the time integration step size of the simulation model is inaccessible; this is used only if return_traject is True and dt is None

objective_function_gradient(state, data_misfit_only=False, checkpoints=None, checkpointed_states=None)

Evaluate the gradient of the 4D-Var objective function given the passed state as an estimate of the true model state/parameter

Parameters:

• state – state giving an estimate of the model initial state/parameter

• checkpoints – used if checkpointed states are passed, otherwise ignored

• checkpointed_states – trajectory of states evaluated at the passed checkpoints, starting at the passed state. This is useful to reduce redundancy in the optimization procedure

Returns:

the gradient of the 4D-Var cost function

Raises:

ValueError if checkpoints and/or checkpointed_states are passed but invalid or incompatible

apply_forward_operator(state, save_states=False, checkpoints=None, scale_by_noise=True)

Apply \(F\), the forward operator to the passed state. The forward operator \(F\) here, refers to the simulation model followed by the observation operator. The result is a data point (an observation), or a dictionary of observations indexed by the time (for time-dependent models)

For time-dependent simulations with multiple observation points (e.g., in 4D-Var settings), the observations are evaluated at the passed checkpoints (if they correspond a registered observation time) If there is an observation at the initial time, the forward operator at the initial time corresponds to applying the observation operator without applying model simulation/dynamics.

Parameters:

• state – data structure holding the model state at the initial time of the registered assimilation window

• checkpoints – times at which to store the computed observations, must be sorted and lie within the registered window. If None (default), use the registered observation times

• save_states – if True states will be checkpointed and saved as well as observations

• scale_by_noise (bool) – if True, the observations are scaled by the observation error precision

Returns:

• checked_observations a dictionary holding times as keys, and observations as values. The observations result from forward integration of the simulation model, and applying observation operator (with/without) scaling by observation error precisions

• checked_states: a dictionary holding model states evaluated at the checkpoints; set to None if save_states is False, otherwise

apply_forward_operator_adjoint(obs, obs_time, eval_at=None, checkpointed_states=None, checkpointed_states_times=None, save_all=False, scale_by_noise=True)

Apply F^*, the adjoint of the forward operator, to the passed observation at the given observation time obs_time. Here, F is a composition of the simulation model and the observation operator (e.g., simulation followed by restriction) applied at all observation times. Thus, the adjoint is applied recursively backward in time. If the passed obs_time coincides with the initial time of the assimilation window, the adjoint corresponds to applying the adjoint of the observation operator only, and the model dynamics in this case correspond to applying an identiy operator.

Parameters:

• obs – the vector in observation space to apply the adjoint to

• obs_time – time at which obs it given, if no observation time is given, it is set to the last registered observation time

• save_all (bool) – (default False) if True, adjoint evaluated at all time points is returned as a dictionary indexed by observation times, otherwise, only resulting adjoint at the initial time of the assimilation window.

• scale_by_noise (bool) – if True, the observations are scaled by the observation error precision

• eval_at – state to linearize the forward operator around, used if checkpointed states are not given to checkpoint forward trajectory

• checkpointed_states – trajectory of states evaluated at all checkpoints (t0 + observation times),

• checkpointed_states_times – times at which checkpointed_states are given This is useful to reduce redundancy in the optimization procedure. Both checkpointed_states and checkpointed_states_times are ignored for linear problems as they are not needed for adjoint evaluation

Returns:

adjoint, checked_adj where:

• adjoint: the adjoint evaluated at the initial time of the assimilation window

• checked_adj a dictionary (indexed by times) of checkpointed states

full_Hessian(eval_at=None, data_misfit_only=False)

Construct the Hessian by evaluating the tangent linear model

full_Hessian_inv(eval_at=None, data_misfit_only=False)

Construct the inverse of the Hessian by evaluating the tangent linear model

posterior_covariance(eval_at=None, data_misfit_only=False)

Construct the inverse of the Hessian by evaluating the tangent linear model

Hessian_matvec(state, eval_at=None, data_misfit_only=False)

Return the product of the Hessian of 4D-Var objective with a vector state.

Parameters:

• state – initial state/guess of the model state/parameter

• eval_at – state at which the Hessian is evaluated (passed to the model adjoint)

• data_misfit_only (bool) – if True the prior term of the Hessian is discarded, and only the data-misfit term is evaluated

Returns:

product of the inverse of the posterior covariance (of the linearized problem) with a state vector

Notes:

This method fixes and replace old code now moved to _incorrect_Hessian_matvec

Hessian_inv_matvec(state, eval_at=None, data_misfit_only=False, precond=None, solver='scipy-minres', solver_options={'callback': None, 'maxiter': None})

Return the product of the Hessian inverse of 4D-Var objective with a vector ‘state’

Parameters:

• state – initial state/guess of the model state/parameter

• eval_at – state at which the Hessian is evaluated (passed to the model adjoint)

• data_misfit_only (bool) – if True the prior term of the Hessian is discarded, and only the data-misfit term is evaluated

Returns:

product of the posterior covariance (of the linearized problem) with a state vector

Remarks:

this function is very important and resembles the most epensive part of linear/linearized Bayesian algorithms. Preconditioning, and/or reduced-order approxmiations are very important for performance here.

property supported_optimizers