\documentclass[a4paper,12pt]{article}
%\newcommand{\figs}{.}
%\usepackage[eng]

\title{The COSTA Workbench}
\author{Julius Harry Sumihar and Nils van Velzen}
%\memonum{CTA memo200802}
\date{\today}

\begin{document}

%\memotitlepage

%\begin{vtlogsheet}
%\vtlogentry{1.0}{JS}{2008-01-18}{Initial version}{CvV}
%\vtfilelocation{$<$COSTA\_DIR$>$/doc/costawb}
%\end{vtlogsheet}

\tableofcontents

%===============================================================================

\begin{abstract}
COSTA is a generic problem solving environment for data assimilation and
model calibration. It aims to provide, among others, reduction of software
development cost, re-usability of methods and models, and ability to quickly
try out alternative methods. Although COSTA is still an ongoing project, it
is important to evaluate the existing version of COSTA. This study aims to
review COSTA by implementing a number of algorithms. Three parameter
estimation algorithms and two data assimilation algorithms are implemented
in this study. These algorithms are applied to a number of models, which
already comply the COSTA environment. The experience in implementing these
algorithms show that once a method is ready, it is very easy to apply it to
different models. Moreover, the COSTA functions reduce the number of codes
lines used in developing the programs and make the programs more readable.
\end{abstract}

\section{Input of the COSTA Workbench}
The COSTA Workbench brings together the three main parts in data
assimilation and model calibration:
\begin{itemize}
\item Observations
\item Model
\item Data Assimilation / Calibration method
\end{itemize}

The main structure of the XML-input file of the workbench defines these
three parts. 




The following example shows the top level of an input file:

\begin{verbatim}
<?xml version="1.0" encoding="UTF-8"?>
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
   <!-- Observations used for data assimilation or calibration -->
   <CTA_SOBS id="obs_assim" database="obs_lorenz96.sql"/>
   
   <!-- Observations used for validation  -->
   <CTA_SOBS id="obs_valid" database="obs_lorenz96.sql"/>
   
   <!-- Model that is used -->
   <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />

   <!-- Data assimilation or calibration method that is used -->
   <method name="ensemble">
      <!-- Configuration of the method -->
  </ensemble_filter>
  </method>
</costa>
\end{verbatim}

The section {\tt obs\_assim} specifies the observations
that are used in the calibration process or are assimilated in the data
assimilation method. The section {\tt obs\_valid} is optional and not
supported by all data assimilation methods. This will specify the
observations that are compared to the model predictions but these are only
used for output purposes and do not influence the simulation/calibration
results.

The section {\tt modelclass} specifies which model is used. In this example
we use a model that is generated by the SP-modelbuilder.

The section {\tt method} specifies the data assimilation or calibration method
that is used. The tag name specifies the used method. The available methods
are: {\tt ensemble}, {\tt rrsqrt}, {\tt coffee}, {\tt ensrf},  
{\tt lbfgs}, {\tt conjugrad} and {\tt simplex}. The configuration of the
method depends on the data assimilation and calibration method.

In the workbench, all sensible combinations of models and methods are
available. Each combination is completely described using an XML-file of the
name {\tt <method>}\_{\tt <model>.xml}.

Thus, to run the workbench for a certain application, use \vspace{1cm}

{\tt \% costawb < xmlfile > } 
\vspace{1cm}

in the example above: 

{\tt \% costawb ens\_lorenz96.xml} 

 

 In this document, we will briefly describe the
models and then  explain the theory and practical use of each method.

   \section{Models used in the workbench}
   In the COSTA workbench, we use five models which are already available
   within the present COSTA. These are small models for testing and teaching
   purposes. All models are represented as a set of differential equations
   where the solution is obtained numerically by using Runge-Kutta method or
   Forward Euler.

The models are the following:
      \subsection{Oscillator model}
      The oscillator model is a simple mass-spring model with friction. It has two describing variables, which are the location of the mass $x$ and its velocity $u$. The two variables are related according to the following equations:
      \begin{eqnarray}
         \frac{dx}{dt}=u \\
         \frac{du}{dt}=-\omega^2 x - \frac{2}{T_{d}} u
      \end{eqnarray}
      where $\omega$ is the oscillation frequency, which depends on the mass and the spring constant, while $T_{d}$ is the damping time.

     \subsection{Lorenz model}
     Edward Lorenz (1963) developed a very simplified model of convection called the Lorenz model. The Lorenz model is defined by three differential equations giving the time evolution of the variables $x,y,z$:
      \begin{eqnarray}
         \frac{dx}{dt}=\sigma(y-x) \\
         \frac{dy}{dt}=\rho x - y -x z \\
         \frac{dz}{dt}=x y - \beta z
      \end{eqnarray}
      where $\sigma$ is the ratio of the kinematic viscosity divided by the thermal diffusivity, $\rho$ the measure of stability, and $\beta$ a parameter which depends on the wave number.

      This model, although simple, is very nonlinear and has a chaotic nature.
      Its solution is very sensitive to the parameters and the initial
      conditions: a small difference in those values can lead to a very different solution.

      \subsection{Lorenz96 model} %source: Miyoshi,T.[2004]
      The Lorenz96 (see Lorenz and Emanuel, 1998) is defined by the following equation
      \begin{equation}
         \frac{dx_i}{dt}=x_{i-1}(x_{i+1}-x_{i-2})-x_{i}+F
      \end{equation}
      where $i=1,...,N$, $N=40$ and the boundary is cyclic, i.e.
      $x_{-1}=x_{N-1}$, $x_{o}=x_{N}$, and $x_{N+1}=x_{1}$, and $F=8.0$. The
      first term of the right hand side simulates ``advection'', and this model
      can be regarded as the time evolution of an arbitrary one-dimensional
      quantity of a constant latitude circle; that is, the subscript $i$
      corresponds to longitude. This model also behaves chaotically in the case
      of external forcing $F=8.0$.

      \subsection{Heat transfer model} %source: wikipedia http://en.wikipedia.org/wiki/Heat_equation
      The model represents a special case of heat propagation in an isotropic and homogeneous medium in the 2-dimensional space. The equation can be written as follows:
      \begin{equation}
         \frac{\partial T}{\partial t}=k(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2})
      \end{equation}
      for $x \in [0,X]$ and $y \in [0,Y]$ where $T$ is the temperature as a
      function of time and space and $k$ is a material-specific quantity
      depending on the thermal conductivity, the density and the heat capacity.
      Here $k$ is set to 1. Neumann and Dirichlet boundary conditions  are
      used.

      \subsection{1-dimensional Advection model}
      %source: Martin's code; Chapter 11 Ad-Diff equations and turbulence; Grima and Newman [2004]
      In this study, a 1-dimensional advection model is also used and can be written as follows
      \begin{eqnarray}
         \frac{\partial c}{\partial t}=u \frac{\partial c}{\partial x}
      \end{eqnarray}
      where $c$ typically describes the density of the particle being studied
      and $u$ is the velocity. On the left boundary $c$ is specified as
      $c_b(t)=1+sin(\frac{2\pi}{10}t)$.

      \subsection{Stochastic Extension}
      The previous subsections mention the deterministic models used in this
      study. Especially for the implementation of Kalman filtering, we need to
      extend the models into a stochastic environment. This is done, for the
      oscillation, Lorenz, and Lorenz96 models, by adding a white noise process
      to each variable. On the other hand, for the heat transfer and 1-d
      advection models, this is done by adding a colored noise process,
      represented by an AR(1) process, to the boundary condition. For the heat
      model the noise is also spatially correlated.




