text:
\section{Introduction}
\label{sec:Introduction}
  Recently, interest has increased in developing scalable data assimilation (DA) and uncertainty quantification methodologies for solving large-scale inverse problems.
  An inverse problem refers to the retrieval of a quantity of interest (QoI) associated with or stemming
  from a physical phenomenon underlying partial noisy experimental or observational data of that physical
  system~\cite{aster2018parameter,vogel2002computational,stuart2010inverse}.
  The QoI could be, for example, the model state, initial condition, or other physics quantity.
  Inverse problems are prominent in a wide spectrum of applications including
  power grids and atmospheric numerical weather prediction~\cite{ghil1991data,smith2013uncertainty}.
  In these problems, the prediction of the physical phenomena is often formulated as an initial value problem,
  while the initial condition of the simulator is corrected by fusing all available information.
  %
  Algorithmic approaches for solving inverse problems seek either a single-point estimate of the target QoI
  or a full probabilistic description of the knowledge about the QoI given all available information.
  The underlying principle of these methods is that information collected from observational systems
  is fused into computational models, along with associated uncertainties, to produce an accurate estimate
  of the ground truth of the physical phenomena of interest.
  %
  In the former approach seeking a single QoI estimate, the solution of an inverse problem is obtained
  by solving an optimization problem with an objective to minimize the mismatch between observational data
  and model simulations, possibly regularized by prior knowledge and uncertainty models.
  The latter approach, commonly known as Bayesian inversion,  seeks to characterize the
  probability distribution of the QoI through the posterior formulated by applying Bayes' rule,
  that is, the probability distribution of the QoI conditioned by all available information.

  DA methods~\cite{bannister2017review,daley1993atmospheric,navon2009data,smith2013uncertainty,ghil1991data,attia2015hmcfilter,attia2015hmcsampling,attia2015hmcsmoother,attia2017reduced,attia2019dates}
  aim to solve large- to extreme-scale inverse problems. They work by fusing
  information obtained from multiple sources, such as the dynamical model, prior knowledge,
  noisy and incomplete measurements, and error models, in order to better estimate the state and parameters
  of the physical system.
  This estimate  improves the predictability of the simulation systems developed to make future predictions
  about the physical phenomena of interest.

  The quality of DA systems, and hence the accuracy of their predictions, is heavily influenced by the extent to
  which the mathematical assumptions reflect reality and depends on the quality of the collected measurements.
  Optimal data acquisition is the problem of determining the optimal observational strategy,
  for example, from a large set of candidate observational schemes.
  This problem is widely formulated as an optimal experimental design (OED)
  problem~\cite{fedorov2000design,pukelsheim2006optimal}, where the design parameterizes and thus determines
  the observational configuration.
  In an OED problem, a design is defined to characterize a candidate configuration or a control, and the
  quality of the design is quantified by using a utility function.
  The optimal design is then defined as the one that maximizes this utility function or, equivalently,
  minimizes some OED criterion~\cite{attia2022optimal}.
  Since the aim of Bayesian inference is to estimate the QoI posterior, Bayesian OED seeks an observational
  configuration that, when combined with the underlying dynamics, would maximize information gain from the
  data or minimize the posterior uncertainty.
  Thus, an optimal design is found by solving an optimization problem with the objective to maximize a
  utility function that quantifies the quality of the design and its influence on the solution
  of the inverse problem.
  OED for inverse problems has experienced a recent surge in interest by the scientific computing community;
  see, for example,~\cite{alexanderian2021optimal} and references therein.

  Numerical testing and experiments are critical for developing efficient
  OED formulations and algorithms. This process is elementary for successful scientific research in general.
  Although statisticians have developed a plethora of mathematical formulations and algorithmic approaches
  for general-purpose OED algorithms, most of the available and publicly accessible OED software tools are
  limited to idealized formulations and specific applications such as finding optimal collocation points for
  regression problems or designing clinical experiments.
  In addition, they are written in the \textrm{R} programming language or MATLAB, thus limiting code reutilization
  and accessibility by a wider audience;
  see, for example,~\cite{foracchia2004poped,wheeler2019package,rasch2011optimal,tian2021autooed,rasch2009software}.

  Unfortunately, these tools do not align well with the interests of the increasingly large computational
  science community for developing new OED formulations and algorithmic approaches for inverse problems,
  DA, and model-constrained  OED~\cite{attia2018goal,alexanderian2016fast,huan2013simulation,bui2013computational,flath2011fast,haber2008numerical,haber2009numerical,pukelsheim2006optimal,pronzato2013design}.
  As a first step in alleviating this limitation, we present and describe {PyOED}\xspace,
  a highly extensible open-source Python package written mainly to enable computational scientists to formulate and rapidly test OED---as well as DA---formulations and algorithmic approaches.

  {PyOED}\xspace is unique in several ways.
  %
  First, to the best of our knowledge, it is the first open-source package for scientific computing that allows implementing and testing the individual components of DA as well as OED systems in a unified and streamlined environment.
  %
  Second, it is written in Python, which is arguably the most popular and adopted programming language for recent algorithmic developments in the computational science disciplines.
  It has a huge user-support community, and the learning curve is relatively smoother than that of other lower-level programming languages such as C/C++.
  %
  Third, {PyOED}\xspace is designed in an object-oriented programming (OOP) fashion, which enables practitioners to reconfigure and reuse the individual building blocks.
  Moreover, it is easy to combine {PyOED}\xspace with other user-defined routines, such as numerical integration of simulation models and optimization routines, which makes {PyOED}\xspace highly extensible and adaptable to a wide range of applications.
  %
  Fourth, {PyOED}\xspace is not limited to a specific inverse problem formulation. Thus new DA and OED methods can be implemented and even interface with other inversion packages, for example, hIPPYlib~\cite{villa2018hippylib}.
  %
  Fifth, {PyOED}\xspace leverages best practices in software development, including detailed documentation with hands-on examples and robust unit-testing techniques.

  The rest of this paper is organized as follows.
  Section~\ref{sec:Background} provides the mathematical formalism of inverse problems, DA, and OED.
  Section~\ref{sec:PyOED} describes the structure and the philosophy of the {PyOED}\xspace package.
  In Section~\ref{sec:Test_cases} we provide a list of numerical experiments to demonstrate
  the general workflow and usage of {PyOED}\xspace.
  Concluding remarks are given in Section~\ref{sec:Conclusions}.

