# SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
# SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
# SPDX-FileContributor: 2003-24 Bradley M. Bell
# ----------------------------------------------------------------------------
{xrst_begin ipopt_nlp_ode_problem dev}
{xrst_spell
  ny
  nz
}

An ODE Inverse Problem Example
##############################

Notation
********
The table below contains
the name of a variable,
the meaning of the variable value,
and the value for this particular example.
If the value is not specified in the table below,
the corresponding value in ``ipopt_nlp_ode_problem.hpp``
can be changed and the example should still run
(with no other changes).

.. csv-table::
   :widths: auto

   **Name**,**Meaning**,**Value**
   :math:`Na`,number of parameters to fit,3
   :math:`Ny`,number components in ODE,2
   :math:`Nz`,number of measurements,4
   :math:`N(i)`,# of grid points between *i-1*-th and *i*-th measurement,
   :math:`S(i)`,# of grid points up to an including *i*-th measurement,

Forward Problem
***************
We consider the following ordinary differential equation:

.. math::
   :nowrap:

   \begin{eqnarray}
      \partial_t y_0 ( t , a ) & = & - a_1 * y_0 (t, a )
      \\
      \partial_t y_1 (t , a ) & = & + a_1 * y_0 (t, a ) - a_2 * y_1 (t, a )
   \end{eqnarray}

with the initial conditions

.. math::

   y_0 (0 , a) = F(a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)

where :math:`Na` is the number of parameters,
:math:`a \in \B{R}^{Na}` is an unknown parameter vector.
The function and :math:`F : \B{R}^{Na} \rightarrow \B{R}^{Ny}`
is defined by the equation above
where :math:`Ny` is the number of components in :math:`y(t, a)`.
Our forward problem is stated as follows: Given :math:`a \in \B{R}^{Na}`
determine the value of :math:`y ( t , a )`, for :math:`t \in R`, that
solves the initial value problem above.

Measurements
************
We use :math:`Nz` to denote the number of measurements.
Suppose we are also given a measurement vector :math:`z \in \B{R}^{Nz}`
and for :math:`i = 1, \ldots, Nz`,
we model :math:`z_i` by

.. math::

   z_i = y_1 ( s_i , a) + e_i

where :math:`s_i \in \B{R}` is the time for the *i*-th measurement,
:math:`e_i \sim {\bf N} (0 , \sigma^2 )` is the corresponding noise,
and :math:`\sigma \in \B{R}_+` is the corresponding standard deviation.

Simulation Analytic Solution
============================
The following analytic solution to the forward problem is used
to simulate a data set:

.. math::
   :nowrap:

   \begin{eqnarray}
      y_0 (t , a) & = & a_0 * \exp( - a_1 * t )
      \\
      y_1 (t , a) & = &
      a_0 * a_1 * \frac{\exp( - a_2 * t ) - \exp( -a_1 * t )}{ a_1 - a_2 }
   \end{eqnarray}

Simulation Parameter Values
===========================

.. list-table::
   :widths: auto

   * - :math:`\bar{a}_0 = 1`
     - initial value of :math:`y_0 (t, a)`
   * - :math:`\bar{a}_1 = 2`
     - transfer rate from compartment zero to compartment one
   * - :math:`\bar{a}_2 = 1`
     - transfer rate from compartment one to outside world
   * - :math:`\sigma = 0`
     - standard deviation of measurement noise
   * - :math:`e_i = 0`
     - simulated measurement noise, :math:`i = 1 , \ldots , Nz`
   * - :math:`s_i = i * .5`
     - time corresponding to the *i*-th measurement,
       :math:`i = 1 , \ldots , Nz`

Simulated Measurement Values
============================
The simulated measurement values are given by the equation

.. math::
   :nowrap:

   \begin{eqnarray}
   z_i
   & = &  e_i + y_1 ( s_i , \bar{a} )
   \\
   & = &
   e_i + \bar{a}_0 * \bar{a}_1 *
      \frac{\exp( - \bar{a}_2 * s_i ) - \exp( -\bar{a}_1 * s_i )}
         { \bar{a}_1 - \bar{a}_2 }
   \end{eqnarray}

