Cycling
=======
DUBFI does not provide a cycling loop.
But if provides the option to initialize from previous output.
This shall make it easy for external tools to implement a cycling loop.

When cycling the inversion, we distinguish the initial prior from which the cycling starts, and the prior of one cycling step.
The latter is constructed from the combination of the initial prior and the posterior of the previous cycling step.

Prior scaling factors
---------------------
DUBFI assumes that the initial prior scaling factors are 0.
For the prior in a cycling step, we assume that :math:`s^\text{prior} = \alpha s^\text{init} + (1 - \alpha) s^\text{prev. posterior}`
where :math:`\alpha` denotes the weight of the initial prior (configuration entry :code:`cycle.initial_prior_weight`) and :math:`s^\text{init}=0`.

Prior B matrix
--------------
We need to find the error covariance matrix :math:`B` of :math:`s^\text{prior}`.
Assume that :math:`s^\text{prev. posterior}` (covariance matrix :math:`B^p`) and :math:`s^\text{init}` (:math:`B^i`) have correlated uncertainties.
We can write the joint error covariance matrix as

.. math::
   :label: B-extended

   \begin{pmatrix}
     B^i & B^{ip} \\
     B^{pi} & B^p
   \end{pmatrix}

where :math:`B^{pi}={B^{ip}}^\top` describes the correlations.
We obtain :math:`B = \alpha^2 B^i + (1-\alpha)^2 B^p + \alpha (1 - \alpha) (B^{ip} + B^{pi})`.

As an error covariance matrix, :eq:`B-extended` must be real-symmetric and positive semi-definite.
Since each cycling step adds information to the initial state, the information in :math:`s^\text{init}` and :math:`s^\text{prev. posterior}` is maximally correlated.
This maximal correlation is obtained by choosing :math:`B^{ip}=\sqrt{B^i}\sqrt{B^p}` such that

.. math::

   \begin{pmatrix}
     B^i & B^{ip} \\
     B^{pi} & B^p
   \end{pmatrix}
   =
   \begin{pmatrix}
     \sqrt{B^i} & 0 \\
     0 & \sqrt{B^p}
   \end{pmatrix}
   \begin{pmatrix}
     I & I \\
     I & I
   \end{pmatrix}
   \begin{pmatrix}
     \sqrt{B^i} & 0 \\
     0 & \sqrt{B^p}
   \end{pmatrix}

where :math:`I` denotes the identity matrix.
This leads to the relation

.. math::
   B = \alpha^2 B^i + (1-\alpha)^2 B^p + \alpha (1 - \alpha) (\sqrt{B^i}\sqrt{B^p} + \sqrt{B^p}\sqrt{B^i})

that is implemented in :py:func:`dubfi.fluxes.util.mix_b` and used in :py:func:`dubfi.fluxes.util.get_prior_cycling`.

.. note::

   Choosing :math:`B=\alpha B^i + (1-\alpha) B^p` is in general not consistent with positivity of :eq:`B-extended`.
   Consider :math:`\alpha=\tfrac{1}{2}`, :math:`B^i=3`, and :math:`B^p=1`.
   Choosing :math:`B=\alpha B^i + (1-\alpha) B^p = 2` implies that :math:`B^{ip}=2`, but the matrix :math:`\begin{pmatrix} 3 & 2 \\ 2 & 1 \end{pmatrix}` is not positive semidefinite.
   Thus, :math:`B=2` is inconsistent with the uncertainties assumed for :math:`s^\text{init}` and :math:`s^\text{prev. posterior}`.