\section{Mathematical Background}
\label{sec:Background}
  In this section we provide a brief overview of the mathematical background of the {PyOED}\xspace core,
  which is important for approaching the OED problem for Bayesian inversion.
  In this presentation we focus on sensor placement for Bayesian inversion as a modal OED formulation.
  We start by discussing the forward problem in~\ref{subsec:forward_problem}; then we introduce the
  inverse problem in~\ref{subsec:inverse_problem} and
  the OED formalism in~\ref{subsec:OED}.

  \subsection{The forward problem}
  \label{subsec:forward_problem}

    The forward problem maps the model parameters (e.g., the initial condition) onto the observation space.
    Consider the forward problem described by

    \begin{subequations}\label{eqn:forward_problem}
    \begin{equation}\label{eqn:forward_problem_base}
        \obs  = \Fcont(\param) + \obsnoise \,,
    \end{equation}
    %
    where $\param$ is the model parameter of interest,
    $\obs \in \Rnum^{\Nobs}$ is the observation, and
    $\obsnoise \in \Rnum^{\Nobs}$ is a noise term that accounts for the inaccuracy of the observational system.
    %
    The forward operator $\Fcont$ is occasionally referred to as the ``parameter-to-observable'' map
    and generally represents a composition of a simulation/solution model $\Sol$ and an observation
    operator $\mat{O}$.
    The simulation model $\Sol$ describes the evolution of the physical phenomena, for example,
    space-time advection and diffusion of a contaminant simulated over a predefined model grid.
    %
    The observation operator $\mat{O}$ projects the simulated state onto the observational grid,
    for example, by interpolation or restriction to the observation grid.
    Thus, the forward problem~\eqref{eqn:forward_problem_base} can be rewritten as
    %
    \begin{equation}\label{eqn:forward_problem_decomposed}
        \obs  = \mat{O} \circ \Sol(\param) + \obsnoise \,,
    \end{equation}
    %
    \end{subequations}
    %
    where $\circ$ is the composition operator, that is, $\mat{O} \circ \Sol(\param)\equiv \mat{O} \left( \Sol(\param)\right)$.

    It is impossible to find a simple unique formalism of the simulation model that accurately represents
    all possible dynamical systems.
    In this work and in the implementation of {PyOED}\xspace, we differentiate two types of simulation models:
    time-independent and time-dependent simulations.
    %
    A wide range of time-dependent simulation models for dynamical systems governed by
    partial differential equations (PDEs) can be described as
    %
    \begin{equation}\label{eqn:dudt}
      \dfrac{\partial \boldsymbol{u}}{\partial t}
        = \boldsymbol{f}(\boldsymbol{u}(\boldsymbol{x}, t, \mu)),
    \end{equation}
    %
    where $\boldsymbol{u}$ represents the prognostic variable(s) (e.g., physics),
    $\boldsymbol{x}$ denotes the spatial coordinates,
    and $\mu$ defines the physics parameters of the model.
    %
    To numerically solve the simulation model $\Sol$ equations, we utilize spatial discretization and temporal integration
    routines.
    If we follow a \emph{discretize-then-optimize} approach, we can rewrite~\eqref{eqn:dudt} in the form
    %
    \begin{equation}\label{eqn:dudt_discretized}
        \x_n := \x(t_{n})

            = \Sol_{t_{n-1}\rightarrow{t_n}}(\x_0, \mu),
    \end{equation}
    %
    where $\x_n := \x(t_{n}) \equiv \x(t_{n},\boldsymbol{x}, \mu) $ is the
    model state at time instance $t_n$ and $\Sol_{t\rightarrow t+\Delta t}$ is a one-time step
    transition mapping that results from the application of a standard spatial discretization
    method (e.g., finite difference, finite volume, or finite element) and
    time integration scheme (e.g., Runge--Kutta routine) with a step size $\Delta t$.
    %
    Thus, the prediction at time $t_n$ can be related to the initial condition $\x(t_0)$ and
    the model parameters using the recursive application of the mapping $\Sol$
    over a time interval $[t_0, t_n]$ as
    %
    \begin{equation}\label{eqn:umap}
      \x_n

        = \Sol_{t_{n-1}\rightarrow {t_n}} \circ \dots \circ \Sol_{t_0 \rightarrow {t_1}}(\x_0,\mu) \,.
    \end{equation}
    %

    In space-time formulations, one can stack the model state in one long vector
    $
      \x := \left( \x_n\tran, \ldots,\x_0\tran\right)\tran \,
    $
    and define the solution operator $\Sol$ as a block operator (e.g., a block matrix) that operates
    recursively on the the components of $\x$.
    This would enable unifying the formulation---to some extent---of the forward problem into
    the form~\eqref{eqn:forward_problem_base}.

    Since the observational measurements (e.g., sensory data) are not necessarily the same as the model
    state, we define the state-to-observable mapping $\mat{O}_n(\cdot)$ to map the
    state $\x_n$ onto the observation space, for example, by restricting that state
    onto the observational grid points, as
    %
    \begin{equation}\label{eqn:ymap}
        \obs(t_n) = \mat{O}_n(\x(t_{n})) \,,
    \end{equation}
    %
    which then can be used to construct a general observational vector
    $\obs \in \Rnum^{\Nobs}$ representing spatiotemporal data, for example, by stacking observations
    at multiple time points.

    Note that while we focused the discussion above on time-dependent problems, the case of
    time-independent simulations can be thought of as a special case of~\eqref{eqn:umap} where
    $\Sol$ maps, for example, the physics parameter $\mu$ to a model state $\x$ and
    the time index is dropped.
    %
    In most inverse problem formulations, based on the application of interest, the inversion parameter
    $\iparam$ stated in~\eqref{eqn:forward_problem_base} stands for the model physics $\mu$, the model initial
    condition $\x_0$, or both.
    %
    Note also that we have omitted details including adaptive time stepping where $\Delta t$ is adaptively
    adjusted to guarantee stability and accuracy of the time integration methodology.
    We intentionally remove such details from the discussion here to simplify the presentation and
    focus more on OED.
    In fact, the OED routines in {PyOED}\xspace are designed to solve OED problems where the utility function
    is regarded as a black box, thus abstracting the OED capabilities from the inverse problem definition.
    For these reasons, we take~\eqref{eqn:forward_problem} to be an acceptable simplification that
    describes a forward problem setup that is general enough for our purposes in this paper.

    Both the simulation model $\Sol$ and the observation operator $\mat{O}$ are imperfect and generally
    include sensory noise and representativeness errors characterizing imperfection of the map between the
    model space and observation space.
    The fact that model observations $\Fcont(\iparam)$ are not perfectly aligned with observational data
    ($\obs$) is modeled---assuming additive noise---by adding the noise term $\delta$ to the
    simulated observations $\Fcont(\iparam)$.
    %
    In most applications, the observational noise follows a Gaussian distribution
    $\obsnoise \sim \GM{\vec{0}}{\Cobsnoise}$, where $\Cobsnoise$ is the observation error covariance
    matrix that captures uncertainty stemming from sensory noise and representativeness errors.
    In this case, the data likelihood is
    %
    \begin{equation} \label{eqn:Gaussian_likelihood}
    \CondProb{\obs}{ \param } \propto
      \exp{\left( - \frac{1}{2}
        \sqwnorm{ \Fcont(\param) - \obs }{ \Cobsnoise\inv } \right) } \,,
    \end{equation}
    %
    where the matrix-weighted norm in~\eqref{eqn:Gaussian_likelihood} is defined as
    $\sqwnorm{\vec{x}}{\mat{A}} = \vec{x}\tran \mat{A} \vec{x} $ for a vector $\vec{x}$ and
    a square symmetric matrix $\mat{A}$ of conformable sizes.

  \subsection{The inverse problem}
  \label{subsec:inverse_problem}

    An inverse problem refers to the retrieval of the model parameter $\param$ from noisy
    observation $\obs$, conditioned by the model dynamics.
    %
    This can be achieved by finding a point estimate or by building a complete probabilistic
    description as discussed in Section~\ref{sec:Introduction}.
    In the former, an optimization problem is solved to minimize the mismatch between the expected
    observations (through simulation models) and real data.
    This is typically employed in variational DA methods where the estimate of the true
    parameter is obtained by minimizing a regularized log-likelihood objective, where regularization
    is employed to enforce smoothness or background information on the parameter.
    In this case, a  point estimate of the true $\param$ is obtained by solving
    %
    \begin{equation}\label{eqn:fdvar}
      \argmin_{\param} \obj(\param)
        := \frac{1}{2} \sqwnorm{ \Fcont(\param) - \obs }{ \Cobsnoise\inv }
         + \frac{1}{2} \sqwnorm{ \param - \paramprior }{ \Cparamprior\inv } \,,
    \end{equation}
    %
    where $\paramprior$ is an initial guess of the unknown true value of $\param$.
    In general, the second term is added to enforce regularization or prior knowledge on the solution, for example,
    if the solution is assumed a priori to follow a Gaussian distribution $\param\sim\GM{\paramprior}{\Cparamprior}$.

    %
    Uncertainty envelopes around the single-point estimate obtained by solving~\eqref{eqn:fdvar} can be
    developed, for example, by using Laplacian approximation~\cite{tierney1986accurate} where the posterior is approximated
    by a Gaussian distribution.
    This approach has been successfully employed in infinite-dimensional Bayesian inversion
    problems~\cite{stuart2010inverse}.
    %
    A fully Bayesian approach, on the other hand, aims to provide a consistent probabilistic description
    of the unknown parameter along with the associated uncertainties and is not limited to Gaussian
    distributions.
    This is achieved by describing the posterior, that is, the probability distribution of
    the model parameter $\param$ conditioned by the available simulations and noisy data $\obs$, and is obtained by applying a form of Bayes' theorem
    \begin{equation}\label{eqn:Bayes}
        \CondProb{\param}{\obs} \propto \CondProb{\obs}{\param} \Prob(\param) \,,
    \end{equation}
    %
    where $\Prob(\param)$ is the prior, $\CondProb{\obs}{\param}$ is the data likelihood, and
    $\propto$ indicates removal of a normalizing constant in the right-hand side of~\eqref{eqn:Bayes}.
    For further details on Bayesian inversion see, for example,,~\cite{smith2013uncertainty,stuart2010inverse}.
    Given the posterior~\eqref{eqn:Bayes}, one can use the maximum a posteriori (MAP) point as an estimate
    of the true unknown QoI or follow a Monte Carlo approach to sample the posterior, thus building a complete probabilistic picture;
    see, for example,~\cite{attia2017reduced}.
    %
    Both the variational and the Bayesian inference approaches provide a plethora of techniques for
    statistical data analysis in general, and specifically for solving inverse problems.

    The Bayesian perspective provides a formal mathematical ground for estimating the physical QoI,
    for example, the model parameter $\iparam$, along with the associated uncertainties given
    the available sources of information.
    In many cases, however, this inversion is an intermediate step, and the goal QoI is a function
    of the model parameter, that is, $\pred :=\PredOper(\param)$.
    A goal-oriented approach is followed in this case where one aims to inspect the posterior
    of the QoI conditioned by the available data~\cite{LiebermanWillcox13,LiebermanWillcox14}.

    \paragraph{The ideal case: linear Gaussian problems}

      If the forward operator $\Fcont$ is linear (or linearized), and assuming Gaussian observational noise
      and a Gaussian prior
      $\GM{\iparb}{\Cparampriormat}$, then the posterior is Gaussian $\GM{\ipara}{\Cparampostmat}$ with
      %
      \begin{equation}\label{eqn:Gaussian_Posterior_Params}
        \Cparampostmat = \left(\F \adj \Cobsnoise\inv \F
          + \Cparampriormat\inv \right)\inv \,, \quad
      %
      \ipara = \Cparampostmat \left( \Cparampriormat\inv \iparb
        + \F\adj \Cobsnoise\inv \, \obs \right) \,,
      \end{equation}
      %
      where $\F \equiv \Fcont$ is the forward model and $\F\adj$ is the associated adjoint.
      %
      Despite being simple, this setup~\eqref{eqn:Gaussian_Posterior_Params} is of utmost importance in
      the Bayesian inversion and OED literature and is elementary for testing implementations of
      new DA and OED approaches, mainly because the posterior can be formulated exactly.
      %
      Moreover, in many large-scale applications, the posterior can be approximated,
      to an acceptable degree, by a Gaussian distribution obtained by
      linearizing the nonlinear operator $\Fcont$ around the MAP estimate. The linearized model is also known as the tangent linear model (TLM) obtained by differentiating $\Fcont$.

      As mentioned earlier, inversion for the parameter $\param$ is often an intermediate stage, and
      the end-goal is to describe the posterior of a general QoI that is not the model parameter $\iparam$
      but rather a goal quantity $\pred$ that depends on the inversion parameter $\iparam$.
      %
      Specifically, goal-oriented inversion seeks the posterior
      $
        \CondProb{\pred}{\obs} \propto \Like{\obs}{\pred, \iparam} \Prob(\pred)\,,
      $
      %
      where $\Prob(\pred)$ is a prior on the goal QoI and
      $
        \Like{\obs}{\iparam} = \Like{\obs}{\pred,\iparam}
      $
      is the data-likelihood~\eqref{eqn:Gaussian_likelihood}, where $\pred$ is determined completely
      based on $\iparam$.
      %
      We focus the discussion here on the case of linear prediction operators $\Predmat$.
      That is, we consider prediction quantities of the form
      %
      \begin{equation}\label{eqn:GOOED_linear_QoI}
        \pred = \Predmat \iparam,
      \end{equation}
      %
      where $\Predmat$ is a linear prediction operator.
      Within the Gaussian linear setting, the prior of the goal QoI $\pred$ is
      $\GM{\predb}{\Cpredpriormat}$ with
      %
      \begin{equation}\label{eqn:prior_prediction_PDF}
        \predb = \Predmat \iparamb,
        \qquad
        \Cpredpriormat = \Predmat \Cparampriormat \Predmat^*,
      \end{equation}
      %
      where $\Predmat^*$ is the adjoint of the prediction operator $\Predmat$.
      %
      The posterior distribution of the prediction $\pred$, conditioned by the observations $\obs$,
      is also Gaussian and is given by $\GM{\preda}{\Cpredpostmat}$, where

      \begin{equation}\label{eqn:posterior_prediction}
      	\begin{aligned}
          \preda = \Predmat \iparama\,, \qquad \Cpredpostmat = \Predmat \Cparampostmat \Predmat^*
          %
          = \Predmat\, \left( \F \adj \Cobsnoise\inv \F + \Cparampriormat^{-1} \right) \, \Predmat^*.
        \end{aligned}
      \end{equation}
      %

      Note that the goal-oriented Bayesian inversion falls back to the standard formulation of a Bayesian
      inverse problem if the prediction operator $\PredOper$ is an identity operator.

  \subsection{Optimal experimental design}
  \label{subsec:OED}

    Here we outline the basics of an OED problem for Bayesian inversion.
    An excellent review of recent advances on model-constrained OED can be found
    in~\cite{alexanderian2020optimal}. An OED optimization problem takes the general form
    %
    \begin{equation}\label{eqn:OED_optimization}
    \design\opt
    = \argmax_{\design} \,
      \mathcal{U}(\design) \,,
    \end{equation}
    %
    where $\mathcal{U}$ is a predefined utility function that quantifies the quality of the design $\design$.
    %
    The nature of  $\design$ depends on the application at hand,
    and the utility function $\mathcal{U}$ is chosen to enable defining what an ``optimal'' design is.
    %
    The optimization problem~\eqref{eqn:binary_OED_optimization} is often associated with
    an auxiliary sparsity-enforcing term $ - \regpenalty \penaltyfunction{\design}$ to prevent dense designs
    and to reduce the cost associated with deploying observational sensors.
    The utility function can be rewritten as
    %
    \begin{equation}
        \mathcal{U}(\design) = \Psi(\design) - \regpenalty \penaltyfunction{\design} \,,
    \end{equation}
    %
    where $\Psi(\cdot)$ is an OED optimality objective, referred to hereafter as the ``optimality criterion,'' which is defined based
    on the inverse problem at hand and on a chosen criterion (e.g., from the well-known OED alphabetic criteria).
    The function $\penaltyfunction{\design}$ asserts regularization or sparsity on the design.
    For example, this could be a resource constraint:
    $\sum_{i}^{\Nsens}\design_i = \wnorm{\design}{} \leq k  $\,, \,\, or
    $\sum_{i}^{\Nsens}\design_i = \wnorm{\design}{} = k \,; k\in \mathtt{Z}_{+}\,,$
    for example,  an upper bound (or exact budget) on the number of sensors (in OED).
    It also could be a sparsifying (possibly nondifferentiable) function,
    for example, $\wnorm{\design}{0}$,  $\wnorm{\design}{1}$.
    %

    Generally speaking, we seek a design that maximizes the utility function $\mathcal{U}$. Other
    formulations, however, involve minimization of an OED optimality criterion; see, for example,~\cite{attia2018goal}.
    Both formulations are adopted in the OED literature and in {PyOED}\xspace and often are even equivalent, as explained below.
    %
    In inverse problems, a design can be associated with the observational configuration
    and thus can be used to optimally select an optimal observational policy.
    For example, a design can be defined to select sensor location or temporal observation frequency
    that can help provide accurate prediction with minimum uncertainty, or it can be defined to select
    an observational configuration that guarantees maximum information gain from the data.

    \paragraph{OED for sensor placement}

      In sensor placement we associate a binary design variable $\design_i$  with the $i$th
      candidate sensor location with $1$ meaning activating the sensor and deactivating it otherwise.
      This defines the design as a binary vector $\design\equiv\binarydesign\in\{0,1\}^{\Nsens}$,
      which collectively define the observational configuration.
      In this case, the OED problem~\eqref{eqn:OED_optimization} takes the form
      %
      \begin{equation}\label{eqn:binary_OED_optimization}
        \design\opt
          = \argmax_{\design \in \{0, 1\}^{\Nsens}} \mathcal{U}(\design) :=
          \Psi(\design) - \regpenalty \penaltyfunction{\design}  \,.
      \end{equation}
      %

      In Bayesian OED for sensor placement, the observation covariance $\Cobsnoise$ is replaced with a weighted version
      $\wdesignmat(\design)$, resulting in the weighted data-likelihood

      \begin{subequations}
      %
      \begin{equation}\label{eqn:weighted_joint_likelihood}
        \Like{\obs}{ \param; \design}
          \propto \exp{\left( - \frac{1}{2} \sqwnorm{ \F(\param) - \obs}{\wdesignmat(\design)} \right) } \,,
      \end{equation}
      %
      where the weighted observation error covariance matrix takes the form
      %
      \begin{equation}\label{eqn:Shur_weighted_Precision}
      %
      \wdesignmat(\design)
        :=  \pseudoinverse{ \designmat(\design) \Cobsnoise \designmat(\design) }
        = \proj\tran(\design) \left(
             \proj(\design) \Bigl(\designmat(\design) \Cobsnoise \designmat(\design) \Bigr) \proj\tran(\design)
           \right)\inv \proj(\design) \,,
      \end{equation}
      %
      \end{subequations}

      where $\dagger$ denotes the Moore--Penrose (pseudo) inverse
      and $\designmat(\design) :=\diag{\design}$ is a diagonal matrix with the binary design
      $\design\in\{0, 1\}^{\Nsens}$ on its diagonal.
      %
      $\proj(\design)$ is a sparse matrix that extracts nonzero rows/columns from the
      design matrix $\wdesignmat$; see~\cite{attia2022optimal} for further details.

    \paragraph{The utility function }

      In the case of linear Bayesian inversion, the posterior is Gaussian with the covariance being
      independent from the actual realizations of the data, as shown by~\eqref{eqn:Gaussian_Posterior_Params}
      and~\eqref{eqn:posterior_prediction}.
      This fact enables designing observational policies before actually deploying the observational sensors.
      Specifically, in linear Bayesian OED, we set the objective to minimize a scalar summary of the posterior uncertainty,
      that is, the posterior covariance matrix.
      %
      This is the underlying principle of the alphabetical criteria~\cite{Pukelsheim93}.
      For example, an A-optimal design is the one that minimizes the trace of the posterior covariance matrix,
      and a D-optimal design is the one that minimizes its determinant (or equivalently the log-determinant).
      Note that in the case of a linear model $\F$, the Fisher information matrix $\mathsf{FIM}$ is equal to
      the inverse of the posterior covariance matrix, that is,
      %
      $
      \mathsf{FIM} = \Cparampostmat\inv(\design)
        = \F\adj \wdesignmat(\design) \F + \Cparampriormat\inv \,.
      $

      %
      Thus, in this case the utility function---discarding the penalty term---is set to
      $\mathcal{U}(\design):=\Trace{\mathsf{FIM}(\design)}$ for
      A-optimal designs and $\mathcal{U}(\design):=\logdet{\mathsf{FIM}(\design)}$ for
      D-optimal designs, and then the utility function is maximized.

      When the model $\Fcont$ is nonlinear, the $\mathsf{FIM}$ requires evaluating the
      TLM at the true parameter, that is,
      $\F=\partial\Fcont|_{\iparam=\iparam_{\rm true}}$. Thus, to obtain an optimal design,
      one can iterate over finding the MAP estimate of $\iparam$ and solving an OED
      problem with Gaussian approximation around that estimate.
      Other utility functions employed in nonlinear OED problems or non-Gaussian distributions include
      the Kullback--Leibler divergence between the posterior and the prior~\cite{huan2014gradient}.

    \paragraph{Popular solution approaches}

      The OED problem~\eqref{eqn:binary_OED_optimization} can be viewed as a mixed-integer program
      and can be solved by using branch-and-bound~\cite{Gerdts05,Leyffer01}. However, this type of research has not yet been applied to model-constrained-OED.
      A common approach to solving~\eqref{eqn:binary_OED_optimization} is to replace the binary optimization
      with the following relaxation:
      %
      \begin{equation}\label{eqn:relaxed_oed}
        \design\opt
          = \argmax_{\design \in [0, 1]^{\Nsens}} \,
        \mathcal{U}(\design) \,,
      \end{equation}
      %
      where the design variables are relaxed to take any values in the interval $[0, 1]$ rather than
      only the binary values $\{0, 1\}$. This approach, if carried out properly, has the effect of generating
      a continuous relaxation surface that connects the values of the objective evaluated at the
      binary designs; see~\cite{attia2022optimal} for further details.
      %
      A gradient-based optimization approach is generally used to solve~\eqref{eqn:relaxed_oed}, which requires
      developing the gradient of both the optimality criterion $\Psi$ and the penalty $\penaltyfunc$
      with respect to the design $\design$.

      As mentioned earlier, the penalty term is generally chosen to promote sparsity on the design.
      A popular penalty function is based on the $\ell_0$ norm to promote design sparsity
      and is thus nonsmooth and consequently is nondifferentiable.
      This difficulty can be alleviated, for example, by approximating the effect of $\ell_0$
      with a sequence of differentiable functions that converge in effect to the $\ell_0$;
      see, for example,~\cite{AlexanderianPetraStadlerEtAl14}.
      %
      In order to guarantee continuity of the relaxation surface, the weighted precision matrix is defined in
      the general form
        %
        \begin{equation}\label{eqn:pointwise_weighted_precision}
        \begin{aligned}
          \wdesignmat(\design)&:=\pseudoinverse{\designmat(\design) \odot \Cobsnoise}\,, \\
          \designmat_{i,j}(\design)
            &:=
            \begin{cases}
              \weightfunc_i\, \weightfunc_j &;\,  i \neq j \\
              \begin{cases}
                0     &; \, \weightfunc_i =0 \\
                \frac{1}{\weightfunc_i^2}  &; \weightfunc_i \neq 0
              \end{cases} &;\, i=j \\
            \end{cases} \,;\,\,
              \begin{matrix}
                & i,j = 1,2,\ldots,\Nsens \,,  \\[-1.0pt]
                & m,n = 1,2,\ldots,\nobstimes  \,,
              \end{matrix}
          \end{aligned}
        \end{equation}
        %
        where $\odot$ is the Hadamard (Schur) product of matrices and $\weightfunc_i\in[0, 1]$ is a weight
        calculated by using $\design_i$,
        for example, $\weightfunc_i:=\design_i$; see~\cite{attia2022optimal} for additional details.
        %
        The formulation~\eqref{eqn:pointwise_weighted_precision} guarantees that for $\design_i\in[0, 1]^{\Nsens}$ it holds
        that
        $
          \lim_{\design \rightarrow \binarydesign} \wdesignmat(\design)
            = \wdesignmat(\binarydesign)
        $
        for a binary design $\binarydesign\in\{0, 1\}^{\Nsens}$ and thus guarantees continuity of the relaxation surface.
        %
        Thus, the solution of the relaxed OED optimization problem~\eqref{eqn:relaxed_oed} is guaranteed to match
        the solution of the original binary OED optimization problem.

        The A- and D-optimal design relaxed optimization problems~\eqref{eqn:relaxed_oed} take the following
        respective forms:
        %
        \begin{subequations}\label{eqn:A_and_D_optimality_Schur}
        \begin{align}
        \design^{\rm A\!-\!opt}
          &= \argmax_{\design\in [0, 1]^{\Nsens}}
            := \Trace{\Predmat \left( \F\adj
              \pseudoinverse{ \Cobsnoise  \odot \designmat(\design) }
            \F
              + \Cparampriormat\inv \right)\inv \Predmat\adj
            }\inv
            - \regpenalty \penaltyfunction{\design}
            \,, \label{eqn:A_optimality_Schur}\\
        %
        \design^{\rm D\!-\!opt}
          &= \argmax_{\design\in [0, 1]^{\Nsens}}
            := \logdet{\Predmat \left( \F\adj
              \pseudoinverse{ \Cobsnoise  \odot \designmat(\design) }
            \F
              + \Cparampriormat\inv \right)\inv \Predmat\adj
            }\inv
            - \regpenalty \penaltyfunction{\design}
            \,.\label{eqn:D_optimality_Schur}
        \end{align}
        \end{subequations}
        %

        The most important piece of information for solving the relaxation~\eqref{eqn:relaxed_oed}
        is the gradient of the utility function; see, for example,~\cite{attia2022optimal} for a detailed derivation of
        the gradients of the objective functions in~\eqref{eqn:A_and_D_optimality_Schur}.
        Gradient formulation, however, is mathematically involved and can be extremely computationally demanding
        because it requires numerous evaluations of the forward operator, the goal operator, and the corresponding adjoints.
        %
        Moreover, the penalty function $\penaltyfunction{\cdot}$ is required to be differentiable.

        A stochastic learning approach to binary OED has been recently presented in~\cite{attia2022stochastic},
        to solve the binary optimization problem~\eqref{eqn:binary_OED_optimization} without the need for relaxation.
        This approach does not require differentiability of the utility function $\mathcal{U}$.

        In this approach the optimal design is defined as
        %
        \begin{equation}\label{eqn:stochastic_OED_optimization}
          \hyperparam\opt
            = \argmax_{\hyperparam\in[0,1]^{\Nsens}}
              \Expect{\design\sim\CondProb{\design}{\hyperparam}}{\mathcal{U}(\design) - b} \,,
        \end{equation}
        %
        where $\CondProb{\design}{\hyperparam}$ is a multivariate Bernoulli distribution with parameter
        $\hyperparam$ specifying probabilities of success/activation of each entry of $\design$,
        that is, $\hyperparam_i\in[0, 1]$.
        Here $b$ is a constant ``baseline'' used to minimize variability of the stochastic estimate
        of the gradient; see~\cite{attia2022stochastic} for further details.
        Algorithm~\ref{alg:REINFORCE_baseline} summarizes the procedure followed to
        solve~\eqref{eqn:stochastic_OED_optimization}.

        \begin{algorithm}[htbp!]
          \caption{Stochastic optimization for binary OED with the optimal baseline.}
          \label{alg:REINFORCE_baseline}
          \begin{algorithmic}[1]

            \Require{Initial distribution parameter $\hyperparam^{(0)}$,
                      step size schedule $\eta^{(n)}$,
                      sample sizes $\Nens,\, m$,
                      baseline batch size $b_m$}
            \Ensure{$\design\opt$}

            \State{initialize $n = 0$}

            \While{Not Converged}
              \State {
                Update $n\leftarrow n+1$
              }

              \State{
                Sample $\{\design[j]; j=1,2,\ldots,\Nens \}\sim
                  \CondProb{\design}{\hyperparam^{(n)}}$
              }

              \State{
                Calculate  $b$ = \Call{OptimalBaseline}{$\hyperparam^{(n)}$,
                  $\Nens$, $b_m$}
              }

              \State\label{algstep:REINFORCE_baseline:grad}{
                Calculate $ \vec{g}^{(n)}=\frac{1}{\Nens} \sum_{j=1}^{\Nens}
                  \left(\obj(\design[j]-b)\right) \sum_{i=1}^{\Nsens}
                  \left(
                    \frac{\design_i[j]}{\hyperparam_i}
                      + \frac{\design[j]_i-1}{1-\hyperparam_i}
                  \right) \,\vec{e}_i $
              }

              \State \label{algstep:REINFORCE_baseline:proj}{
                Update $\hyperparam^{(n+1)}
                  = \Proj{}{\hyperparam^{(n)} - \eta^{(n)} g^{(n)} } $
              }

            \EndWhile

            \State{
              Set $\hyperparam\opt = \hyperparam^{(n)}$
            }

            \State\label{algstep:REINFORCE_baseline:sampling} {
              Sample $\{\design[j];j=1,2,\ldots,m \} \sim
                \CondProb{\design}{\hyperparam\opt}$,
              and calculate $\obj(\design[j])$
            }

            \State \Return{
                $\design\opt $: the design $\design$ with smallest
                  value of $\obj$ in the sample.
            }

            \Function{OptimalBaseline}{$\theta$, $\Nens$, $b_m$}

              \State{Initialize $b \gets 0$}

              \For {$e$ $\gets 1$ to $b_m$}

                \For {$j$ $\gets 1$ to $\Nens$}
                  \State{
                    Sample $\design[j] \sim \CondProb{\design}{\hyperparam}$
                  }
                  \State\label{algstep:REINFORCE_baseline:logprob}{
                    Calculate $\vec{r}[j] = \sum_{i=1}^{\Nsens}
                      \left(
                        \frac{\design_i[j]}{\hyperparam_i}
                          + \frac{\design[j]_i-1}{1-\hyperparam_i}
                      \right) \,\vec{e}_i $
                  }
                \EndFor

                \State{
                  Calculate $\vec{d}[e] = \frac{1}{\Nens} \sum_{j=1}^{\Nens}
                  \vec{r}[j] \,,\,$
                  and

                  $\, \vec{g}[e] = \frac{1}{\Nens} \sum_{j=1}^{\Nens}
                    \obj(\design[j]) \, \vec{r}[j] $
                }
                \State{
                  Update $b \gets b + \left(\vec{g}[e] \right)\tran \vec{d}[e]$
                }
              \EndFor

              \State{
                Update
                  $b \gets

                    b \, \Nens \, / \left( b_m \,
                     \sum_{i=1}^{\Nsens} \frac{1}{\hyperparam_i- \hyperparam_i^2} \right)
                  $
                    %

              }

              \State \Return $b$
            \EndFunction

          \end{algorithmic}
        \end{algorithm}

        %
        Note that we do not provide an exclusive set of formulations or solution approaches in this study. We provide here only an exemplary set of formulations and algorithms used to inspire the development of {PyOED}\xspace, which itself can be used to test further formulations and algorithmic approaches.

