#import "@preview/cetz:0.3.2"
#import "@preview/cetz-plot:0.1.1"
#import "@preview/finite:0.4.1"

#set page(paper: "a4")
#set par(justify: true)
#let breakl(b) = {
  raw(b.text.replace(regex(" *\\\\\n *"), " "), lang: b.lang, block: b.block)
}

#show raw: set text(font: "CaskaydiaCove NF")
#show link: underline

= Notes on hybrid system simulation

== Questions

=== Internal state abstraction

In `LazySim(S).accelerate`, we need access to the type of the simulation's
internal state in order to use the functions provided by `S`. This requires
another definition of the discrete and continuous node types where the inner
state is visible. Is the existential type really needed in the original
definitions? What guarantees does it provide? Could we do with explicit type
parameters for the state? On the other hand, could we somehow implement
`LazySim(S).accelerate` without neeeding this access?

Should the simulation state abstraction be parameterized by a solver and model
state abstractions? What does it mean to have a functional simulation state when
the solver and/or model states are modified in place?

We provide a state copy function which is assumed to ensure that the closures
returned by the solver remain valid after another call to the solver (i.e., that
no underlying data structure has been modified). In particular, the state copy
function has no obligation to copy the entire state, only what is needed to
ensure the preservation of this property. Its signature `sstate -> sstate` does
not accurately reflect this. Should the state abstraction be decomposed into two
parts, one relevant for this problem, and which should be copied by the `copy`
function, and another, which does not matter?

=== Simulation semantics

What does providing an input whose starting date is not equal to the end date of
the previous input mean (and in particular, in the case where the domains
overlap: $"end"_((n-1)) > "start"_n$)?

We immediately reset the solvers when obtaining a new input. This is (most
likely) very inneficient, as adaptive solvers usually start with a very small
step size, and increase this step size as needed as they gain more knowledge of
the underlying function. Testing this, and thinking about how we might avoid
this automatic reset, would be beneficial.

All tests have, up to now, been made only with "stateless"
#smallcaps([Runge-Kutta]) solvers that do not really suffer loss of precision
from resets. Attempting to run some examples with solvers such as
#smallcaps([Sundials CVODE]) would be beneficial to measure the impact that
resets have on the accuracy of the results.

== The problem

The overarching problem is _transparent assertions_. We want to simulate
assertions with their own solvers, because simulating assertions with the same
solver as their parent system changes the final results, when assertions are
supposed to be transparent.

Assertions are themselves hybrid systems, and can contain their own assertions.
We thus need a semantics that is both *independent* of the solver, which becomes
a parameter, and *recursive*, in that the simulation of an assertion is itself
the simulation of a hybrid system with assertions.

The main idea is as follows: a solver can be seen as a synchronous function (an
inner state, a step function, and a reset function).
```ocaml
type ('p, 'a, 'b) dnode = DNode :
  { state: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } ->
  ('p, 'a, 'b) dnode
```
An ODE solver is, in fact, a special case of a discrete node, which is
initialized with an initial value problem, takes as input a time horizon to
reach, and returns as output an actual time reached and a dense solution.
```ocaml
type ('y, 'yder) ivp = { init : 'y; fder : time -> 'y -> 'yder; stop : time }
type ('y, 'yder) csolver = (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
```
Similarly, a zero-crossing solver is initialized with a zero-crossing problem,
takes as input a horizon and dense function, and returns an actual time reached
and optional zero-crossing.
```ocaml
type ('y, 'zout) zc = { fzer : time -> 'y -> 'zout; yc : 'y; size: int }
type ('y, 'zin, 'zout) zsolver =
  (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
```
A complete solver is then a composition of these two: it is initialized with
both an `ivp` and `zc`, takes in a horizon, and returns an actual time reached,
dense solution and optional zero-crossing.