\section{Calibration methods used in the workbench}
As already mentioned, three parameter estimation algorithms and four sequential data
assimilation methods have been implemented and coupled with the models described
earlier. We first  discuss the three parameter estimation algorithms.

In the parameter estimation algorithms, the basic idea is to find
the set of model parameters which minimizes the cost function measuring the
distance between the observation and the model prediction. Two different
cost functions are implemented. The first one is similar to equation
(\ref{eq.costfunction1}), while in the second one we add the background
component as follows
\begin{equation}
\label{eq.costfunction2}
  J(x_o) = (x_o-x^b)'(P^b)^{-1}(x_o-x^b)+\sum_{k=1}^N (y^o(k)-Hx(k))R^{-1}(y^o(k)-Hx(k))'
\end{equation}
Where $x^b$ and $P^b$ are the background or initial estimate of $x_o$ and its covariance respectively. This additional component ensures that the solution will not be too far from the initial guess.

\subsection{Parameter estimation with the Simplex method}
\label{sssec.simplex}
The simplex method that is implemented in COSTA is the one due to Nelder and
Mead (1965). It is a systematic procedure for generating and testing
candidate vertex solutions to a minimization problem. It begins at an
arbitrary corner of the solution set. At each iteration, the simplex method
selects the variable that will produce the largest change towards the
minimum solution. That variable replaces one of its compatriots that is
most severely restricting it, thus moving the simplex to a different corner
of the solution set and closer to the final solution.

A simplex is the geometrical figure consisting, in \emph{N} dimensions, of
$N+1$ points and all their interconnecting line segments, polygonal faces,
etc. In two dimensions, a simplex is a triangle. In three dimensions it is
a tetrahedron. The simplex method must be started with $N+1$ points of
initial guess, defining the initial simplex. The simplex method now takes a
series of steps. The first step is to move the vertex where the cost is
largest through the opposite face of simplex to a lower point. This step is
called the \emph{reflect} step. When the cost of the new vertex is even
smaller than all the remaining, the method expands the simplex even
further, called the \emph{expand} step. If none of these steps produce a 
better vertex, the method will contract the simplex in the same direction
of the previous step and take a new point. This step is called the
\emph{contract} step. When the cost of the new point is still not better
than the previous one, the method will take the last step called
\emph{shrink}. In this step all of the simplex points, except the one with
the lowest cost, are 'shrinked' toward the best vertex.

As it basically only tries and compares the solution of several different
sets of parameters, the method requires only function evaluations and not
derivatives.

\subsubsection{Using the Simplex method}
For using the implemented simplex method the user needs to specify in the
input file the \emph{initial-guess} of the set of parameters to calibrate
as well as the \emph{initial-step} for creating the initial simplex. The
initial-step consists of $N$ entries, where $N$ is the number of parameters
to calibrate. An initial vertex is obtained by adding an element of the
initial-step to the corresponding parameter in the initial-guess.
Performing this one by one to all the parameters in the initial-guess,
there are $N+1$ vertices forming the initial simplex. Although it can be
easily extended, at the moment the program only supports the model with
a maximum of 10 parameters.

The stopping criteria for the iteration are the maximum number of iteration
and the maximum cost difference between the worst and the best vertices,
where worst vertex refers to the one with the biggest cost value and best
vertex is the one with the lowest. When the maximum cost difference is very
small we may expect that the (local) minimum of the cost function has been
reached. The optimum solution is the vertex with the lowest cost. However,
as output we display all the final vertices with their respective costs.

Since the model parameters are stored as COSTA state vector, most variable
operations are performed in term of COSTA state vector. The operations are
usually to compute new vertex. This can easily be carried out by using the
combination of functions like cta\_state\_axpy and cta\_state\_scal, with
the aid of functions like cta\_state\_duplicate and cta\_state\_copy for
assigning values to working variables. These COSTA functions reduce
significantly the lines of codes required to perform the operation.

An important function required to implement the simplex method is the sorting
function. This function is used to sort the vertices with respect to their
cost values in descending order. While the vertices are stored as COSTA
state vectors, their cost values are stored as a Fortran array. The sorting
is done by using an external routine, which works only with Fortran
variables. The output of this routine are the sorted array of cost values
and an integer array containing the index of the sorted array.

In the section above  we learnt that the simplex method consists
of several steps which are taken to find a new vertex with lower cost. When
all the steps do not produce a better vertex the program will give a
warning. For some cases this may indicate that the solution does not exist.
This occurs for example with the Lorenz model. Since it is a chaotic model,
small difference in the parameters yields very different cost values. This
makes the Lorenz model not a very good model for calibration tests in the
present setup. Perhaps shortening the interval over which the optimization
is carried out can make the cost function better behaved.
\subsubsection{The configuration of the Simplex method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt iteration}: Number of iterations and convergence tolerance
\item {\tt model<n>}: Model configuration for the initial model state at every
starting vertex
\end{itemize}

\subsubsection{XML-example}

\begin{verbatim}
  <parameter_calibration>
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="60" tol="0.001">    </iteration>
    <!-- 1st VERTEX:-->
    <model1>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="8.95" avg_omega="13.5"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model1>
    <!-- 2nd VERTEX:-->
    <model2>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="6.0" avg_omega="15.9"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model2>
    <!-- 3rd VERTEX:-->
    <model3>
       <modelbuild_sp>
         <xi:include href="models/functions_oscill.xml"/>
          <model>  
             <parameters avg_t_damp="9.0" avg_omega="12.5"> </parameters> 
          </model> -->  
       </modelbuild_sp> 
    </model3>
  </parameter_calibration>
\end{verbatim}
 %\caption{Example of input file for simplex method\label{Fig:Simplex_input}}