\section{{PyOED}\xspace: Structure and Philosophy}
\label{sec:PyOED}
  {PyOED}\xspace aims to provide a unified platform for implementing and testing model-constrained OED algorithmic approaches
  including formalisms~\eqref{eqn:binary_OED_optimization},~\eqref{eqn:relaxed_oed},
  ~\eqref{eqn:A_and_D_optimality_Schur}, and~\eqref{eqn:stochastic_OED_optimization}.
  %
  Solving model-constrained OED and inverse problems requires proper understanding and formulation of
  the underlying dynamical system, observational configuration, uncertainty models,
  DA and inversion algorithms, and the OED objective~\cite{ucinski2000optimal} and the selected utility function.
  %
  {PyOED}\xspace is a stand-alone, yet extensible, Python package that provides users and researchers in
  the computational science and engineering disciplines with a testing suite that effectively
  glues these components in an OOP fashion.
  %
  For example, {PyOED}\xspace provides a variety of time-dependent and time-independent simulation models.
  These include systems governed by linear algebraic equations, ordinary differential equations,
  and PDEs.
  {PyOED}\xspace is also equipped with a set of classes implementing various observational operators,
  probabilistic uncertainty models, and DA and OED methods.
  A high-level overview of the {PyOED}\xspace major components and their coupling for solving DA and OED problems
  is provided in Figure~\ref{fig:DA_and_OED}.
  In Section~\ref{subsec:code_structure} we briefly describe the main components of {PyOED}\xspace
  and outline the functionality they provide in correspondence with the diagram~\ref{fig:DA_and_OED}.
  %
  \begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.75\textwidth]{Plots/DA_and_OED}
    \caption{High-level overview of the main components of {PyOED}\xspace.}
    \label{fig:DA_and_OED}
  \end{figure}
  %

  \subsection{Code structure}
  \label{subsec:code_structure}

    Figure~\ref{fig:PyOED_Subpackages} shows the main subpackages (ordered alphabetically)
    shipped with the current version of {PyOED}\xspace (v1.0).
    %
    The rest of this section~\ref{subsec:code_structure} provides a high-level description of the packages/subpackages
    of~{PyOED}\xspace as displayed in Figure~\ref{fig:PyOED_Subpackages}.
    %
    \begin{figure}[!htbp]
      \centering
      \includegraphics[width=0.45\textwidth]{Plots/PyOED_Subpackages}
        \caption{Main subpackages (ordered alphabetically) available in the current version of {PyOED}\xspace (v1.0).}
      \label{fig:PyOED_Subpackages}
    \end{figure}
    %

    \paragraph{\code{pyoed.models}}

      Following the convention in the {DATeS}\xspace package~\cite{attia2019dates}, we use the word ``model'' to refer
      to three entities:  \emph{the simulation model}, \emph{the observation model (or operator)}, and \emph{the error models}.

      The simulation model provides a prediction about the behavioral pattern of the physical phenomena of concern.
      In the current version of {PyOED}\xspace (v1.0) we provide various simulation models under the
      \code{pyoed.models.simulation\_models}, including several versions of the Lorenz system~\cite{lorenz1996predictability} and advection-diffusion models. The structure of these prototypical simulation models should provide clear guidelines to practitioners willing to adopt {PyOED}\xspace for their particular applications.

      The observation operator maps the model state onto the observation grid, thus providing a functional mapping
      between the model state and observational data.
      Two of the most prominent observation operators in experimental settings are the identity operator and an interpolator.
      {PyOED}\xspace provides implementations of several observation operators including these two, with an observational design properly incorporated to enable
      altering observational configurations at any point in the DA or OED solution process.
      Observation operators are provided in the \code{pyoed.models.observation\_operators} subpackage.

      The error models quantify the uncertainty associated with the model parameter, model state, and observational data.
      An experimental design can be associated with any of these pieces.
      For example, in sensor placement, an experimental design is associated with the observational grid; thus, modifying the
      observational design affects the observational error model.
      For example, in the relaxation approach~\eqref{eqn:A_and_D_optimality_Schur}, the design weights scale the entries of the covariance matrix,
      and the stochastic approach~\eqref{eqn:stochastic_OED_optimization} works by removing rows/columns of the observation
      error covariance matrix corresponding to zero design variables.
      %
      {PyOED}\xspace provides various implementations of error models suitable for modeling priors, as well as  observational errors in
      Bayesian inversion, where a design variable is consistently implemented to enable modifying the experimental design during any
      step of the DA and/or OED solution process.
      The current version of {PyOED}\xspace (v1.0) provides various error model implementations through
      the subpackage \code{pyoed.models.error\_models}, including a Gaussian model and Laplacian model.

    \paragraph{\code{pyoed.assimilation}}

      {PyOED}\xspace provides a set of DA tools that include algorithms for ``filtering''
      and ``smoothing.'' These two terms are widely used in the DA literature.
      The former algorithm solves inverse problems that involve time-independent or time-dependent
      simulation models, while the latter algorithm is restricted to time-dependent models.
      Filtering involves prediction (of observation) using the parameter-to-observable map,
      followed by a correction procedure to correct knowledge of the QoI given the observational data.
      In filtering for time-dependent simulations, the observational data is assimilated sequentially, with one observation time per assimilation window/cycle.
      Examples of filtering DA methods include three-dimensional variational  DA, and Kalman filtering~\cite{evensen2009data,asch2016data}.
      %
      Smoothing, on the other hand, is concerned with history matching; these algorithms try
      to find the QoI that best matches multiple spatiotemporal observations (a trajectory)
      and is usually defined as an initial value problem.

      Examples include space-time Bayesian inversion~\cite{Stuart10}
      and four-dimensional variational DA~\cite{asch2016data},
      for which vanilla implementations are provided in~{PyOED}\xspace.
      Implementations of filtering DA algorithms are provided through \code{pyoed.assimilation.filtering},
      and smoothing algorithms are provided in \code{pyoed.assimilation.smoothing}.

    \paragraph{\code{pyoed.optimization}}

      Numerical optimization routines are elementary for solving OED optimization problems, as well as the
      variational approaches for solving DA problems.
      A variety of optimization software packages  can be used for solving numerical optimization problems
      including those  described in this work.
      {PyOED}\xspace enables using external optimization packages, including  Python's \code{Scipy} package,
      to solve DA and OED optimization problems.
      {PyOED}\xspace, however, provides specific implementations of optimization procedures not available in popular optimization
      packages, such as the stochastic algorithm described by Algorithm~\ref{alg:REINFORCE_baseline},
      various versions of the stochastic average approximation (SAA) algorithm, and robust optimization~\cite{attia2023robust}.

    \paragraph{\code{pyoed.ml}}

      This subpackage is intended to provide implementations of machine learning  algorithms useful for DA and OED applications.
      For example, the stochastic learning approach to OED~\eqref{eqn:stochastic_OED_optimization} can be seen as
      a reinforcement learning (RL) approach to solving the OED problem.
      %
      The module \code{pyoed.ml.reinforcement\_learning} under this package provides implementation of RL components,
      including an agent, a policy, transition probability, actions, and utility functions.

    \paragraph{\code{pyoed.stats}}

      This package aims to collect statistical procedures used by other parts of the package, such as
      sampling routines, and implementation of random variables and their probabilistic utility functions including
      density evaluation and log-probabilities.
      This version of {PyOED}\xspace (v1.0) provides an exemplary implementation of a multivariate Bernoulli distribution required by
      the RL algorithm~\ref{alg:REINFORCE_baseline}.
      Since statistical tools are crucial for various DA and OED algorithms, we chose to keep the subpackage \code{pyoed.stats} rather
      than moving these implementations to other parts of the package.

      This approach is advantageous because we continuously extend the package with various statistical tools, for example,
      for randomized approximation methods for Bayesian inversion.

    \paragraph{\code{pyoed.oed}}

      OED is the main component of {PyOED}\xspace that provides implementations of various
      algorithmic approaches for solving OED problems, including relaxation~\eqref{eqn:A_and_D_optimality_Schur} and
      stochastic learning~\eqref{eqn:stochastic_OED_optimization},
      as well as recent developments including robust OED~\cite{attia2023robust}.
      Most implementations in this package take an inverse problem (DA object) as input and use it to access all
      the underlying components, thus gaining access to the simulation model, error models, and observation operator as well
      as the experimental design.
      This approach enables the user to modify an experimental design, solve the DA problem if needed, and solve the underlying OED optimization problem.
      The core OED functionalities in most {PyOED}\xspace routines, however, can be used with black-box utility functions, waiving the need for an inverse problem if needed.

    \paragraph{\code{pyoed.utility}}

      This subpackage aims to collect implementations of general-purpose functionality, such as file I/O and visualization, as well as
      general mathematical and statistical procedures.
      The subpackage includes  matrix-free implementations of expensive operations such as evaluating the trace and log-determinant of a matrix.
      It also provides routines to approximate matrix trace using statistical randomization~\cite{AlexanderianSaibaba17,saibaba2016randomized}.

    \paragraph{\code{pyoed.examples}}

      This subpackage provides various example scripts that users can follow to learn how to effectively use various pieces of the package.
      The modules in this subpackage explain how to load all pieces of the subpackage independently
      and explain how to properly coordinate these components to design a consistent DA and/or OED experiment.

    \paragraph{\code{pyoed.tutorials}}

      Given the popularity of Jupyter Notebooks in the computational science community,
      we converted some of the examples in the subpackage \code{pyoed.examples} to Jupyter Notebooks
      and provided them in this subpackage \code{pyoed.tutorials}. We employ them in the test cases presented in Section~\ref{sec:Test_cases}.
      We can add more tutorials on reasonable demand.

  \subsection{{PyOED}\xspace utilization workflow}
  \label{subsec:workflow}

    While the components of DA and OED problems can be used independently of each other, some level
    of ordering is mandatory for proper utilization.
    For example, an inverse problem (DA) object cannot be instantiated before a simulation model,
    an observation operator, and error model objects.
    Similarly, an OED problem for Bayesian inversion cannot be solved before creating an
    inverse problem.
    The general workflow for utilizing {PyOED}\xspace components is displayed in
    Figure~\ref{fig:PyOED_workflow}.
    A practical guide that illustrates how to follow this simple workflow is described in
    Section~\ref{sec:Test_cases}.
    %
    \begin{figure}[!htbp]
    \centering
      \includegraphics[width=0.45\textwidth]{Plots/PyOED_workflow}
      \caption{Workflow  describing initialization order and access level of {PyOED}\xspace components.}
      \label{fig:PyOED_workflow}
    \end{figure}
    %

  \subsection{Extending and contributing to {PyOED}\xspace}
  \label{subsec:extending_PyOED}

    {PyOED}\xspace is meant to be extensible. Thus, we continuously interface other software tools that
    provide efficient implementations of the components of DA and OED problems.
    For example, hIPPYlib~\cite{villa2018hippylib} is a software package for solving
    high-dimensional inverse problems following an optimize-then-discretize approach.
    It has been employed to empirically verify several inversion and
    OED algorithmic developments recently.
    Instead of rebuilding the functionality of hIPPYlib and similar packages,
    we have interfaced with some of its components to show how
    easily and efficiently {PyOED}\xspace can extend other successful packages.
    For example, {PyOED}\xspace interfaces with finite-element (FE) implementations of advection-diffusion
    and Poisson models from hIPPYlib as well-as point observation operators.
    \commentout{
        Some modifications to the interfaced hIPPYlib components, however, had to be
        introduced to make the code more general rather than
        being tied to a specific inverse problem formulation.
    }
    However, such extension does not hinder the functionality or limit the extensibility of {PyOED}\xspace.
    Specifically, interfacing with such external packages is optional and is not provided in the core of {PyOED}\xspace,
    mainly because the backbone of these packages is not guaranteed not to be quickly outdated or be unmaintained.
    Thus, these extensions (e.g., interfacing with hIPPYlib) are made optional during the import process of {PyOED}\xspace subpackages,
    and dependent functionality is used only when properly installed and available on the current architecture.

  \subsection{Code availability}
  \label{subsec:code_availability}

    The development version of {PyOED}\xspace is available from \url{https://gitlab.com/ahmedattia/pyoed}.
    \commentout{
        As soon as license is granted by ANL, update this section.
        Software location, license, download from repository,
        and documentation webpage, contact, support and collaboration.
    }

\section{Test Cases}
\label{sec:Test_cases}
  {PyOED}\xspace comes with a set of prototypical test problems with increasing complexity for both DA and OED.
  An ideal case typically used in scientific publications is the linear Gaussian setup,
  where the simulation model and the observation operator are both linear
  and the error models (observation noise and the prior) are both Gaussian; see Section~\ref{sec:Background}.
  In this case, the posterior is also a Gaussian; the solution of the inverse
  problem is unique---the posterior mean and mode (MAP) are identical;
  and posterior moments (mean and covariance) both have closed forms that can be obtained by applying the Kalman filter theory.
  Such a simplified setup can be used for testing new formulations in both DA and OED, and thus it is provided in {PyOED}\xspace.
  We discuss this formulation and in Section~\ref{subsec:ToyLinear} show how it can be utilized.
  In Section~\ref{subsec:results:standard_oed_examples} we discuss in further detail a standard experiment
  widely used in OED scientific research and offered by {PyOED}\xspace.

  \subsection{An ideal setup: linear Gaussian toy problem}
  \label{subsec:ToyLinear}

    Consider a time-dependent forward problem defined
    at time instances $t_0+ i \Delta t, i=0, 1, \ldots, \nobstimes$,
    for a fixed step size $\Delta t$,  as follows:
    %
    \begin{align}
        \x_n = \mat{A}\, \x_{n-1}\,, \qquad
        \obs_n = \mat{I} \x_n + \obsnoise\,,\quad n=1, 2, \ldots ,
    \end{align}
    %
    where $\x_{n}\in\Rnum^{\Nstate}$ is the discrete model state at time instance $t_n$,
    $\mat{A}\in \Rnum^{\Nstate \times \Nstate}$ is a matrix representing model evolution over
    time interval $[t_{n-1}, t_{n}]$, and $\mat{I}$ is the identity observation operator/matrix.
    %
    If we assume $\obsnoise\sim \GM{\vec{0}}{\mat{R}}$ and assume
    a Gaussian prior $\x_{0}\sim\GM{\x_{0}^{\rm pr}}{\mat{\Gamma}_{\rm pr}}$,
    then the posterior is  Gaussian $\GM{\x_{0}^{\rm post}}{\mat{\Gamma}_{\rm post}}$ with
    %
    \begin{equation}\label{eqn:Gaussian_Posterior_Toy}
      \mat{\Gamma}_{\rm post} = \left(\mat{A} \tran \mat{R}\inv \mat{A}
        + {\mat{\Gamma}_{\rm pr}}\inv \right)\inv \,,
      \qquad
      %
        \x_{0}^{\rm post} = \mat{\Gamma}_{\rm post}
            \left( {\mat{\Gamma}_{\rm pr}}\inv \x_{0}^{\rm pr}
                + \sum_{i=1}^{\nobstimes}{\mat{A}\tran \mat{R}\inv \, \obs }
                \right) \,.
    \end{equation}
    %

    Since~\eqref{eqn:Gaussian_Posterior_Toy} is a closed form of the posterior,
    we can use it to test and debug new DA and OED implementations.
    This fact is highly utilized in the unit tests developed in {PyOED}\xspace.
    %
    To create a proper experiment, we will follow the workflow described by
    Figure~\ref{fig:PyOED_workflow}.
    %
    In the rest of this section (\ref{subsec:ToyLinear}) we describe how
    to initialize an inverse problem in {PyOED}\xspace with the
    settings~\eqref{eqn:Gaussian_Posterior_Toy}, and we provide a simple scheme
   that can be followed to initialize other experiments.
    The code summarized here is provided in the \code{pyoed.examples.fourDVar\_driver} module with additional comments,
    details, and capabilities that can help the user understand the workflow for creating and solving an inverse problem.
    A Jupyter Notebook \code{pyoed.tutorials.toy\_linear} is also available and can be used to regenerate the numerical
    results presented in this section (\ref{subsec:ToyLinear}).

    \paragraph{Creating the models}

    Assuming {PyOED}\xspace is already in the Python path, the first step is to import/load
    the simulation model (that describes $\mat{A}$), the observation operator
    (here an identity operator), and the error models to create the prior and
    the observation error model.
    This can be done as described in the code snippet~\ref{snipp:toy_linear_models_import}.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_models_import},
        firstline=2,lastline=4,
        caption={Import essential modules to create the simulation model object,
            the observation operator, the prior, and the observation error model.}]{./Code_Snippets/toy_linear.py}

    Note that we have  imported only the classes we need in this example. However, {PyOED}\xspace provides several
    other implementations of the simulation models, observation operators, and error models.
    In order to create a simulation model, an object of \code{ToyLinearTimeDepndent}
    is instantiated as described by snippet~\ref{snipp:toy_linear_simulation_model}.
    %
    This generates an internal two-dimensional array of size $5 \times 5$ that represents the forward model
    $\mat{A}$, which integrates the model state forward by a timestep $dt=0.1$.
    That internal array can be reproduced by setting the \code{random\_seed} parameter
    in the passed configurations dictionary.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_simulation_model},
        firstline=8,lastline=8,
        caption={Instantiate the simulation model object.}]{./Code_Snippets/toy_linear.py}

    Since each simulation model has its own configurations, many of which are
    assigned default values, we follow the strategy of {DATeS}\xspace~\cite{attia2019dates}
    and use dictionaries to pass model arguments. {PyOED}\xspace aggregates and validates
    the passed dictionary against the default values and initiates the model accordingly.
    For example, in snippet~\ref{snipp:toy_linear_simulation_model} we
    specify a \code{random\_seed} argument that guarantees reproducibility of
    any randomly generated data inside the model object.
    This is done by keeping an internal random state inside the model object that
    is independent from other objects and is initialized to the passed random seed.
    Thus, if no random seed is passed, each time the same model object is instantiated,
    completely random sequences will be generated if requested.
    Implementations of simulations model must provide an implementation of a class-level
    method \code{get\_default\_configs} that returns a dictionary with all default
    values used if not passed upon instantiation.
    In order to ensure that, all {PyOED}\xspace simulation models inherit the class \code{pyoed.models.simulation\_models.SimulationModel} that guarantees enforcing the implementation of mandatory methods required for the seamless integration of various components in {PyOED}\xspace.
    The final configurations of a simulation model is a combination of those in
    the passed configurations dictionary and the default values, with precedence
    given to the passed configurations.
    The simulation model's configurations (a copy of it, in fact) can be accessed through
    the attribute \code{configurations}. We generally choose to return a copy to guarantee that all settings are validated before modification.
    For example, one cannot modify the time step $dt$ without verifying whether the
    timestepping implementation is tied to that time step or not and updating dependencies accordingly.

    \paragraph{The observation operator}

    Similar to simulation models, an observation operator is created by passing the
    settings in the configurations dictionary to the observation operator class constructor as shown in
    snippet~\ref{snipp:toy_linear_obs_oper}. Thus, the observation operator
    has access to the model grid and other useful attributes to create and manipulate
    data structures (such as the model state) without having to provide any new
    implementations. This is mainly because observations in this case are the same as
    the corresponding model states (discarding observation noise).

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_obs_oper},
        firstline=11,lastline=11,
        caption={Create an identity observation operator.}]{./Code_Snippets/toy_linear.py}

    \paragraph{The prior and the data noise models}

    The error models (prior and observation noise) are created in this example as shown in
    snippet~\ref{snipp:toy_linear_error_models}.
    Note that the \code{random\_seed} configurations variable is used to set the random number generator for reproducibility.
    This enables us regenerate a set of experiments and generate proper benchmarks for a fair comparison
    between various implementations. As mentioned earlier, if no random seed is passed, each instance of
    the error model is assigned a randomly generated state that guarantees
    that each instance has its own different sequence of random numbers/vectors realizations.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_error_models},
        firstline=14,lastline=15,
        caption={Create the prior and the observation error model.}]{./Code_Snippets/toy_linear.py}

    \paragraph{The inverse problem (DA) object}

    The next step is to put these models together in action and use them to create an
    inverse problem.
    We illustrate the utilization of a DA object to solve the inverse problem
    following a 4DVar formulation.
    The literature provides a plethora of variants of the general 4DVar scheme.
    {PyOED}\xspace provides a few implementations; however, the most basic (vanilla)
    implementation is used here for illustration.
    Two approaches are followed in {PyOED}\xspace for instantiating a DA object. The first is to pass
    all configurations (upon initialization) in the configurations dictionary
    \code{configs} similar to the case of simulation models, error models, and
    observation operators. The second approach is to use the proper registration methods
    associated with the created object after instantiation.
    The latter can be also used to update components of the DA
    object after initialization. For example, one might want to change the settings of the
    assimilation time window, register new observations or remove the
    old ones, or modify or even change the prior.
    Since the first approach has already been explained with the simulation and
    error models, we demonstrate the second approach here.
    Specifically, a 4DVar assimilation object is created as in
    snippet~\ref{snipp:toy_linear_4dvar}.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_4dvar},
        firstline=18,lastline=22,
        caption={Create the inverse problem object with default settings,
          and then add (register) all the pieces created above, that is, the simulation model,
          the observation operator, the prior, and the observation error model.}]{./Code_Snippets/toy_linear.py}

    The next step is to register observational data (along with observation times).
    A standard strategy in experimentation is to create synthetic data from
    a ground truth (known as a twin experiment). This is explained by snippet~\ref{snipp:toy_linear_synthethize},
    where we define the assimilation timespan (window) to be the interval $[0, 0.3]$,
    and the observations are taken at $3$ time instances $0.1, 0.2, 0.3$.
    The observational data is mimicked by adding random noise (using the observation error model) to the observed ground truth at the corresponding observation time instance.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_synthethize},
        firstline=24,lastline=36,
        caption={Create synthetic noisy observations, and register all observation time points
          and observational data to the inverse problem object.}]{./Code_Snippets/toy_linear.py}

    The final step in the DA procedure is to solve the inverse problem and assess the quality of the solution.
    For this setup, we know the ground truth, and thus one can evaluate the root mean squared error (RMSE),
    which is a standard error metric in statistics in the DA literature.
    In order to solve the inverse problem, the \code{solve\_inverse\_problem} method of the 4DVar
    DA object is called (snippet~\ref{snipp:toy_linear_solve_inversion}).
    This method will raise an instructive error if any of the essential elements, for example, the simulation model, are not registered properly. Note that this function is flexible and allows the posterior covariance to be constructed if needed. It also allows waiving finding the MAP estimate, which can be advantageous in OED applications due to associated computational savings. For example, one might want to estimate the posterior covariance in the linear Gaussian case without evaluating the MAP.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_solve_inversion},
        firstline=39,lastline=39,
        caption={Solve the inverse problem.
          The optimization initial point is set by default to the prior mean; however, it can be modified if a better initial guess is known.
            Here, the posterior covariance is evaluated, and consequently the posterior
            is updated with both the mean and the covariance matrix.}]{./Code_Snippets/toy_linear.py}

    {PyOED}\xspace provides several utility functions to evaluate statistics, such as the RMSE, which can be used to quantify the accuracy of the inverse problem solution.
    Snippet~\ref{snipp:toy_linear_rmse} shows how to call the utility function
    \code{calculate\_rmse} and use it to evaluate the prior and the analysis (posterior)
    RMSE, which are then printed.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_rmse},
        firstline=42,lastline=46,
        caption={Calculate and print the RMSE values associated with the prior mean (initial guess here) and the posterior mean.}]{./Code_Snippets/toy_linear.py}

    The same procedure can be easily followed to inspect the RMSE results over the whole
    assimilation timespan as described by Snippet~\ref{snipp:toy_linear_rmse_traject} with
    results plotted in Figure~\ref{fig:toy_linear_rmse_traject}.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_rmse_traject},
        firstline=49,lastline=53,
        caption={Generate RMSE over the whole assimilation window.}]{./Code_Snippets/toy_linear.py}

    %
    \begin{figure}[!htbp]
      \centering
        \includegraphics[width=0.45\textwidth]{Plots/ToyLinear4DVar_RMSE}
        \caption{RMSE results of the solution of the inverse problem presented  in Section~\ref{subsec:ToyLinear}.
          The RMSE results of both the prior and the posterior trajectories plotted here are
          obtained by running the code in snippet~\ref{snipp:toy_linear_rmse_traject}.
          }
      \label{fig:toy_linear_rmse_traject}
    \end{figure}
    %

    One can also analyze the posterior covariances, for example, by generating and plotting the posterior covariance matrix.
    Given the linear Gaussian settings in the present setup, one can validate the generated posterior covariance matrix against
    the exact formula~\eqref{eqn:Gaussian_Posterior_Toy}.
    One way to construct the posterior covariance matrix is to invoke the posterior model
    as shown in snippet~\ref{snipp:toy_linear_posterior_covariance}.

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_posterior_covariance},
        firstline=74,lastline=74,
        caption={Construct and retrieve the posterior covariance matrix}]{./Code_Snippets/toy_linear.py}

    Note, however, that one should avoid constructing the covariance matrix for high-dimensional
    error models.
    Alternatively, matrix-free implementations of covariance (and precision)
    matrix-vector product should be used.
    For example, to multiply the prior covariance by \code{state},
    one should call \code{prior.covariance\_matvec(state)}. The error models provide many
    attributes to efficiently access the statistics of the model, such as covariance diagonal, and trace.
    The posterior covariance matrix, constructed by employing the posterior functionality
    as in snippet~\ref{snipp:toy_linear_posterior_covariance}
    and the covariance matrix evaluated by applying~\eqref{eqn:Gaussian_Posterior_Toy}
    along with the mismatch errors
    are plotted in Figure~\ref{fig:toy_linear_posterior_covariance}.
    %
    \begin{figure}[!htbp]
      \centering
        \includegraphics[width=0.65\textwidth]{Plots/ToyLinear4DVar_PostCov}
        \caption{Entries of the posterior covariance matrix and the associated errors.
          Left: the posterior covariance matrix generated by snippet~\ref{snipp:toy_linear_posterior_covariance}.
          Middle: the closed-form posterior covariance matrix given by~\eqref{eqn:Gaussian_Posterior_Toy}.
          Right: RMSE obtained by pointwise comparison of the covariance matrices obtained by solving
          the inverse problem (left) and by using the closed form (middle).}
      \label{fig:toy_linear_posterior_covariance}
    \end{figure}
    %

    \paragraph{An OED experiment}

    The simulation model is instantiated with a model grid of size \code{nx=5}, and the
    observation operator copies the model state.
    In fact, one can inspect the model array representation $\mat{A}$ for this toy linear model
    by calling \code{model.get\_model\_array()}.
    Both the model state and the observation vector sizes here are $5$. Thus, there are actually $5$ candidate sensor locations
    (observation gridpoints), and one can try to find the optimal subset of sensors by using an OED implementation.

    Here we briefly illustrate utilizing an OED object to find the A-optimal design for the toy linear example discussed above; see
    snippet~\ref{snipp:toy_linear_A_oed}.
    First, the proper OED module (here following~\cite{attia2022optimal}) is imported and is used to create the
    \code{oed\_problem} instance.
    The A-optimality criterion is registered (which can be changed later by registering a proper OED criterion), and the OED problem is then solved.
    %
    The results of solving the OED problem, for example, the optimal observational design, are then stored in
    \code{oed\_result}, which is an instance of (or derived from) the \code{pyoed.oed.OEDResults} class. This class
    provides access to various attributes of the OED problem and the solution process (such as the optimization trajectory and
    brute force solution if requested).
    This gives a taste of how simple it is to create and test DA and OED problems in {PyOED}\xspace.
    Further details on OED implementations in {PyOED}\xspace are discussed in the following
    section (\ref{subsec:results:standard_oed_examples}).

    \lstinputlisting[language=Python,
        label={snipp:toy_linear_A_oed},
        firstline=103,lastline=106,
        caption={Create an OED object, and solve the OED optimization problem for the toy linear model.}]{./Code_Snippets/toy_linear.py}

  \subsection{A standard model-constrained OED experiment}
  \label{subsec:results:standard_oed_examples}

    Parameter identification for an advection-diffusion (AD) model is the foundation of
    an experiment widely used in the model-constrained OED literature for validating
    theoretical developments; see, for example,~\cite{huan2013simulation,bui2013computational,alexanderian2016fast,attia2018goal}.
    Comparing independent scientific OED developments is admittedly hard, mainly because of the lack of availability of
    open software packages developed for OED. This is one of the main goals and features of {PyOED}\xspace.
    %
    Specifically, {PyOED}\xspace will enable OED researchers to compare the performance of new OED algorithmic approaches with other methods.
    Moreover, it enables comparison with solution by brute force search for small- to moderate-dimensional problems.
    %
    In this section~\ref{subsec:results:standard_oed_examples} we describe in detail the steps required to construct and solve
    an OED problem in {PyOED}\xspace with an AD simulation model.
    This problem has been utilized independently in several OED developments;
    see, for example,~\cite{bui2013computational,alexanderian2016fast,attia2018goal,attia2022optimal,attia2022stochastic}.
    Here, we show how {PyOED}\xspace can be used to solve and benchmark this OED problem, thus providing a starting point for utilizing and developing
    multiple approaches for solving OED problems in general in {PyOED}\xspace.
    Following the same approach as in~\ref{subsec:ToyLinear},
    we start by describing the components of the inverse problem and briefly show how they are initialized in {PyOED}\xspace;
    then we initialize and solve the OED problem using the efficient stochastic approach summarized by Algorithm~\ref{alg:REINFORCE_baseline}.
    Additionally, we discuss the steps that should be modified to utilize other solution formulations
    and methods such as the relaxation approach~\eqref{eqn:relaxed_oed}.

    The code summarized here is provided in the \code{pyoed.examples.OED\_AD\_FE} module with additional comments,
    details, and capabilities.
    A Jupyter Notebook \code{pyoed.tutorials.OED\_AD\_FE} is also available and can be used to regenerate the numerical
    results presented in this section (\ref{subsec:ToyLinear}).

    \paragraph{The simulation model}

      The governing equation of the contaminant field
      $u = \xcont(\mathbf{x}, t)$ is modeled by the following AD model equations
      with the associated boundary conditions:
      %
      \begin{equation}\label{eqn:advection_diffusion}
        \begin{aligned}
          \xcont_t - \kappa \Delta \xcont + \vec{v} \cdot \nabla \xcont &= 0
            \quad \text{in } \domain \times [0,T],   \\
          \xcont(x,\,0) &= \theta \quad \text{in } \domain,  \\
          \kappa \nabla \xcont \cdot \vec{n} &= 0
            \quad \text{on } \partial \domain \times [0,T],
        \end{aligned}
      \end{equation}
      %
      where $\kappa>0$ is the diffusivity, $T$ is the simulation final time and
      $\vec{v}$ is the velocity field.

      The domain is $\domain:=(0, 1) \times (0, 1)$ with two rectangular regions
      modeling two buildings inside the domain.
      The velocity field $\vec{v}$ is known exactly and is obtained by solving a steady Navier--Stokes
      equation, with the side walls driving the flow, as detailed in~\cite{PetraStadler11,attia2022optimal}.

      To create a simulation model object implementing~\eqref{eqn:advection_diffusion}, ground truth of the initial condition,
      and plot the domain (with finite elements discretization) as well as the velocity field, one can use Snippet~\ref{snipp:AD_FE_Model};
      the output is shown in Figure~\ref{fig:AD_domain_velocity}.

      \lstinputlisting[language=Python,
        label={snipp:AD_FE_Model},
        firstline=9,lastline=19,
        caption={Create an object representing the simulation model~\eqref{eqn:advection_diffusion}.}]{./Code_Snippets/OED_AD_FE.py}

      %
      \begin{figure}[!htbp]
        \centering
          \includegraphics[width=0.30\textwidth]{Plots/AD_Domain}
          \qquad
          \includegraphics[width=0.30\textwidth]{Plots/AD_VelocityField}
          \caption{
            Left: finite elements discretization of the domain $\domain$ of the AD problem~\eqref{eqn:advection_diffusion}.
            Right: the velocity field $\vec{v}$.
          }
          \label{fig:AD_domain_velocity}
      \end{figure}
      %

    \paragraph{The prior}

      In this setup, following~\cite{attia2018goal,PetraStadler11,VillaPetraGhattas2016},
      we choose a Laplacian prior of the parameter $\iparam$ is $\GM{\iparb}{\Cparampriormat}$,
      with $\Cparampriormat$ being a discretization of $\mathcal{A}^{-2}$,
      where $\mathcal{A}$ is a Laplacian operator.
      In {PyOED}\xspace, a Laplacian prior can be created as described in snippet~\ref{snipp:Laplacian_prior}.

      \lstinputlisting[language=Python,
        label={snipp:Laplacian_prior},
        firstline=22,lastline=28,
        caption={Create a Laplacian prior.}]{./Code_Snippets/OED_AD_FE.py}

    \paragraph{The observation operator}

      A common observational configuration is to consider uniformly distributed candidate sensor locations
      and solve an OED problem to choose the optimal subset of candidate sensor locations.
      A uniform observation operator can be created and incorporated in this problem as described in
      snippet~\ref{snipp:AD_uniform_observations}.

      \lstinputlisting[language=Python,
        label={snipp:AD_uniform_observations},
        firstline=32,lastline=36,
        caption={Create a uniform observation operator with $10$ candidate locations.}]{./Code_Snippets/OED_AD_FE.py}

      Assuming Gaussian observational noise model, a Gaussian observation error model is created as described
      in~\ref{snipp:AD_Gaussian_Noise}

      \lstinputlisting[language=Python,
        label={snipp:AD_Gaussian_Noise},
        firstline=40,lastline=43,
        caption={Create a Gaussian noise model.}]{./Code_Snippets/OED_AD_FE.py}

    \paragraph{The inverse problem: 4DVar}

      As with the case of the toy linear problem described above in Section~\ref{subsec:ToyLinear}, the elements of the inverse problem here can be created
      as described by snippet~\ref{snipp:AD_inverse_problem}.
      Similarly, synthetic observations (data) can be created  as in snippet~\ref{snipp:AD_synthetic_data}.
      Note that all steps followed so far in this example are  similar to those followed in the case
      of the toy linear model discussed in Section~\ref{subsec:ToyLinear}.

      \lstinputlisting[language=Python,
        label={snipp:AD_inverse_problem},
        firstline=46,lastline=54,
        caption={Create the DA object.}]{./Code_Snippets/OED_AD_FE.py}

      \lstinputlisting[language=Python,
        label={snipp:AD_synthetic_data},
        firstline=57,lastline=66,
        caption={Create synthetic observations, and associate them to the inverse problem object.}]{./Code_Snippets/OED_AD_FE.py}

    \paragraph{The OED problem}

      The discussion above is valid for all model-constrained OED approaches in {PyOED}\xspace. In what follows, we describe the steps needed
      to create an OED object that follows the stochastic approach~\eqref{eqn:stochastic_OED_optimization}.
      An OED object is created in {PyOED}\xspace by following the steps in snippet~\ref{snipp:AD_OED_object}.
      Here, we seek to activate only $4$ sensors out of the candidate $10$, and we use the trace of the Fisher information matrix ($\mathsf{FIM}$)
      as the utility function. To enforce the budget, we use an $\ell_0$ penalty term as detailed in~\cite{attia2022stochastic}.

      \lstinputlisting[language=Python,
        label={snipp:AD_OED_object},
        firstline=69,lastline=81,
        caption={Create an OED object implementing the stochastic approach~\ref{alg:REINFORCE_baseline}.}]{./Code_Snippets/OED_AD_FE.py}

      The OED problem can be solved as described by snippet~\ref{snipp:AD_OED_solve}.

      \lstinputlisting[language=Python,
        label={snipp:AD_OED_solve},
        firstline=85,lastline=90,
        caption={Solve the OED problem.}]{./Code_Snippets/OED_AD_FE.py}

     The resulting \code{oed\_results} object can be used to generate several analysis plots as
     described by the code in snippet~\ref{snipp:AD_OED_plot}.
     This generates multiple standard plots, including the
     performance of the optimization algorithm over consecutive iterations in Figure~\ref{fig:AD_FE_Plots}(left),
     the optimal sensor locations generated by the optimization algorithm in Figure~\ref{fig:AD_FE_Plots}(middle),
     and comparison of the quality of the solution with respect to brute force search shown in Figure~\ref{fig:AD_FE_Plots}(right).

      \lstinputlisting[language=Python,
        label={snipp:AD_OED_plot},
        firstline=92,lastline=93,
        caption={Create standard plots for assessing the performance of the optimization routine
          and the quality of the generated design.}]{./Code_Snippets/OED_AD_FE.py}

     %
     \begin{figure}[!htbp]
       \centering
         \includegraphics[width=0.35\textwidth]{Plots/AD_OED_optimization_iterations}
         \quad
         \includegraphics[width=0.22\textwidth]{Plots/AD_OED_optimal_sensors}
         \quad
         \includegraphics[width=0.35\textwidth]{Plots/AD_OED_bruteforce}
       \caption{Subset of the plots generated by running the code in snippet~\ref{snipp:AD_OED_plot}.
         Left: value of the utility (objective) function, i.e., the penalized OED criterion, over consecutive iterations of the optimization algorithm.
         Middle: optimal solution, showing optimal sensor locations in the domain.
         Right: value of the objective of the optimal solution (red star) returned by algorithm~\ref{alg:REINFORCE_baseline},
         compared with the global optimum solution (black $x$ mark),
         and all possible solutions marked as blue circles; the x-axis shows the indexes of all possible binary designs from $1$ to $2^{\Nsens=10}=1024$,
         and the y-axis shows the corresponding values of the optimization objective.
       }
       \label{fig:AD_FE_Plots}
     \end{figure}
     %

     To use other OED formulations to solve the same problem, the user only needs to update the code in the snippet~\ref{snipp:AD_OED_object}
     with the proper OED implementation.
     For example, the relaxation approach~\eqref{eqn:relaxed_oed} can be used as illustrated in the case of the linear toy model above
     in~\ref{subsec:ToyLinear}; see snippet~\ref{snipp:toy_linear_A_oed}.
     %
     Specifically, the relaxation approach~\eqref{eqn:relaxed_oed} can be used to solve the present optimal sensor placement problem by replacing the code in
     snippet~\ref{snipp:AD_OED_object} with the following code in snippet~\ref{snipp:AD_OED_relaxed_object}, which demonstrates the simplicity of {PyOED}\xspace interface.
     Results of snippet~\ref{snipp:AD_OED_relaxed_object} are omitted from the presentation here for clarity and because the main goal here is to discuss
     usage of the approaches in {PyOED}\xspace rather than assessing the quality of the solution approach, which is left for interested users of the package and for future benchmarking research.

      \lstinputlisting[language=Python,
        label={snipp:AD_OED_relaxed_object},
        firstline=96,lastline=110,
        caption={Create an OED object implementing the relaxation approach~\eqref{eqn:relaxed_oed}.}]{./Code_Snippets/OED_AD_FE.py}

     Note, however, that we had to change the penalty function in snippet~\ref{snipp:AD_OED_relaxed_object} because
     the $\ell_0$ penalty function used in snippet~\ref{snipp:AD_OED_object}
     is non differentiable, while the relaxation approach requires the OED objective function to be differentiable.
     For details, see, for example,~\cite{attia2022stochastic}.

