feat (doc): some notes

doc/notes.typ

@@ -0,0 +1,97 @@
+#import "@preview/cetz:0.3.2"
+#import "@preview/cetz-plot:0.1.1"
+
+#set page(paper: "a4")
+#set par(justify: true)
+
+= Notes on hybrid system simulation
+
+== The problem
+
+The initial problem is _transparent assertions_. We want to simulate assertions 
+with their own solvers, because simulating assertions with the same solver as 
+the 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 :
+  { s: '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) 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, '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.
+```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.
+
+== Output format
+
+In "lazy" mode, the output has format ```ocaml 'b signal = 'b value option```.
+
+In "greedy" mode, the output has format ```ocaml 'b value list```.
+
+Indeed, 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:
+
+#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({
+    cetz-plot.plot.plot(
+      size: (1, 1), x-tick-step: 1, y-tick-step: 1, axis-style: "left", {
+        cetz-plot.plot.add(((0,0),(1,1)))
+        cetz-plot.plot.add(((1,2),), mark: "o", mark-size: .03)})}))
+  #colbreak()
+  #linebreak()
+  $ f(t) = cases(
+    t   & "if" 0.0 <= t <= 1.0,
+    2.0 & "if" t = 1.0 + partial,
+    bot & "otherwise",
+  ) $
+])
+
+How do we represent infinitesimal steps in floating-point numbers ?
+
+A possible representation of the above could be
+#align(center, ```ocaml
+f = [{ start = 0.0; length = 1.0; u = fun t -> t   };
+     { start = 1.0; length = 0.0; u = fun _ -> 2.0 }]
+```)
+but returning a list is problematic, however, as it cannot be passed easily as 
+input to the underlying models. Unless the input can itself be considered as a 
+list, in which case the simulation simply operates over the whole list, 
+resetting the solver as needed. This has the advantage of not having to wait 
+several steps, providing nothing, until the solver is done integrating the 
+input, before we can provide it with another input.

src/lib/hsim/sim.ml

@@ -60,6 +60,10 @@ module LazySim (S : SimState) =
                 let (h, f, z), sstate = solver.step sstate stop in
                 let mstate = model.cset mstate (f h) in
                 let h' = input.start +. h in
+                let fout t =
+                  model.fout mstate (input.u (now +. t)) (f (now +. t)) in
+                let out =
+                  { start = input.start +. now; length = h -. now; u = fout } in
                 let state = match z with
                 | None ->
                     let status =
@@ -70,11 +74,6 @@ module LazySim (S : SimState) =
                     let status = S.running ~mode:Discrete ~now:h' status in
                     let mstate = model.zset mstate z in
                     S.update ~status ~mstate ~sstate state in
-                let fout t =
-                  model.fout mstate (input.u (now +. t)) (f (now +. t)) in
-                let out =
-                  { start = input.start +. now; length = h -. now; u = fout }
-                in
                 Some out, state in
       let reset m s =
         let mstate = model.reset m (S.mstate s) in
@@ -92,20 +91,59 @@ module LazySim (S : SimState) =
 module GreedySim (S : SimState) =
   struct
 
+    (* TODO: greedy simulation *)
+
     let sim
       (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
       : ('p, 'a, 'b) sim
     = let state = S.init ~mstate:model.state ~sstate:solver.state in
       let rec step state input =
+        let status = S.status state and mstate = S.mstate state
+        and sstate = S.sstate state in
         match input, S.is_running state with
         | Some input, _ ->
             let mode = Discrete and now = 0.0 and stop = input.length in
             let status = S.running ~mode ~input ~now ~stop (S.status state) in
             let state = S.update ~status state in
-            None, state
+            step state None
         | None, false -> None, state
-        | None, true -> assert false
+        | None, true ->
+            let input = S.input state and now = S.now state
+            and stop = S.stop state in
+            match S.mode state with
+            | Discrete ->
+                let o, mstate = model.step mstate (input.u now) in
+                let state =
+                  let h = model.horizon mstate in
+                  if h <= 0.0 then S.update ~mstate state
+                  else if now >= stop then raise Common.Utils.TODO
+                  else if model.jump mstate then
+                    let y = model.cget mstate in
+                    let fder t = model.fder mstate (offset input now t) in
+                    let fzer t = model.fzer mstate (offset input now t) in
+                    let ivp = { fder; stop = stop -. now; init = y } in
+                    let zc = { yc = y; fzer } in
+                    let sstate = solver.reset (ivp, zc) sstate in
+                    let status = S.running ~mode:Continuous status in
+                    S.update ~status ~mstate ~sstate state
+                  else
+                    let status = S.running ~mode:Continuous status in
+                    S.update ~status state in
+                let start = input.start +. now in
+                Some { start; length = 0.0; u = fun _ -> o }, state
+            | Continuous ->
+                let (h, f, z), sstate = solver.step sstate stop in
+                let mstate = model.cset mstate (f h) in
+                let h' = input.start +. h in
+                let fout t =
+                  model.fout mstate (input.u (now +. t)) (f (now +. t)) in
+                let out =
+                  { start = input.start +. now; length = h -. now; u = fout } in
+                match z with
+                | None ->
+                    let status =
+                      if h >= stop then S.running ~mode:Discrete ~now:h' status
       in
       let reset = assert false in
       DNode { state; step; reset }