\subsection{Parameter estimation with the Conjugate Gradients method}
\label{sssec.conjugrad}
The problem of minimization of multivariable function is usually solved by
determining a \emph{search direction} vector and solve it as a line
minimization problem. If $x$ is a vector containing the variables to be
determined and $h$ is the vector of search direction, at each iteration
step the minimization problem of a function $f$ is formulated as to find
the step size $\lambda$ that minimizes $f(x+\lambda h)$. At the next
iteration, $x$ is replaced by $x+\lambda h$ and a new search direction is
determined. Different methods basically propose different ways of finding
the search direction.

The conjugate gradient method is an algorithm for finding the nearest local
minimum of a function which uses \emph{conjugate directions} for going
downhill. Two vectors $u$ and $v$ are said to be conjugate (with respect to
a matrix $A$) if
\begin{equation}
   u'Av=0
\end{equation}
where in the minimization problem, $A$ is typically the \emph{Hessian}
matrix of the cost function. In the conjugate gradient methods, the search
direction is somehow constructed to be conjugate to the old gradient.

The two most important conjugate gradient methods are the
\emph{Fletcher-Reeves} and the \emph{Polak-Ribierre} methods (see Press
et.al., 1989). These algorithms calculate the mutually conjugate directions
of search with respect to the Hessian matrix of the cost function directly
from the function and the gradient evaluations, but without the direct
evaluation of the Hessian matrix. The new search direction $h_{i+1}$ is
determined by using
\begin{equation}
   h_{i+1} = g_{i+1}+\gamma_i h_i
\end{equation}
where $h_i$ is the previous search direction, $g_{i+1}$ is the negative of
local gradient at iteration step $i+1$, while $\gamma_i$ is determined by
using the following equations for \emph{Fletcher-Reeves} and the
\emph{Polak-Ribierre} methods respectively:

\begin{eqnarray}
   \gamma_i=\frac{g_{i+1} \cdot g_{i+1}}{g_{i} \cdot g_{i}} \\
   \gamma_i=\frac{(g_{i+1}-g_{i}) \cdot g_{i+1}}{g_{i} \cdot g_{i}}
\end{eqnarray}

If the vicinity of the minimum has the shape of a long, narrow valley, the
minimum is reached in far fewer steps than would be the case using the
\emph{steepest descent} method, which makes use of minus of the local
gradient as the search direction.

In this study, the line minimization to find the step size $\lambda$ that
minimizes $f(x+\lambda h)$ at every iteration step is done by using the
\emph{golden section search} algorithm. This is an elegant and robust
method of locating a minimum of a line function by bracketing it with three
points: if we can find three points $a,b,$ and $c$ where $f(a)>f(b)<f(c)$
then there must exist at least one minimum point in the interval $(a,c)$.
The points $a,b,$ and $c$ are said to \emph{bracket} the minimum. This
algorithm involves evaluating the function at some point $x$ in the larger
of the two intervals $(a,b)$ or $(b,c)$. If $f(x)<f(b)$ then $x$ replaces
the midpoint $b$ and $b$ becomes an end point. If $f(x)>f(b)$ then $b$
remains the midpoint with $x$ replacing one of the end points. Either way
the width of the bracketing interval will reduce and the position of the
minimum will be better defined. The procedure is then repeated until the
width achieves a desired tolerance. It can be shown that if the new test
point, $x$, is chosen to be a proportion $(3-\sqrt{5})/2$ (hence Golden
Section) along the larger sub-interval, measured from the mid-point $b$,
then the width of the full interval $(a,c)$ will reduce at an optimal rate.

In the absence of an adjoint formalism in COSTA a numerical approximation
was used to compute the gradient of the cost function.


\subsubsection{Using the Conjugate gradient method}
In this COSTA-method both the Fletcher-Reeves and Polak-Ribiere methods are
implemented. For comparison purpose, we also implement a steepest-descent
method. The user can choose the method to use by specifying the field
method \emph{id} in the input file, where 1 refers to Fletcher-Reeves, 2 to
Polak-Ribiere, and 3 to steepest-descent.

The stopping criteria are the maximum number of iterations, the tolerance
for the step size $\lambda$ as described in subsection
\ref{sssec.conjugrad} and the tolerance for the local gradient. When the
line-minimization does not move the parameters significantly to the new
ones, the cost function may have reached its (local) minimum. On the other
hand, the gradient around the (local) minimum is also close to zero. Hence
it is possible to use the norm of the gradient vector as one of the
stopping criteria. The same tolerance is also used within the golden
section search routine.

Besides specifying the method to use and the stopping criteria parameters,
the user also needs to specify the initial guess of the parameters to
calibrate as well as the initial two points $AX$ and $BX$ for the golden
section search minimization.

Since the method requires the computation of the gradient, while there are no
adjoint models available, the gradient of the cost function is computed
with finite difference. This is fine for small number of parameters like
the ones used in this study. However it will become computationally
expensive for many parameters. The size of the perturbation, $delta$, must
also be specified in the input file.

In our implementation, the gradient is also stored as COSTA state to make
it easier to assign values between the variables representing model
parameters to calibrate. The norm of the gradient is computed using the
function cta\_state\_nrm2.

Besides the state vector operations, like scaling and adding, in this
method we also need to compute the dot product between two state vectors.
This can easily be done by using the function cta\_state\_dot.


\subsubsection{The configuration of the CG method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt maxit}: the maximum number of iterations
\item {\tt tol\_step, tol\_grad}: Stopping criterium  according to the step
  size or the gradient norm
%\item {\tt output}: output specification
\item {\tt AX, BX}: Starting point for golden section search
\item {\tt delta}: Perturbation size 
\item {\tt iteration}: Number of iterations and convergence tolerance
\item {\tt model}: Model configuration for the initial model
\end{itemize}

\subsubsection {XML-example}

   \begin{verbatim}
  <parameter_calibration>
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <iteration maxit="160" tol_step="0.0001" tol_grad="0.0001">
    </iteration>
    <method id="1" AX="0.0" BX="0.02" delta="1E-8">
    <!-- id=1: Fletcher-Reeves, id=2: Polak-Ribiere, id=3: Steepest-descent-->
    </method>
    <model>
       <modelbuild_sp>
       <functions>
          <create>oscillparam_create</create>
          <covariance>oscillparam_covar</covariance>
          <getobsvals>oscillparam_obs</getobsvals>
          <compute>oscillparam_compute</compute>
       </functions>
       <model>
          <parameters avg_t_damp="8.5" avg_omega="12.9">
          </parameters>
       </model>
       </modelbuild_sp>
    </model>
    <paramstd param1="3" param2="3">
    </paramstd>
  </parameter_calibration>
   \end{verbatim}