\section{Concluding Remarks}
\label{sec:Conclusions}
  This work describes {PyOED}\xspace, a highly extensible high-level software package for OED in inverse problems and DA.
  {PyOED}\xspace aims to be a comprehensive Python toolkit for model-constrained OED.
  The package targets scientists and researchers interested in understanding the details of OED formulations and approaches.
  It is also meant to enable researchers to experiment with standard and innovative OED technologies within external test problems (e.g., simulations).
  %
  The mathematical formulations of OED, inverse problems, and DA overlap significantly, and thus, we plan to extend {PyOED}\xspace with a plethora of Bayesian inversion, DA, and OED implementations as well as new scientific simulation models,
  observation error models, and observation operators.
  %

  While we focused the discussions in this paper on specific OED approaches, the current version {PyOED}\xspace (v1.0) provides several other implementations and emphasizes
  implementing the essential infrastructure that enables combininig DA and OED elements with other parts of the package.
  The main limitation of the initial version of {PyOED}\xspace is scalability.
  Specifically, the concept is developed without parallelization capability.
  In future versions of {PyOED}\xspace, scalability will be achieved by adding message passing interface (MPI) support, for example using the \code{mpi4py} package,
  and by supporting \code{PETSc}~\cite{balay2022petsc}.
  Performance will also be enhanced by converting or rewriting suitable parts of the package in \code{Cython}.