for :math:`k = 1, \ldots , Nz`.

Inverse Problem
***************
The maximum likelihood estimate for :math:`a` given :math:`z`
solves the following inverse problem

.. math::
   :nowrap:

   \begin{eqnarray}
   {\rm minimize} \;
      & \sum_{i=1}^{Nz} H_i [ y( s_i , a ) , a ]
      & \;{\rm w.r.t} \; a \in \B{R}^{Na}
   \end{eqnarray}

where the functions
:math:`H_i : \B{R}^{Ny} \times \B{R}^{Na} \rightarrow \B{R}` is
defined by

.. math::

   H_i (y, a) = ( z_i - y_1 )^2

Trapezoidal Approximation
*************************
This example uses a trapezoidal approximation to solve the ODE.
This approximation procedures starts with

.. math::

   y^0 = y(0, a) = \left( \begin{array}{c} a_0 \\ 0 \end{array} \right)

Given a time grid :math:`\{ t_i \}` and
an approximate value :math:`y^{i-1}` for :math:`y ( t_{i-1} , a )`,
the a trapezoidal method approximates :math:`y ( t_i , a )`
(denoted by :math:`y^i` ) by solving the equation

.. math::

   y^i  =  y^{i-1} +
   \left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{t_i - t_{i-1} }{ 2 }

where :math:`G : \B{R}^{Ny} \times \B{R}^{Na} \rightarrow \B{R}^{Ny}` is the
function representing this ODE; i.e.

.. math::

   G(y, a) = \left(  \begin{array}{c}
      - a_1 * y_0
      \\
      + a_1 * y_0 - a_2 * y_1
   \end{array} \right)

This :math:`G(y, a)` is linear with respect to :math:`y`, hence
the implicit equation defining :math:`y^i` can be solved
inverting the a set of linear equations.
In the general case,
where :math:`G(y, a)` is non-linear with respect to :math:`y`,
an iterative procedure is used to calculate :math:`y^i`
from :math:`y^{i-1}`.

Trapezoidal Time Grid
=====================
The discrete time grid, used for the trapezoidal approximation, is
denoted by :math:`\{ t_i \}` which is defined by:
:math:`t_0 = 0` and
for :math:`i = 1 , \ldots , Nz` and for :math:`j = 1 , \ldots , N(i)`,

.. math::
   :nowrap:

   \begin{eqnarray}
      \Delta t_i & = & ( s_i - s_{i-1} ) / N(i)
      \\
      t_{S(i-1)+j} & = & s_{i-1} + \Delta t_i * j
   \end{eqnarray}

where :math:`s_0 = 0`,
:math:`N(i)` is the number of time grid points between
:math:`s_{i-1}` and :math:`s_i`,
:math:`S(0) = 0`, and
:math:`S(i) = N(1) + \ldots + N(i)`.
Note that for :math:`i = 0 , \ldots , S(Nz)`,
:math:`y^i` denotes our approximation for :math:`y( t_i , a )`
and :math:`t_{S(i)}` is equal to :math:`s_i`.

Black Box Method
****************
A common approach to an inverse problem is to treat the forward problem
as a black box (that we do not look inside of or try to understand).
In this approach, for each value of the parameter vector :math:`a`
one uses the
:ref:`ipopt_nlp_ode_problem@Trapezoidal Approximation`
(on a finer grid that :math:`\{ s_i \}`)
to solve for :math:`y_1 ( s_i , a )` for :math:`i = 1 , \ldots , Nz`.

Two levels of Iteration
=======================
As noted above, the trapezoidal approximation
often requires an iterative procedure.
Thus, in this approach, there are two levels of iterations,
one with respect to the parameter values during the minimization
and the other for solving the trapezoidal approximation equation.

Derivatives
===========
In addition, in the black box approach, differentiating the ODE solution
often involves differentiating an iterative procedure.
Direct application of AD to compute these derivatives
requires a huge amount of memory and calculations to differentiate the
iterative procedure.
(There are special techniques for applying AD to the solutions of iterative
procedures, but that is outside the scope of this presentation).