\subsection{Parameter estimation with the LBFGS method}
For the problem of minimizing a multivariable function quasi-Newton methods
are widely employed. This methods involve the approximation of the Hessian
(or its inverse) matrix of the function. The LBFGS (Limited
memory-Broyden-Fletcher-Goldfarb-Shanno) method is basically a method to
approximate the Hessian matrix in the quasi-Newton method of optimization.
It is a variation of the standard BFGS method, which is given by (Nocedal,
1980; Byrd,et.al.,1994)
\begin{equation}
  \label{eq.LBFGS_linemin}
    x_{i+1}=x_i-\lambda_iH_ig_i, \space \space i=0,1,2,\cdots
\end{equation}
where $\lambda_i$ is a steplength, $g_i$ is the local gradient of the cost
function, and $H_i$ is the approximate inverse Hessian matrix which is
updated at every iteration by means of the formula
\begin{equation}
    H_{i+1}=V'_i H_i V_i + \rho_i s_i s'_i
\end{equation}
where
\begin{eqnarray}
    \rho_i=\frac{1}{y'_i s_i} \\
    V_i=I-\rho_i y_i s'_i
\end{eqnarray}
and
\begin{eqnarray}
    s_i=x_{i+1}-x_i \label{eq.LBFGS_s} \\
    y_i=g_{i+1}-g_i \label{eq.LBFGS_y}
\end{eqnarray}

Using this method, instead of storing the matrices $H_i$, one stores a
certain number of pairs, say $m$, of pairs $\{s_i,y_i\}$ that define them
implicitly. The product of $H_i g_i$ is obtained by performing a sequence
of inner products involving $g_i$ and the $m$ most recent vector pairs
$\{s_i,y_i\}$ to define the iteration matrix.

Like in conjugate gradient methods, the line minimization for determining
$\lambda$ in equation (\ref{eq.LBFGS_linemin}) is implemented by using the
Golden Section Search algorithm.

\subsubsection{Using the LBFGS method}
In the implementation of LBFGS method we use COSTA matrices to store the
$s_i$ and $y_i$ vectors in equation (\ref{eq.LBFGS_s}) and
(\ref{eq.LBFGS_y}). These vectors are stored as the column of the COSTA
matrices. The elimination of the oldest vectors, however, requires the
combination of a function for getting a column of a COSTA matrix and
another function for setting the column. The latter function is already
available in the present COSTA, i.e. cta\_matrix\_setcol. However, the
former one is not yet available. We resorted to worked around this by
creating a loop consisting a combination of cta\_matrix\_getval and
cta\_vector\_setvals for performing this operation. 

Vectors $s_i$ and $y_i$ are stored as a COSTA state for the same reason as in
the implementation of the conjugate gradient and simplex methods. The
operations with a COSTA matrix, however, require the variable to be stored as a 
COSTA vector. Therefore at some points in the program we were required to
use cta\_state\_getvec and cta\_state\_setvec to get the values of a COSTA
state and assign it to a COSTA vector and vice versa.

\subsubsection{The configuration of the LBFGS method}

Like in the conjugate gradient method, the implemented LBFGS method
requires maximum number of iterations, step-size tolerance, tolerance for
gradient, and the perturbation size, $delta$, for computing the gradient.
However for using this method the user also needs to specify $nstore$,
which is the number of vectors $s_i$ and $y_i$ to store.

\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt iteration}: Number of iterations and convergence tolerance, {\tt
  maxln} is the maximum number of iteration steps in the line search part.
\item {\tt model}: Model configuration for the initial model
\item {\tt method}: Method configuration.
    {\tt nstore} specifies the storage size, i.e. max number of vector s and y to store.
    {\tt c1} and {\tt c2} are the parameters used in the 1st and 2nd Wolfe
    conditions (these are used in the line search part). 
    {\tt delta} specifies the perturbation size in computing the gradient.
\end{itemize}


\subsubsection{XML-example}
\begin{verbatim}

  <parameter_calibration>
    <!-- Simulatie timespan en stapgrootte via CTA_Time -->
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="10" maxln="20" tol_step="0.0001" tol_grad="0.0002">
    </iteration>
    <method nstore="3" c1="1E-4" c2="0.5" delta="1E-8">

    </method>
       <model>
         <xi:include href="models/lorenz_sp.xml"/>
       </model>
  </parameter_calibration>


 
\end{verbatim}

%--------------------------------------------------------------------

\section{Data assimilation methods used in the workbench}