redpajama_set_name: RedPajamaArXiv


Okay I just got around to finally downloading the 401 patch today.
I'm a graphics artist in the industry, and that actually means I get to spend LESS time playing games!! Go figure.
Anyways, I have to say, wow, its soo much fun!
I'm flying with the royal navy and man has this gotten wack, it took me like 10 trys to get of the carrier once lol.
dogfights are insanly fun to! thankfully my fighter has a ton of guns.. but ammo seems to be lacking lol.
I'm not exactly a "ww2 plane expert" but this game seems to have insane detail.
My first dogfight rocked, after finally getting off the carrier I headed in full bore against the enemy, thankfully this map seems to avoid the massive trips of the others and I was in the action in no time, a few head on passes and some turns and I got my first kill!!
then I heard a THUNK THUNK THUNK!!
Before I knew it smoke was everywhere and my plane was shaking like hell, I switched to the outside few and I could see parts falling off and some guy on my *** blowing it apart lol!!
I was also just comming up to the crest in my climb, the plane stalled out and ahhhhh soooo gracefully fell back to the earth it was awsome... then..
then wow, just like in those ww2 movies it started doing an insane nose dive!! I couldnt even pull out in time so I hit the silk!!
Great work guys! this rocks!
Now I'm off to see if I can take off again!
Yep. I can stay away for a little while, but once I fly again, I'm in love with the sim all over again.
When I first got the game, this is exactly how I was like. Its awesome and the best way to approach your first outings.
Glad you are enjoying yourself..... welcome aboard.
the FM improvement is staggering!
cant wait to get the full merged install.