Simultaneous Method
*******************
The simultaneous forward and inverse method
uses constraints to include the solution of
the forward problem in the inverse problem.
To be specific for our example,

.. math::
   :nowrap:

   \begin{eqnarray}
   {\rm minimize}
   & \sum_{i=1}^{Nz} H_i ( y^{N(i)} , a )
   & \; {\rm w.r.t} \; y^1 \in \B{R}^{Ny} , \ldots , y^{S(Nz)} \in \B{R}^{Ny} ,
     \; a \in \B{R}^{Na}
   \\
   {\rm subject \; to}
      & y^j  =  y^{j-1} +
   \left[ G( y^{j-1} , a ) + G( y^j , a ) \right] * \frac{ t_j - t_{j-1} }{ 2 }
      & \; {\rm for} \; j = 1 , \ldots , S(Nz)
   \\
      & y^0 = F(a)
   \end{eqnarray}

where for :math:`i = 1, \ldots , Nz`,
:math:`N(i)` is the number of time intervals between
:math:`s_{i-1}` and :math:`s_i` (with :math:`s_0 = 0`)
and :math:`S(i) = N(1) + \ldots + N(i)`.
Note that, in this form, the iterations of the optimization procedure
also solve the forward problem equations.
In addition, the functions that need to be differentiated
do not involve an iterative procedure.
{xrst_toc_hidden
   cppad_ipopt/example/ode_problem.hpp
}
Source
******
The file ``ipopt_nlp_ode_problem.hpp`` contains
source code that defines the example values and functions defined above.

{xrst_end ipopt_nlp_ode_problem}
-----------------------------------------------------------------------------
{xrst_begin ipopt_nlp_ode_simple dev}
{xrst_spell
  fg
  ny
  nz
}

ODE Fitting Using Simple Representation
#######################################

Purpose
*******
In this section we represent the objective and constraint functions,
(in the simultaneous forward and reverse optimization problem)
using the :ref:`cppad_ipopt_nlp@Simple Representation`
in the sense of ``cppad_ipopt_nlp`` .

Argument Vector
***************
The argument vector that we are optimizing with respect to
( :math:`x` in :ref:`cppad_ipopt_nlp-name` )
has the following structure

.. math::

   x = ( y^0 , \cdots , y^{S(Nz)} , a )

Note that :math:`x \in \B{R}^{S(Nz) + Na}` and

.. math::
   :nowrap:

   \begin{eqnarray}
      y^i & = & ( x_{Ny * i} , \ldots ,  x_{Ny * i + Ny - 1} )
      \\
      a   & = & ( x_{Ny *S(Nz) + Ny} , \ldots , x_{Ny * S(Nz) + Na - 1} )
   \end{eqnarray}

Objective Function
******************
The objective function
( :math:`fg_0 (x)` in :ref:`cppad_ipopt_nlp-name` )
has the following representation,

.. math::

   fg_0 (x) = \sum_{i=1}^{Nz} H_i ( y^{S(i)} , a )

Initial Condition Constraint
****************************
For :math:`i = 1 , \ldots , Ny`,
we define the component functions :math:`fg_i (x)`,
and corresponding constraint equations, by

.. math::

   0 = fg_i ( x ) = y_i^0 - F_i (a)

Trapezoidal Approximation Constraint
************************************
For :math:`i = 1, \ldots , S(Nz)`,
and for :math:`j = 1 , \ldots , Ny`,
we define the component functions :math:`fg_{Ny*i + j} (x)`,
and corresponding constraint equations, by

.. math::

   0 = fg_{Ny*i + j } = y_j^{i}  -  y_j^{i-1} -
      \left[ G_j ( y^i , a ) + G_j ( y^{i-1} , a ) \right] *
         \frac{t_i - t_{i-1} }{ 2 }

{xrst_toc_hidden
   cppad_ipopt/example/ode_simple.hpp
}
Source
******
The file ``ipopt_nlp_ode_simple.hpp``
contains source code for this representation of the
objective and constraints.