#breakl(```ocaml
type ('y, 'yder, 'zin, 'zout) solver =
  (('y, 'yder) ivp * ('y, 'zout) zc, time, time * (time -> 'y) * 'zin option) \
  dnode
```)

The construction of a full solver from an ODE and a zero-crossing solver is
available in the file `src/lib/hsim/solver.ml`.

The simulation of a hybrid system with a solver can itself be seen as a discrete
node, operating on streams of functions defined on successive intervals.

```ocaml
type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
```

== The simulation

=== Steps

A simulation step is either a discrete step or a continuous step.

Discrete steps call the `model.step` function, build a constant output function
defined on a singleton interval `[now, now]`, and end in one of four possible
cases:
- We need to cascade (perform another discrete step)
- We reached the horizon (reset the solver? idle mode?)
- We can still integrate, but there has been a state jump (reset solver)
- We can still integrate, and there has been no state jump (perform a continuous
  step)

#columns(4, [#align(center, {
  import cetz: *
  import cetz-plot.plot: *
  let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
  canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((1,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((4, none),), y-grid: "both",
      style: p, {
      add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline")
      add(((1,4),), mark: "o", mark-size: .03)
      add(((1,5),(1.3,4.7),(1.7,3.2),(2,3)), line: "spline")
    })
  })
  colbreak(); canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
      add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
    })
  })
  colbreak(); canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((1, none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((5, none),), y-grid: "both", {
      add(((0,2),(0.5,2.3),(1,3)), line: "spline")
      add(((1,5),(1.5,4.7),(2,4)), line: "spline")
    })
  })
  colbreak(); canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
      add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
    })
  })
})])

Continuous steps call the `solver.step` function, parameterized by `model.fder`
and `model.fzer`, build an output function using `model.fout` defined on the
interval `[now, now + h]`, and end in one of three possible cases:
- We found a zero-crossing (perform a discrete step)
- We found no zero-crossing, but have reached the horizon (perform a discrete
  step)
- We found no zero-crossing and have not reached the horizon (perform another
  continuous step)

#columns(3, [#align(center, {
  import cetz: *
  import cetz-plot.plot: *
  let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
  canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((.7,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: -2, y-ticks: ((0,$0$),), y-grid: "both", {
      add(((0,-1),(.5,-.6),(.7,0),(1,2.5),(2,1),(3,.5)), line: "spline")
    })
  })
  colbreak(); canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((3,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((3, none),), y-grid: "both", {
      add(((0,2),(.5,3.2),(1,5.5),(2,4),(3,3)), line: "spline")
    })
  })
  colbreak(); canvas({
    plot(
      size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
      x-ticks: ((2,none),), x-grid: "both", y-tick-step: none, y-label: none,
      y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
      add(((0,2),(0.5,2.3),(1,3),(1.3,5),(2,4),(3,3)), line: "spline")
    })
  })
})])

=== Modes and output format

There are two simulation modes, depending on whether the solver provides a way
to copy its internal state as part of its interface, or more precisely, whether
the dense solution returned by the solver's step function remains valid after we
modify the solver's state (through another integration step or a full reset).
Both modes output values representing functions defined on real intervals
```ocaml
type continuity = Continuous | Discontinuous
type 'a value   = { start: time; length: time; u: time -> 'a; cont: continuity }
type 'a signal  = 'a value option
```
where `u` is defined on the interval `[start, start + length]` and `cont`
represents the relationship of a value with the next one in the stream: if
`cont = Continuous`, the next value in the stream simply extends the current
one; if `cont = Discontinuous`, the next value does not extend the current one,
either due to a jump, a discrete step or the horizon being reached.

In *_lazy_* mode, the output has format ```ocaml 'b signal = 'b value option```.
The input stream first provides a value for the simulation to run, and must then
provide as many `None` values as needed for the solver to finish integrating as
much as possible (that is, until it reaches a zero-crossing or the horizon
provided by the current input value). This mode is the one used when the solver
modifies its solution in place (for example, #smallcaps([Sundials CVODE])). Each
step of the simulation runs _one_ step of the solver, and then runs all
subsystems as much as they need to in order for them to finish using the output
of the solver. The solver output is then not needed anymore, and we can call the
solver again to progress with the simulation.

We can also optionally "accelerate" the simulation if given a solver supporting
state copies. That is, we can keep calling the solver as long as we keep making
continuous steps. When we reach a zero-crossing or the horizon, we stop calling
the solver and return the concatenation of all the intermediate functions. The
caller can then keep making simulation steps, which will run through all the
discrete steps, until we reach another continuous step sequence, where the next
simulation step will again run through as many continuous steps as it can, and
so on and so forth.

#columns(3, {
  import cetz: *
  import cetz-plot.plot: *
  align(right, canvas({
    plot(
      size: (2, 1), axis-style: "left", x-tick-step: none, x-label: none,
      y-tick-step: none, y-label: none, x-ticks: ((1,none),(2, none)),
      x-grid: "both", {
      add(((0,2),(1,3)))
      add(((1,3),(2,4)))
      add(((2,3.5),), mark: "o", mark-size: .03)
      add(((2,3),(3,3)))
    })
  }))
  colbreak(); linebreak()
  align(center, $ ="accelerate"=> $)
  colbreak(); align(left, canvas({
    plot(
      size: (2, 1), axis-style: "left", x-tick-step: none, x-label: none,
      y-tick-step: none, y-label: none, x-ticks: ((2, none),),
      x-grid: "both", {
      add(((0,2),(2,4))); add(((0,2),))
      add(((2,3.5),), mark: "o", mark-size: .03)
      add(((2,3),(3,3)))
    })
  }))
})

In *_greedy_* mode, the output has format ```ocaml 'b value list```. The input
stream provides a value for the simulation to run, and the simulation then calls
the solver as much as it needs to until it reaches the end of the possible
integration period (zero-crossing or horizon). We concatenate the list of
results (all functions returned from the solver steps are continuous at their
edges and defined on more than a single instant, and we can thus create a large,
piecewise-defined function from all the results up to this point). The solver
then performs as many discrete steps as needed (i.e. we cascade). These are kept
in the output list _as is_. We cannot concatenate a discrete step (i.e. a
function defined on a singleton interval $[t,t]$) with anything else, as it may
be hidden by the (inclusive) border of the other function. For instance,