today was one of those lazy saturdays - the sun was hiding + the gentle pitter patter of rain accompanied the afternoon. this past week has been pretty crazy, so i took the liberty of lounging in my bed, enjoying of voice of Lewis Watson [ so. good. ] + finally getting some writing done.


text:
New England Mutiny get back in win column with victory over Worcester
Posted Jun 4, 2018
Sydney Brackett notched two goals and an assist in New England's win Saturday. (George Skovera)
By Conor Ryan | cryan@masslive.com
After dropping their first loss of the season to the expansion Worcester Smiles last weekend on the road, the New England Mutiny exacted revenge Saturday night at Agawam High School.
New England captain Kate Howarth scored in the fifth minute of play and the Mutiny never looked back -- as they cruised to a comfortable 4-0 win over the Smiles in the team's first home game at Agawam H.S. since the 2011 season.
Syracuse University product Sydney Brackett notched two goals and an assist, including her second tally just two minutes into the second half, to ice the result. Howarth scored the final goal in the victory while midfielder Brenna Brown was credited with two assists.
New England goalkeeper Lysianne Proulx continued her rock solid play in net, making seven saves.
With a 2-1-1 record, the Mutiny are now in second place behind conference leader Connecticut Fusion (3-1) in the UWS East race. The top two teams in the East will battle it out for the spot in the four-team UWS Championship (July 20-22). New England is back in action at Agawam H.S. on Saturday, June 9th (5:30 PM kickoff) when they take on the New York Surf. The two teams split their series last season as both won their respective home games. This will be the only match between them in 2018.