{xrst_end ipopt_nlp_ode_simple}
-----------------------------------------------------------------------------
{xrst_begin ipopt_nlp_ode_fast dev}
{xrst_spell
  fg
  ny
  nz
}

ODE Fitting Using Fast Representation
#####################################

Purpose
*******
In this section we represent a more complex representation of the
simultaneous forward and reverse ODE fitting problem (described above).
The representation defines the problem using
simpler functions that are faster to differentiate
(either by hand coding or by using AD).

Objective Function
******************
We use the following representation for the
:ref:`ipopt_nlp_ode_simple@Objective Function` :
For :math:`k = 0 , \ldots , Nz - 1`,
we define the function :math:`r^k : \B{R}^{Ny+Na} \rightarrow \B{R}`
by

.. math::
   :nowrap:

   \begin{eqnarray}
   fg_0 (x) & = & \sum_{i=1}^{Nz} H_i ( y^{S(i)} , a )
   \\
   fg_0 (x) & = & \sum_{k=0}^{Nz-1} r^k ( u^{k,0} )
   \end{eqnarray}

where for :math:`k = 0 , \ldots , Nz-1`,
:math:`u^{k,0} \in \B{R}^{Ny + Na}` is defined by
:math:`u^{k,0} =   ( y^{S(k+1)} , a )`

Range Indices I(k,0)
====================
For :math:`k = 0 , \ldots , Nz - 1`,
the range index in the vector :math:`fg (x)`
corresponding to :math:`r^k ( u^{k,0} )` is 0.
Thus, the range indices are given by
:math:`I(k,0) = \{ 0 \}` for :math:`k = 0 , \ldots , Nz-1`.

Domain Indices J(k,0)
=====================
For :math:`k = 0 , \ldots , Nz - 1`,
the components of the vector :math:`x`
corresponding to the vector :math:`u^{k,0}` are