\subsection{Introduction to Data assimilation}
   The terminology of 'data assimilation' originates from the field of meteorology where in 1950's and 1960's new methods were developed to improve weather forecasts. Due to rapid development of computers in the 1950's the applications of large numerical models to weather forecasting became possible. In this early phase the initial state of the model was estimated directly from the measurement. The numerical model was then used to produce the forecast. It was soon recognized that the forecasts could be improved if the initial state was not only based on the measurements but also on the forecast produced by the previous model run. The first data assimilation method was direct insertion which was based on replacement of a model variable by its measured value. Although it was simple to implement, it suffered from the lack of smoothness in the assimilation results. Much of the research at the time was concentrated on finding a proper way to introduce the measurements without introducing oscillations.

   In recent years more and more complex methods are being used. These methods are either based on the minimization of a criterion or on statistical principles. This difference divides the data assimilation methods into two classes: the variational and the sequential methods.

   \subsubsection{Variational methods}
      Variational methods aim at adjusting a model solution to all the observations available over the assimilation period (see Talagrand, 1997). In variational methods, one first defines a scalar function which, for any model solution over the assimilation interval, measures the ``distance'' or ``misfit'' between that solution and the available observations. The so-called \emph{cost function} will typically be a sum of squared differences between the observations and the corresponding model values. One will then look for the model solution that minimizes the cost function.

      A common use of variational method in meteorology is to obtain an optimal initial condition for a model forecast. A suitable cost function for this is the following:
      \begin{equation}
      \label{eq.costfunction1}
        J(x_o) = \sum_{k=1}^N (y^o(k)-Hx(k))R^{-1}(y^o(k)-Hx(k))'
      \end{equation}
      where $x_o$ is the initial value of the variable to be determined, $x(k)$ the variable at time $t_k$, $H$ the observation operator, $y^o(k)$ the observation at time $t_k$ and $R$ the representation covariance. The optimal initial condition, $x_o$, is the one that minimizes $J$.

      If a suitable initial state $x_o$ has been obtained, the analysed states are formed with a model forecast:
      \begin{eqnarray}
          x^a_o = x_o \\
          x^a(k) = M[x^a(k-1)]
      \end{eqnarray}
      for $k=1, \cdots, N$. The final analysed state $x^a_N$ is optimised given data from spatial different locations and from different times in the interval $(t_o,t_N)$. This approach is then referred to as \emph{4D-VAR}. %Data from the period before $t_o$ is used implicitly if it was used to form the background state $x^b$.

      The minimization of the cost function is often based on quasi-Newton methods. These methods require computation of the gradient of the cost function. In most situations, it is impossible to establish explicit analytical expressions for the gradient. It is possible to numerically and approximately determine the gradient through explicit finite perturbations of the initial state. But this would be much too costly for practical implementation since it requires to compute the cost function, i.e. to effectively integrate the model over the assimilation period, as many times as there are independent components in the initial states. Therefore to compute the gradient efficiently an adjoint model should be used.

   \subsubsection{Sequential Methods}
      While variational method is based on minimization of the cost function within a time interval, sequential method assimilates the data each time the observation becomes available. In sequential method the adjusted model solution is expressed as a linear combination of the forecast state and the data elements following the equation:
      \begin{equation}
        x^a(k) = x^f(k) + K(k)(y^o(k)-Hx^f(k))
      \end{equation}
      Here $x^f(k)$ represents the forecast state, while $x^a(k)$ the analysis state, i.e. the adjusted state. The gain matrix $K(k)$ describes how elements of the state should be adjusted given the difference between the measurement and the forecast. Different methods are available to determine $K$.

      One of the popular sequential methods is the \emph{optimal interpolation} (see Daley, 1991). This method uses a gain matrix based on an empirical covariance function or matrix. The basic assumption is that both the forecast and the observation error are normally distributed with mean zero and known covariances. The idea of optimal interpolation is to set the analysed state to the conditional mean of the true state given the observations, $x^a(k)=E[x^t(k)|y^o(k)]$. Application of Bayes theorem to Gaussian distribution shows that this could be achieved with a linear gain:
      \begin{eqnarray}
          x^a(k)=x^f(k)+K(k)(y^o(k)-H'x^f(k)) \\
          K(k)=P^f(k) H [H' P^f(k) H + R(k)]^{-1} \label{eq.K_OI}
      \end{eqnarray}
      The gain matrix $K$ in equation (\ref{eq.K_OI}) is known as the \emph{conditional mean gain} or the \emph{minimal variance gain}. A problem is how to choose suitable covariance matrices $P$ and $R$. In this method the user needs to specify the error covariance at each assimilation time.

      Another development in this class is the Kalman filtering. The Kalman filter can be seen as an extension of the optimal interpolation scheme, accounting for the evolution of errors from previous times. The target of the Kalman filter is to obtain a distribution for the true state in terms of a mean $\hat{x}$ and covariance $P$, given the model and the measurements. Like in the optimal interpolation method, the Kalman filter also assumes normally distributed forecast and observation errors. The Kalman filter is originally derived for linear systems, which in state-space form can be written as:
      \begin{eqnarray}
          x(k+1)= M(k) x(k) + \eta(k) \\
          y(k) = H(k) x(k) + \nu(k)
      \end{eqnarray}
      where $x$ is the system state, $A$ the linear model operator, $\eta \sim N(0,Q)$ the system noise, $y$ the predicted observation, $H$ the observation operator, and $\nu \sim N(0,R)$ the observation error. The Kalman filter algorithm consists of two steps:
      \begin{enumerate}
       \item Forecast step:
             \begin{eqnarray}
                x^f(k+1) = M(k) x^a(k) \\
                P^f(k+1) = M(k) P^f(k) M'(k) + Q(k)
             \end{eqnarray}
       \item Analysis step:
             \begin{eqnarray}
                x^a(k)=x^f(k) + K(k) (y^o(k) - H(k) x^f(k)) \\
                P^a(k)=(I-K(k) H(k)) P^f(k) \\
                K(k) = P^f(k) H(k) (H'(k) P^f(k) H(k) + R(k))^{-1}
             \end{eqnarray}
      \end{enumerate}



\subsection{Data assimilation  with the RRSQRT method}

      The Kalman filter gives optimal estimates of $x$ and $P$ for linear models.
      The main problem of applying the Kalman filter directly to environmental
      models is the computation of the covariance matrix $P$. Since such models
      usually have a big number of states (e.g. $O(10^4)$) the covariance will
      also become very big, which causes very expensive computational costs or
      even the impossibility to compute. Another problem is that the real life model is
      usually nonlinear. Therefore methods are proposed and developed to modify
      the Kalman filter to solve these difficulties. Two most popular
      algorithms are the \emph{reduced-rank square-root} (Verlaan and Heemink,
      1995, 1997) filter and the \emph{ensemble Kalman filter} (Evensen, 1994).

      The reduced-rank square-root (RRSQRT) filter algorithm is based on a
      factorization of the covariance matrix $P$ of the state estimate
      according to $P=LL'$, where $L$ is a matrix with the $q$ leading
      eigenvectors $l_i$ (scaled by the square root of the eigenvalues),
      $i=1,...,q$, of $P$ as columns. The algorithm can be summarized as
      follows.

       \begin{enumerate}
        \item Initialization
              \begin{equation}
                 [L^a(0)]=[l_1^a(0),\cdots,l_q^a(0)]
              \end{equation}
              where $l_i^a$ denotes the $q$ leading eigenvectors of the initial covariance matrix $P_o$.

        \item Forecast step
              \begin{eqnarray}
               x^f(k)=M[x^a(k-1)] \\
               l^f_i(k)=\frac{1}{\epsilon} \{ M[x^a(k-1)+ \epsilon  l^a_i(k-1)]-M[x^a(k-1)] \} \\
               \tilde{L}^f(k)=[l_1^f(0),\cdots,l_q^f(0),Q(k-1)^{1/2}] \\
               L^f(k)=\Pi^f(k)\tilde{L}^f(k)
              \end{eqnarray}
              where $\epsilon$ represents a perturbation, often chosen close to 1, $\Pi^f(k)$ is a projection onto the $q$ leading eigenvectors of the matrix $\tilde{L}^f(k)\tilde{L}^f(k)'$. Note that here $M$, the model operator, need not be linear.

        \item Analysis step
              \begin{eqnarray}
                P^f(k)=L^f(k)L^f(k)' \\
                K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
                x^a_k = x^fa_k + K_k [y^o_k - H_k x^f_k] \\
                \tilde{L}^a(k)=\{[I - K_k H_k] L_f, K_k R_k^{1/2}\} \\
                L^a(k)=\Pi^a(k)\tilde{L}^a(k)
              \end{eqnarray}
              where $\Pi^a(k)$ is a projection onto the $q$ leading eigenvectors of the matrix $\tilde{L}^a(k)\tilde{L}^a(k)'$. This reduction step is again introduced to reduce the number of columns in $L^a_k$ to $q$ in $\tilde{L}^a(k)$.

       \end{enumerate}

%\subsubsection{Using the RRSQRT method}


%HIER MOET NOG WAT!!!

\subsubsection{The configuration of the RRSQRT method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of modes
\item {\tt mode}: Model configuration
\end{itemize}

\subsubsection{XML-example}

\begin{verbatim}
<?xml version="1.0" encoding="UTF-8"?>
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <!-- Observations used for assimilation -->
  <CTA_SOBS id="obs_assim" database="obs_lorenz.sql"/>

  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />

  <!-- Filter configuration -->
  <method name="lbfgs">
  <parameter_calibration>
    <!-- Simulatie timespan en stapgrootte via CTA_Time -->
    <CTA_TIME id="simulation_span" start="0" stop="50.0" step=".1"/>
    <output>
        <filename>results_oscill.m</filename> 
       <CTA_TIME id="times_mean" start="0" stop="50.0" step="1"/>
    </output>
    <iteration maxit="10" maxln="20" tol_step="0.0001" tol_grad="0.0002">
    </iteration>
    <method nstore="3" c1="1E-4" c2="0.5" delta="1E-8">
    <!-- nstore specifies the storage size, i.e. max number of vector s and y to store -->
    <!-- c1 and c2 are the parameters used in the 1st and 2nd Wolfe conditions-->
    <!-- delta specifies the perturbation size in computing gradient-->
    </method>
       <model>
         <xi:include href="models/lorenz_sp.xml"/>
       </model>
  </parameter_calibration>
  </method>
</costa>
\end{verbatim}


\subsection{Data assimilation  with the Ensemble method}

      While the RRSQRT represents the covariance matrix $P$ based on the first
      $q$ leading eigenvectors, the ensemble Kalman filter (EnKF) is based on a
      representation of the probability density of the state estimate by a
      finite number $N$ of randomly generated system states. The EnKF algorithm
      can be summarized as follows.

       \begin{enumerate}
        \item Initialization
              
              An ensemble of $N$ states $\xi^a_i(0)$ are generated to represent the uncertainty in $x_o$.

        \item Forecast step
              \begin{eqnarray}
               \xi_i^f(k)=M[\xi_i^a(k-1)] + \eta_i(k-1) \label{eq.EnKF_forecast} \\
               x^f(k)=\frac{1}{N} \sum_{i=1}^N \xi_i^f(k) \\
               E^f(k)=[\xi_1^f(k)-x^f(k),\cdots,\xi_N^f(k)-x^f(k)]
              \end{eqnarray}

        \item Analysis step
              \begin{eqnarray}
                P^f(k)=\frac{1}{N-1}E^f(k)E^f(k)' \\
                K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
                \xi^a_i(k) = \xi^f_i(k) + K(k) [y^o(k) - H(k) \xi^f_i(k) + \nu_i(k)]
              \end{eqnarray}
       \end{enumerate}
       where the ensemble of state vectors are generated with the realizations
       $\eta_i(k)$ and $\nu_i(k)$ of the model noise and observation noise
       processes $\eta(k)$ and $\nu(k)$, respectively.

       For most practical problems the forecast equation
       (\ref{eq.EnKF_forecast}) is computationally dominant. As a result the
       computational effort required for the EnKF is approximately $N$ model
       simulations. The standard deviation of the errors in the state estimate
       are of a statistical nature and converge very slowly with the sample
       size ($\approx N$). Here it should be noted that for many atmospheric
       data assimilation problems the analysis step is also a very time
       consuming part of the algorithm.


%\subsubsection{Using the Ensemble method}

\subsubsection{The configuration of the Ensemble method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of ensembles
\item {\tt model}: The configuration of a model
\end{itemize}

\subsubsection{XML-example}
\begin{verbatim}
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <!-- De invoer m.b.t. de observaties -->
  <CTA_SOBS id="obs_assim" database="obs_advec1d.sql"/>

  <!-- Used model class -->
  <xi:include href="models/advec1d_cls.xml"/>

  <!-- De invoer m.b.t. het filter.  -->
  <!-- Dit filter maakt zelf de model-instantiaties aan -->
  <method name="ensemble"> 
  <ensemble_filter>
    <CTA_TIME id="simulation_span" start="0" stop="5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>
    <modes max_id="50">
       <model>
       </model>
    </modes>
  </ensemble_filter>
  </method> 
</costa>


\end{verbatim}


\subsection{Data assimilation  with the Ensemble Square-Root method (ENSRF)}

%   \subsubsection{Kalman filtering with Ensemble Square-Root filter}
       There are two fundamental problems associated with the use of EnKF.
       First is that the ensemble size is limited by the computational cost of
       applying the forecast model to each ensemble member. The second one is
       that small ensembles have few degrees of freedom available to represent
       errors and suffer from sampling errors that will further degrade
       the forecast error covariance representation. Sampling errors lead to loss
       of accuracy and underestimation of error covariances. This problem can
       progressively worsen, resulting in filter divergence.

       In ensemble square-root filters (ENSRF), the analysis step is done
       deterministically without generating any observation noise realization
       (see Tippett,et.al,2003; Evensen, 2004). Since no random sample is
       generated, this extra source of sampling error is eliminated. Therefore,
       these methods are expected to perform better than the ones with
       perturbed observations for a certain type of applications.

       Dropping the time index $k$ for simplicity, the covariance analysis update of the ENSRF is obtained by rewriting the covariance as 
       \begin{eqnarray}
          P^a = E^a E^{a'}=[I-P^f H'(H P^f H' + R)^{-1} H]P^f \\
          = E^f[I-E^{f'}H'(HE^f E^{f'}H'+R)^{-1}HE^f] E^{f'} \\
          = E^f [I- V D^{-1} V'] E^{f'}
       \end{eqnarray}
       where $V=(HL^f)'$ and $D=V' V + R$. Then from this equation it is clear that the analysis ensemble can be calculated from
       \begin{equation}
       \label{eq.ENSRF_EX}
           E^a = E^f X
       \end{equation}
       where $X X'=(I-V D^{-1} V')$. Therefore one can say that the updated ensemble $E^a$ is a linear combination of the columns of $E^f$ and is obtained by inverting the matrix $D$ and computing a matrix square root $X$ of the matrix $[I- V D^{-1} V']$. Note that $X$ can also be replaced by $XU$, where $U$ is an arbitrary orthogonal matrix, so that $(XU)(XU)'=XX'$.

       As we have seen that in ENSRF the analysis step consists of determining the \emph{transformation matrix}, $X$. A number of methods are available to compute $X$. The method implemented in this study is derived as follows:
       \begin{equation}
       \label{eq.xx1}
          X X'=(I-V D^{-1} V') = (I-V (V' V + R)^{-1} V')
       \end{equation}
       If we write $R=SS'$ then we can rewrite Equation \ref{eq.xx1} as
       \begin{equation}
       \label{eq.xx2}
          X X'= \Psi' (\Psi \Psi' + I)^{-1} \Psi)
       \end{equation}
       where $\Psi=S^{-1}V=S^{-1}HE^f$. When computing the singular value decomposition of $\Psi$, i.e. $\Psi=\Gamma \Sigma \Lambda'$, equation (\ref{eq.xx2}) can be written as
       \begin{eqnarray}
          X X'= (\Gamma \Sigma \Lambda')' ((\Gamma \Sigma \Lambda') (\Gamma \Sigma \Lambda')' + I)^{-1} (\Gamma \Sigma \Lambda')) \\
          = \Lambda (I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma) \Lambda' \\
          = (\Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)})(\Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)})'
       \end{eqnarray}
       Thus, a solution of the transformation matrix $X$ is given by
       \begin{equation}
           X = \Lambda \sqrt{(I-\Sigma' (\Sigma \Sigma'+I)^{-1} \Sigma)}
       \end{equation}
%       Moreover, the transformation matrix, $T$, for updating the ensembles can be shown to be
%       \begin{equation}
%           
%       \end{equation}



\subsubsection{Using the ENSRF method}

         The user specifies the number of ensembles in the input file. Like in
         COFFEE filter, the ensembles in ENSRF filter are also stored as
         instances of COSTA models. An array of COSTA state vectors is also
         used for representing matrix $E$ in equation (\ref{eq.ENSRF_EX}).
         There is always a decision to make whether to store $E$ as a COSTA
         matrix or as an array of COSTA state vectors. We inclined sometimes to
         store it as a COSTA matrix, since the ENSRF algorithm consists of
         mostly matrix operations. However, since the model state is stored as
         COSTA state vector it is easier to store $E$ as an array of COSTA
         state vector. Moreover, there is also a function available in COSTA
         for performing linear algebra operation between COSTA matrix and COSTA
         state, i.e. cta\_state\_gemm. This function is used for example in
         computing the correction for all ensemble members, which is done by
         multiplying the transformation matrix and the forecast ensembles.

         The implementation of ENSRF filter requires a singular value
         decomposition (SVD). At the moment there is no COSTA function, which
         can perform SVD directly to a COSTA matrix nor to an array of COSTA
         state vectors. Therefore in this study this is done by assigning the
         values of the COSTA variables to a Fortran array and use the LAPACK
         function DGESVD to perform the SVD.

         This method also requires the square root of observation error
         covariance $R^{1/2}$. However, as also mentioned in the previous
         subsection, the function cta\_sobs\_getcovmat gives the covariance
         matrix $R$. Moreover, there is no COSTA function yet available for
         computing the square root of a square matrix like for example Cholesky
         decomposition. However since in the example models the observation
         noise is always independent, we can easily compute $R^{1/2}$ by taking
         square root of its diagonal elements.

\subsubsection{The configuration of the ENSRF method}
\begin{itemize}
\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of ensembles
\item {\tt model}: Model configuration
\end{itemize}

\subsubsection{XML-example}
 \begin{verbatim}
<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <CTA_SOBS id="obs_assim" database="obs_oscill.sql"/>
  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />
  <method name="ensrf"> 
<ensrf_filter>
  <!-- Simulatie timespan en stapgrootte via CTA_Time -->
  <CTA_TIME id="simulation_span" start="0" stop="5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>

     <modes max_id="50">
        <model>
          <xi:include href="models/oscill_sp.xml"/>
        </model>        
     </modes>
</ensrf_filter>
  
  </method> 
</costa>
 \end{verbatim}




\subsection{Data assimilation  with the COFFEE method}

%      \subsubsection{Kalman filtering with COFFEE filter}
       COFFEE (Complementary Orthogonal subspace Filter For Efficient
       Ensembles) is a hybrid filter, which combines the RRSQRT filter and the
       EnKF (Heemink, et.al., 2001). One problem of the RRSQRT algorithm is
       that repeated projection on the leading eigenvalues leads to a
       systematic bias in forecast errors. Because of the truncation, the
       covariance matrix is always underestimated, which may result in a filter
       divergence problem. The truncated part of the covariance matrix does not
       contribute to the improvements of the state estimate. The COFFEE filter
       attempts to solve this problem by representing the truncated part of the
       covariance matrix as random ensembles and to add them to the EnKF part.
       The RRSQRT part acts as a variance reductor for the ensemble filter,
       thus reducing the statistical errors of the Monte Carlo approach.
       Moreover, by embedding the reduced-rank filter in an EnKF the covariance
       is not underestimated, eliminating the filter divergence problems of the
       reduced-rank approach (also for very small numbers of $q$).

       The COFFEE algorithm can be summarized as follows:
       \begin{enumerate}
        \item Initialization
              \begin{equation}
              \label{eq.COFFEE_initialmatrix}
                 [L^a(0) E^a(0)]=[l_1^a(0),\cdots,l_q^a(0),\xi_1^a(0),\cdots,\xi_N^a(0)]
              \end{equation}
              where $l_i^a$ denotes the $q$ leading eigenvectors of the initial covariance matrix $P_o$. The random ensemble members are generated only to represent the truncated part of the covariance matrix $P_o-L^a(0)L^a(0)'$.
        \item Forecast step

              The $L^a$ is updated using the RRSQRT update as follows
              \begin{eqnarray}
               x^f(k)=M[x^a(k-1)] \\
               l^f_i(k)=\frac{1}{\epsilon} \{ M[x^a(k-1)+ \epsilon  l^a_i(k-1)]-M[x^a(k-1)] \} \label{eq.COFFEE_lf} \\
               \tilde{L}^f(k)=[l_1^f(0),\cdots,l_q^f(0),Q(k-1)^{1/2}] \label{eq.COFFEE_Lf}\\
               L^f(k)=\Pi^f(k)\tilde{L}^f(k)
              \end{eqnarray}
              where $\epsilon$ represents a perturbation, often chosen close to 1, $\Pi^f(k)$ is a projection onto the $q$ leading eigenvectors of the matrix $\tilde{L}^f(k)\tilde{L}^f(k)'$.

              The ensemble $\xi_i^a$ is updated using similar equations with the EnKF update
              \begin{eqnarray}
               \xi_i^f(k)=M[\xi_i^a(k-1)] \\
               x^f(k)=\frac{1}{N} \sum_{i=1}^N \xi_i^f(k) \\
               E^f(k)=[\xi_1^f(k)-x^f(k)+\eta_1(k),\cdots,\xi_N^f(k)-x^f(k)+\eta_N(k)]
              \end{eqnarray}
              where $\eta_i$ are the ensembles representing the truncated part of the covariance matrix with $E[\eta_i(k)\eta_i(k)']=[I-\Pi^f(k)]\tilde{L}^f(k) \tilde{L}^f(k)'[I-\Pi^f(k)]'$.

        \item Analysis step

              In the analysis step the gain matrix $K$ is computed using
              \begin{eqnarray}
                P^f(k)=L^f(k)L^f(k)' + \frac{1}{N-1}E(k)E(k)' \\
                K(k) = P^f(k)H(k)'[H(k)P^f(k)H(k)'+R(k)]^{-1} \\
                x^a(k) = x^f(k) + K(k) [y^o(k) - H(k) x^f(k)] \\
                \tilde{L}^a(k)=\{[I - K(k) H(k)] L_f(k), K(k) R(k)^{1/2}\} \label{eq.COFFEE_La} \\
                L^a(k)=\Pi^a(k)\tilde{L}^a(k) \\
                \xi^a_i(k) = \xi^f_i(k) + K(k) [y^o(k) - H(k) \xi^f_i(k) + \nu_i(k)]
              \end{eqnarray}
       \end{enumerate}



\subsubsection{Using the COFFEE method}
         The COFFEE filter requires the user to specify in the input file the
         number of modes for the RRSQRT part and the number of ensembles for
         the EnKF part. Moreover, the perturbation $\delta$ (i.e. $\epsilon$ in
         equation (\ref{eq.COFFEE_lf})) needs also to be specified.

         In this implementation, the modes and ensembles are represented as COSTA model instances. Here each model instance represent one mode or ensemble. This is chosen for supporting parallel computation implementation in the future.

         In the implementation of COFFEE filter, the initial matrix $[L^a(0)
           E^a(0)]$ in equation \ref{eq.COFFEE_initialmatrix} is always assumed
         to be zero. This however should be extended to accommodate the option
         where the user can specify this from the input file. Moreover, the size of
         $L$ is limited to a maximum of 200. This is due to the fact that its
         dimension varies according to the number of noise parameters in the
         model as well as to the number of observations. Since we know the
         number of observations only at the analysis-step and that the number
         of observations may vary, it is not possible to determine the
         dimension of $L$ a priori. Furthermore, in the present implementation
         the computation cost increases with number of observations. More study
         is required to develop the implementation for solving this problem. It
         should be noted also that the present implementation does not support
         covariance localization using Schur product either.

         The COFFEE filter also requires the square root of the model error
         covariance matrix $Q^{1/2}$ in equation (\ref{eq.COFFEE_Lf}) as well
         as the observation error covariance matrix $R^{1/2}$ in equation
         (\ref{eq.COFFEE_La}). The covariance matrices $Q$ and $R$ can be
         obtained by using the functions cta\_model\_getnoisecovar and
         cta\_sobs\_getcovmat respectively. However since we mostly work with
         the square root of the covariance matrices, it is better if these
         functions are modified to give the square root matrices. Note that
         cta\_model\_getnoisecovar gives the noise covariance matrix in terms
         of an array of COSTA state vectors, while cta\_sobs\_getcovmat gives
         the observation noise covariance matrix in term of COSTA matrix. This
         may cause confusion for new programmers working with COSTA. Moreover,
         it may be better to use a similar name for the two functions. This
         makes the functions easier to remember.

         Note that for general nonlinear models, the additional columns of
         matrix $L^f$ with matrix $Q^{1/2}$ in equation (\ref{eq.COFFEE_Lf})
         can be computed using finite difference as follows
         \begin{eqnarray}
             \hat{x}(k) = M[x(k-1),u(k-1),0] \\
             x^n_i(k) = M[x(k-1),u(k-1),\epsilon Q^{1/2}(:,i)] \\
             l^n_i(k) =  \frac{1}{\epsilon} (x_i(k) - \hat{x}(k))\\
             \tilde{L}^f(k)=[l_1^f(k),\cdots,l_q^f(k),l^n_1(k),...,l^n_w(k)]
         \end{eqnarray}
         where $M$ is the general nonlinear model operator, which is a function
         of the previous time-step state vector $x$, the input forcing $u$ and
         the noise $\eta$. The state vector $\hat{x}$ is obtained from the
         deterministic model, i.e. when no noise is present, while $x^n_i$ is
         from the stochastic model, with the realization $\eta=\epsilon
         Q^{1/2}(:,i)$, where $Q^{1/2}(:,i)$ refers to the $i^{th}$ column
         vector of matrix $Q^{1/2}$. The implementation requires a function
         which can set the model noise realization. However, since there is no
         such function available yet, it is not possible to
         implement this in COSTA at the moment.



\subsubsection{The configuration of the COFFEE method}
\begin{itemize}

\item {\tt simulation\_span}: the overall timespan to run the model
\item {\tt output}: output specification
\item {\tt modes}: The number of modes
\item {\tt delta}: perturbation, see equation \ref{eq.COFFEE_lf}
\item {\tt model}: Model configuration
\item {\tt ensembles}: the number of ensembles
\end{itemize}

\subsubsection{XML-example}

\begin{verbatim}

<costa xmlns:xi="http://www.w3.org/2001/XInclude">
  <CTA_SOBS id="obs_assim" database="obs_oscill.sql"/>

  <!-- Used model class -->
  <CTA_MODELCLASS id="modelclass"  name="CTA_MODBUILD_SP" />


  <method name="coffee"> 
<coffee_filter>
  <CTA_TIME id="simulation_span" start="0" stop=".5" step="0.0005"/>
    <output>
       <filename>results.m</filename>
       <CTA_TIME id="times_mean" start="0" stop="50" step="0.05"/>
    </output>

     <modes max_id="10">
        <delta id="1E-4"> </delta>
        <model>
          <xi:include href="models/oscill_sp.xml"/>
        </model>        
     </modes>

     <ensembles max_id="5"> </ensembles>
</coffee_filter>
  
  </method> 
</costa>
\end{verbatim}


% \bibliographystyle{plain}
% \bibliography{\styles/costa}

\end{document}