#align(center, ```ocaml
f = concat { start=0.0; length=1.0; u=fun t -> t   }
           { start=1.0; length=0.0; u=fun t -> 2.0 }
```)

#columns(2, [#{
  align(center, cetz.canvas({
    import cetz-plot.plot: *
    plot(size: (1, 0.5), axis-style: "left", x-tick-step: 1, x-label: "time",
      y-tick-step: 1, y-label: none, {
      add(((0,0),(1,1)))
      add(((1,2),), mark: "o", mark-size: .03)
    })
  }))
  colbreak()
  $ f(t) = cases(
    t   & "if" t in [0, 1],
    2   & "if" t = 1 + partial,
    bot & "otherwise",
  ) $
}])

Here, what is `f.u 1.0`? We cannot differentiate between the multiple
infinitesimal steps at a given real time `t` using only floating-point numbers
as input. Therefore, we keep the list of inputs separated, apart from the
adjacent continuous functions, which we can concatenate together. Until we reach
the horizon given by the input function, we keep alternating continuous and
discrete steps. The output is a list of functions, which we can interpret as a
stream, and the subsystems (the assertions, for instance) can then be called
repeatedly until they reach the same horizon as the parent system. This has the
advantage of not needing/having optional values in the input and output streams,
and the calling system does not need to guess how many steps the simulation will
need to reach the requested horizon.

=== Resets

==== Solver reset

When do we need to reset the solver itself?

If we receive a new input that is not continuous w.r.t. the previous input, we
should reset the solver before starting the integration. This implies that an
input value is tagged with information relating it to the previous input value.
If this is not the case, we can reset whenever we receive a new input value,
although this is less satisfactory.

If a discrete step results in a state jump, we should reset the solver before
continuing the integration.

===== Continuity

When defined in terms of limits, the function $f$ is _continuous at some point_
$c$ of its domain if $ lim_(x -> c) f(x) = f(c) $

The Weierstrass and Jordan definition of continuous function is as follows:
given a function $f : D -> bb(R)$ and an element $x_0 in D$, $f$ is said to be
continuous at the point $x_0$ when the following holds:
$ forall epsilon in bb(R)^star_+. exists delta in bb(R)^star_+. forall x in D.
  x_0 - delta < x < x_0 + delta ==> f(x_0) - epsilon < f(x) < f(x_0) + epsilon $

Every differentiable function $f : (a, b) -> bb(R)$ is continuous. The
derivative $f'(x)$ of a differentiable function $f(x)$ need not be continuous.
If $f'(x)$ is continuous, $f(x)$ is said to be _continuously differentiable_.
The set of such functions is denoted $C^1((a,b))$. More generally, the set of
functions $f : Omega -> bb(R)$ (from an open interval $Omega$ to the reals) such
that $f$ is $n$ times differentiable and such that the $n$-th derivative of $f$
is continuous is denoted $C^n (Omega)$.

Every continuous function $f : [a,b] -> bb(R)$ is integrable.

===== Differentiability classes

Consider an open set $U$ on the real line and a function $f$ defined on $U$ with
real values. Let $k$ be a non-negative integer. The function $f$ is said to be
of differentiability class $C^k$ if the derivatives $f', f'', ..., f^((k))$
exist and are continuous on $U$. If $f$ is $k$-differentiable on $U$, then it is
at least in the class $C^(k-1)$ since $f', f'', ..., f^((k-1))$ are continuous
on $U$. The function $f$ is said to be *infinitely differentiable*, *smooth*, or
of *class* $C^infinity$, if it has derivatives of all orders on $U$.

==== Simulation reset

Two possible options for the simulation reset:
- ask for reset parameters for both the model and solver:
  ```ocaml
  'p * (('y, 'yder) ivp *  ('y, 'zin, 'zout) zc
  ```
- ask for only what is needed for the reset of the model, then build a new
  initial value problem and zero-crossing problem using the new model state and
  reset the solver using them. We can build `fder` and `fzer` quite easily. We
  still need to build `init`, `yc` and `time`.

  `init` and `yc` need to use `model.cget` in order to obtain a continuous state
  element. However, the usual simulation process first starts by executing a
  discrete step to initialize the model state. Does `model.cget` make sense
  before this first step has been executed?
  #figure(finite.automaton(
    (D: (D: "cascade", C: "no cascade"),
     C: (C: "no zero-crossing", D: "zero-crossing")),
    initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
  )) <automaton>
  Tagged values that indicate whether they are discontinuous w.r.t. the previous
  input value do not help us here. Indeed, we could consider that the first
  value in the stream is always considered discontinuous, which would force a
  reset as we receive an input value. However, the initialization of the solver
  still depends on the initialization of the model, which (might) need a
  discrete step, and we are back to square one.

  Additionally, what stop time do we use for the initial state?

  If possible, this would have the advantage of making the reset function only
  dependent on the model, since the solver is updated to match the model's
  simulation requirements.
  ```ocaml
  val sim : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode ->
            ('y, 'yder, 'zin, 'zout) solver ->
            ('p, 'a, 'b) sim
  ```

  Unfortunately, this does not work : the fder and fzer functions are of type
  `fder : 'a -> 'y -> 'yder` and `fzer : 'a -> 'y -> 'zout`, and we thus need a
  way to obtain a value of type `'a`. This is usually done through
  `input.u : time -> 'a`, but we have no input available during the reset, which
  makes this impossible. We thus need reset parameters for both the model and
  solver.

== Mathematical model

#link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
simulation loop as follows:

#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
  The simulation loop of a hybrid system is a _synchronous function_ defining
  four streams $t(n)$, $l x(n)$, $y(n)$ and $z(n)$ with $n in bb(N)$.
  $t(n) in bb(R)$ is the increasing sequence of instants at which the solver
  stops. $l x(n)$ is the value at time $t(n)$ of the _continuous state
  variables_, that is, of all variables defined by their derivatives in the
  original model. $y(n)$ is the value at time $t(n)$ of the _discrete state_.
  $z(n)$ indicates any _zero-crossings_ at instant $t(n)$ on signals monitored
  by the solver, that is, any signals that become equal to or pass through zero.

  The synchronous function has two modes: the discrete mode ($D$) contains all
  computations that may change the discrete state or that have side effects. The
  continuous mode ($C$) is where ODEs are solved. The two modes alternate
  according to the execution scheme summarized in @automaton.

  In the _continuous mode_, the solver computes an approximation of the solution
  of the ODEs and monitors a set of expressions for zero-crossings. The solver
  is abstracted by a function $s o l v e(f)(g)$ parameterized by $f$ and $g$
  where:
  - $x'(tau)=f(y(n),tau,x(tau))$ defines the derivatives of continuous state
    variables $x$ at instant $tau in bb(R)$;
  - $u p z(tau)=g(y(n),tau,x(tau))$ defines the current values of a set of
    zero-crossing signals $u p z$, indexed by $i in {1, ..., k}$.
  The continuous mode $C$ computes four sequences $s$, $l x$, $z$ and $t$ such
  that $ (l x, z, t, s)(n+1) = s o l v e(f)(g)(s,y,l x,t,s t e p)(n) $ where
  - $s(n)$ is the internal state of the solver at instant $t(n) in bb(R)$.
    Calling $s o l v e(f)(g)$ updates the state to $s(n+1)$.
  - $x$ is an approximation of a solution of the ODE
    $ x'(tau)=f(y(n),tau,x(tau)) $
    It is parameterized by the current discrete state $y(n)$ and initialized at
    instant $t(n)$ with the value of $l x(n)$, that is, $x(t(n)) = l x(n)$.
  - $l x(n+1)$ is bounded by the horizon $t(n) + s t e p(n)$ that the solver has
    been asked to reach, that is
    $ t(n) <= t(n+1) <= t(n) + s t e p(n) $
  - $z(n+1)$ signals any zero-crossings detected at time $t(n+1)$. An event
    occurs with a transition to the discrete mode $D$ when horizon
    $t(n) + s t e p(n)$ is reached, or when at least one of the zero-crossing
    signals $u p z(i)$, for $i in {1, ..., k}$ crosses zero, which is indicated
    by a true value for the corresponding boolean output $z(n+1)(i)$. $
    e v e n t = z(n+1)(0) or ... or z(n+1)(k) or (t(n+1) = t(n) + s t e p(n)) $
  If the solver raises an error (for example, a division by zero or an inability
  to find a solution), we consider that the simulation fails.

  All discrete changes occur in the _discrete mode_ $D$. It is entered when an
  event is raised during integration. During a discrete phase, the function
  $n e x t$ defines $y$, $l x$, $s t e p$, $e n c o r e$, $z$ and $t$:
  $ (y, l x, s t e p, e n c o r e)(n+1) & = n e x t(y, l x, z, t)(n) \
                                 z(n+1) & = f a l s e \
                                 t(n+1) & = t(n) $
  where
  - $y(n+1)$ is the new discrete state; outside of mode $D$, $y(n+1) = y(n)$.
  - $l x(n+1)$ is the new continuous state, which may be changed directly in the
    discrete mode.
  - $s t e p(n+1)$ is the new step size.
  - $e n c o r e(n+1)$ is true if an additional discrete step must be performed.
    Function $n e x t$ can decide to trigger instantaneously another discrete
    event, causing an _event cascade_.
  - $t(n)$ (the simulation time) is unchanged during discrete phases.
  The initial values for $y(0)$, $l x(0)$ and $s(0)$ are given by an
  initialization function $i n i t$. Finally, $s o l v e(f)(g)$ may decide to
  reset its internal state if the continuous state changes. If
  $i n i t \_ s o l v e$ initializes the solver state, we have
  $ r e i n i t & = (l x(n+1) != l x(n)) \
         s(n+1) & = "if" r e i n i t "then" i n i t \_ s o l v e(l x(n+1),s(n))
                    "else" s(n) $
  Taken together, the definitions from both modes give a synchronous
  interpretation of the simulation loop as a stream function that computes the
  sequences $l x$, $y$ and $t$ at instant $n+1$ according to their values at
  instant $n$ and an internal state. Writing $s o l v e(f)(g)$ abstracts from
  the actual choice of integration method and zero-crossing detection algorithm.
  A more detailed description of $s o l v e(f)(g)$ would be possible (for
  example, an automaton with two states: one that integrates, and one that
  detects zero-crossings) but with no influence on the code generation problem
  which must be independent of such simulation details.

  Given a program written in a high-level language, we must produce the
  functions $i n i t$, $f$, $g$, and $n e x t$. In practice, they are
  implemented in an imperative language like C. Code generation for hybrid
  models has much in common with code generation for synchronous languages.
])