.. math::
   :nowrap:

   \begin{eqnarray}
   u^{k,0} & = & ( y^{S(k+1} , a )
   \\
   & = &
   (   x_{Ny * S(k+1)} \; , \;
      \ldots \; , \;
      x_{Ny * S(k+1) + Ny - 1} \; , \;
      x_{Ny * S(Nz) + Ny } \; , \;
      \ldots \; , \;
      x_{Ny * S(Nz) + Ny + Na - 1}
   )
   \end{eqnarray}

Thus, the domain indices are given by

.. math::

   J(k,0) = \{
      Ny * S(k+1) \; , \;
      \ldots  \; , \;
      Ny * S(k+1) + Ny - 1 \; , \;
      Ny * S(Nz) + Ny \; , \;
      \ldots  \; , \;
      Ny * S(Nz) + Ny + Na - 1
   \}

Initial Condition
*****************
We use the following representation for the
:ref:`ipopt_nlp_ode_simple@Initial Condition Constraint` :
For :math:`k = Nz` we define the function
:math:`r^k : \B{R}^{Ny} \times \B{R}^{Na + Ny}` by

.. math::
   :nowrap:

   \begin{eqnarray}
      0 & = & fg_i ( x ) = y_i^0 - F_i (a)
      \\
      0 & = & r_{i-1}^k ( u^{k,0} ) = y_i^0 - F_i(a)
   \end{eqnarray}

where :math:`i = 1 , \ldots , Ny` and
where :math:`u^{k,0} \in \B{R}^{Ny + Na}` is defined by
:math:`u^{k,0}  = ( y^0 , a)`.

Range Indices I(k,0)
====================
For :math:`k = Nz`,
the range index in the vector :math:`fg (x)`
corresponding to :math:`r^k ( u^{k,0} )` are
:math:`I(k,0) = \{ 1 , \ldots , Ny \}`.

Domain Indices J(k,0)
=====================
For :math:`k = Nz`,
the components of the vector :math:`x`
corresponding to the vector :math:`u^{k,0}` are

.. math::
   :nowrap:

   \begin{eqnarray}
   u^{k,0} & = & ( y^0 , a)
   \\
   & = &
   (   x_0 \; , \;
      \ldots \; , \;
      x_{Ny-1} \; , \;
      x_{Ny * S(Nz) + Ny } \; , \;
      \ldots \; , \;
      x_{Ny * S(Nz) + Ny + Na - 1}
   )
   \end{eqnarray}

Thus, the domain indices are given by

.. math::

   J(k,0) = \{
      0 \; , \;
      \ldots  \; , \;
      Ny - 1 \; , \;
      Ny * S(Nz) + Ny \; , \;
      \ldots  \; , \;
      Ny * S(Nz) + Ny + Na - 1
   \}

Trapezoidal Approximation
*************************
We use the following representation for the
:ref:`ipopt_nlp_ode_simple@Trapezoidal Approximation Constraint` :
For :math:`k = 1 , \ldots , Nz`,
we define the function :math:`r^{Nz+k} : \B{R}^{2*Ny+Na} \rightarrow \B{R}^{Ny}` by

.. math::

   r^{Nz+k} ( y , w , a )
   =
   y - w -  [ G( y , a ) + G( w , a ) ] * \frac{ \Delta t_k }{ 2 }

For :math:`\ell = 0 , \ldots , N(k)-1`,
using the notation :math:`i = Ny * S(k-1) + \ell + 1`,
the corresponding trapezoidal approximation is represented by

.. math::
   :nowrap:

   \begin{eqnarray}
   0 & = & fg_{Ny+i} (x)
   \\
   & = &
   y^i  -  y^{i-1} -
   \left[ G( y^i , a ) + G( y^{i-1} , a ) \right] * \frac{\Delta t_k }{ 2 }
   \\
   & = &
   r^{Nz+k} ( u^{Nz+k , \ell} )
   \end{eqnarray}

where :math:`u^{Nz+k,\ell} \in \B{R}^{2*Ny + Na}` is defined by
:math:`u^{Nz+k,\ell}  = ( y^{i-1} , y^i , a)`.

Range Indices I(k,0)
====================
For :math:`k = Nz + 1, \ldots , 2*Nz`,
and :math:`\ell = 0 , \ldots , N(k)-1`,
the range index in the vector :math:`fg (x)`
corresponding to :math:`r^k ( u^{k,\ell} )` are
:math:`I(k,\ell) = \{ Ny + i , \ldots , 2*Ny + i - 1  \}`
where :math:`i = Ny * S(k-1) + \ell + 1`.

Domain Indices J(k,0)
=====================
For :math:`k = Nz + 1, \ldots , 2*Nz`,
and :math:`\ell = 0 , \ldots , N(k)-1`,
define :math:`i = Ny * S(k-1) + \ell + 1`.
The components of the vector :math:`x`
corresponding to the vector :math:`u^{k,\ell}` are
(and the function :math:`fg (x)` in :ref:`cppad_ipopt_nlp-name` )

.. math::
   :nowrap:

   \begin{eqnarray}
   u^{k, \ell} & = & ( y^{i-1} , y^i , a )
   \\
   & = &
   (   x_{Ny * (i-1)} \; , \;
      \ldots \; , \;
      x_{Ny * (i+1) - 1} \; , \;
      x_{Ny * S(Nz) + Ny } \; , \;
      \ldots \; , \;
      x_{Ny * S(Nz) + Ny + Na - 1}
   )
   \end{eqnarray}

Thus, the domain indices are given by

.. math::

   J(k,\ell) = \{
      Ny * (i-1) \; , \;
      \ldots  \; , \;
      Ny * (i+1) - 1 \; , \;
      Ny * S(Nz) + Ny \; , \;
      \ldots  \; , \;
      Ny * S(Nz) + Ny + Na - 1
   \}

{xrst_toc_hidden
   cppad_ipopt/example/ode_fast.hpp
}
Source
******
The file ``ipopt_nlp_ode_fast.hpp``
contains source code for this representation of the objective
and constraints.

{xrst_end ipopt_nlp_ode_fast}
------------------------------------------------------------------------------