redpajama_set_name: RedPajamaCommonCrawl


For the first time, Golden Guitar winners Travis Collins and Amber Lawrence have collaborated on their debut EP called Our Backyard which will be released on Friday August 4, and to celebrate they will be heading out on tour from August 3 with tickets on sale now.
Brimming with classics that tell quintessentially Australian stories, the music and lyrics on the Matt Fell produced Our Backyard reflect the true Australian experience across life, love and friendship.
With both artists coming off album releases, Amber and Travis agree that now was the best time to combine their talents and dedicate their efforts to creating a collaborative and genuinely enjoyable album. "Now that our solo careers are cemented and we've both come off our own solo albums we can have a little fun and throw caution to the wind," says Amber.
Friends since 2004 when they met at the College of Country Music (now Academy of Country Music), Amber and Travis wanted to create music that was more than simple duets but a true collaboration containing the essence of both artists; two distinct voices coming together.
Prior to this release, Amber has been on the road promoting her fifth studio album Happy Ever After which came out in September 2016. The two singles 'Happy Ever After' and 'Cheers to the Girls' were snapped up by country radio, hitting the top chart radio and video playlists. One of the hardest working artists in the industry, the 2015 & 2016 CMC Oz Female Artist of the Year and 2015 Golden Guitar Female Artist of the Year continues creating authentic and progressive country music that engages her fans like no other.
Travis Collins, a mainstay in the Australian music landscape, has cemented his footprint as one of Australia's best-loved country music artists. Earlier this year he took home not only his first Golden Guitar but a total of three: Male Artist of the Year, Single of the Year (Just Another Girl) and APRA AMCOS Song of the Year (Call Me Crazy). Two months later, Travis was also awarded three CMC Awards: Australian Video of the Year, Male Artist of the Year, and CMC Oz Artist of the Year. As well-known for his masterful guitar playing and powerhouse vocals as he is for his depth of song-writing and high energy live performances, Travis continues to go from strength to strength with his most recent album Hard Light producing three #1 songs and charting in the ARIA top 20 all-genre charts.