feat: start of assertions

doc/notes.typ

@@ -129,13 +129,11 @@ cases:
 #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",
+      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")
@@ -180,7 +178,6 @@ interval `[now, now + h]`, and end in one of three possible cases:
 #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,
@@ -319,6 +316,71 @@ 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.
 
+=== Composing simulations
+
+The composition of two simulations must maintain the following invariant: the
+output is only ever absent if the solvers have finished integrating up to the
+requested horizon. This is because we need to keep providing absent values
+(`None`) to the first simulation until it reaches the requested horizon,
+otherwise we reset beforehand, and the only information we have on whether the
+simulation has reached its horizon is the nature of its output. Naïve
+composition breaks this:
+
+#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
+  import cetz.draw: *
+  rect((0, 0), (1, 1))
+  rect((2, 0), (3, 1))
+  content((0.5, 0.5), $M$)
+  content((2.5, 0.5), $N$)
+  line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
+  line((1, 0.5), (2, 0.5), mark: (end: "straight"))
+  line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
+}), 
+$ M_i & : & a_1 & | & bot                     & |
+\ M_o & : & b_1 & | & b_2                     & |
+\ N_i & : & b_1 & | & b_2                     & |
+\ N_o & : & c_1 & | & #text(fill: red, $c_2$) & |
+$)
+
+Instead, we need to make sure we only provide $N$ with a new value when
+needed. This is done through the following composition:
+#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
+  import cetz.draw: *
+  rect((0, 0), (1, 1))
+  rect((2, 0), (3, 1))
+  content((0.5, 0.5), $M$)
+  content((2.5, 0.5), $N$)
+  line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
+  line((1, 0.5), (2, 0.5), mark: (end: "straight"))
+  line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
+}), 
+$ I   &: & a_1 & | & bot & | & bot &   &     & | & bot & | & bot &   &     & |
+\ M_i &: & a_1 & | &     & | &     & | & bot & | &     & | &     & | & bot & |
+\ M_o &: & b_1 & | &     & | &     & | & b_2 & | &     & | &     & | & bot & |
+\ N_i &: & b_1 & | & bot & | & bot & | & b_2 & | & bot & | & bot & | & bot & |
+\ N_o &: & c_1 & | & c_2 & | & bot & | & c_3 & | & c_4 & | & bot & | & bot & |
+\ O   &: & c_1 & | & c_2 & | &     &   & c_3 & | & c_4 & | &     &   & bot & |
+$)
+
+```ocaml
+let step (ms, ns) = function
+  | Some i ->
+      let v, ms = m.step ms (Some i) in
+      let o, ns = n.step ns v in
+      o, (ms, ns)
+  | None ->
+      let o, ns = n.step ns None in
+      match o with
+      | Some o -> Some o, (ms, ns)
+      | None ->
+          let v, ms = m.step ms None in
+          match v with
+          | None -> None, (ms, ns)
+          | Some v ->
+              let o, ns = n.step ns (Some v) in
+              o, (ms, ns)
+```
+
 === Resets
 
 ==== Solver reset
@@ -411,6 +473,67 @@ Two possible options for the simulation reset:
   makes this impossible. We thus need reset parameters for both the model and
   solver.
 
+=== Assertions
+
+Assertions are themselves hybrid systems, with their own internal state, `fder`,
+`fzer`, `fout` and `reset` functions, etc.
+A hybrid model with assertions is then a recursive datatype
+```ocaml
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
+  HNodeA : {
+    body : ('s, 'p,' a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    assertions : ('p, 's, float, 'y, 'yder, 'zin, 'zout) hnode_a list;
+  }
+```
+
+Assertions in sub-models need to be "lifted" to the parent model in order to be
+visible for the simulation. Since the inner state of sub-models is contained in
+the state of the parent model, these assertions can be composed with projectors
+on the inner state of the parent model in order to provide them with the correct
+state.
+
+What does it mean for an assertion to run in continuous time ?
+There are two possible approaches:
+- Sampling: take "snapshots" of the state at various timestamps and check the
+  assertions at these snapshots.
+  - Problems: expensive memory usage, and requires state copies for the model.
+    Can also impact stop times for the parent model, which defeats the point of
+    transparent assertions.
+- Real-time: consider the solution provided by the solver as a continuous
+  function and use it as input to monitor for zero-crossings
+  - Problems: limits the expressions we can use in continuous assertions to
+    continuous functions; boolean expressions are inherently discontinuous. We
+    can model inequalities using zero-crossing constructs like `up`
+    (for instance, `a <= b` $-->$ `up(a -. b)`).
+
+Building a continuous function of the entire state is possible. Since the entire
+state can be considered as the sum of a discrete and a continuous part, and that
+the discrete part is considered piecewise-constant, we can build a function
+`time -> 's` using the model's `cset` function (which updates the state's
+continuous part) as follows:
+#grid(columns: (1fr, 1fr), align: center + horizon, ```ocaml
+let f_st (ms : 's) (f : time -> 'y) =
+  fun t -> cset ms (f t)
+```, cetz.canvas({
+  import cetz-plot.plot: *
+  plot(
+    axis-style: "left", x-label: none, x-tick-step: none,
+    y-label: none, y-tick-step: none, y-min: 0.1, {
+    add(((0, 1), (1, 1)), style: (stroke: blue), label: $D$)
+    add(((1, 2), (2.5, 2)), style: (stroke: blue))
+    add(domain: (0, 1), t => - calc.sin(t + 1) + 1.5, style: (stroke: red))
+    add(domain: (1, 2.5), t => calc.sin(t) + 2, style: (stroke: red), label: $C$)
+  })
+}))
+
+If `cset` is pure (that is, it returns a copy of the state, rather than modify
+it in-place), this function poses no problem. If `cset` is not pure and modifies
+the state in-place, this is still useable as long as `f_state` is never used in
+parallel: since `cset` only modifies the continuous part of the state, we can
+always recover the original state by calling `f_state` with the appropriate
+horizon. However, the model's continuous state needs to be re-updated to the
+reached horizon after all assertions/submodels have finished running.
+
 == Mathematical model
 
 #link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
@@ -497,6 +620,70 @@ simulation loop as follows:
   models has much in common with code generation for synchronous languages.
 ])
 
+In [Branicky, 1995b], dynamical systems are defined as follows: 
+
+#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
+  A continuous (resp. discrete) dynamical system defined on the topological
+  space $X$ over the semigroup $bb(R)^+$ (resp. $bb(Z)^+$) is a function 
+  $f : X times bb(R)^+ -> X$ (resp. $f : X times bb(Z)^+ -> X$) with the
+  following three properties:
+  1. initial condition: $f(p, 0) = p$,
+  2. continuity on both arguments,
+  3. semigroup property: $f(f(p, t_1), t_2) = f(p, t_1 + t_2)$
+  for any point $p in X$ and any $t_1, t_2 in bb(R)^+$ (resp $bb(Z)^+$).
+
+  Technically, such functions are referred to as semidynamical systems, with the
+  term dynamical system reserved for those which the semigroups $bb(R)^+$,
+  $bb(Z)^+$ above may be replaced by the groups $bb(R)$, $bb(Z)$. However, the
+  more "popular" notion of dynamical system in math and engineering - and the
+  one used here - requires only the semigroup property. Thus, the term
+  _reversible_ dynamical system is used when it is necessary to distinguish the
+  group from semigroup case.
+
+  The shorthand $[X, S, f]$ denotes the dynamical system $f$ defined on $X$ over
+  the semigroup $S$; $X$ is referred to as its _state space_ and points in $X$
+  are called _states_. If a dynamical system is defined on a subset of $X$, we
+  say it is a dynamical system _in_ $X$.
+
+  For every fixed value of the parameter $s$, the function $f(dot, s)$ defines
+  a mapping of the space $X$ into itself. Given $[X, bb(Z)^+, f]$, $f(dot, 1)$
+  is its _transition function_. Thus if $[X, bb(Z)^+, f]$ is reversible, its
+  transition function is invertible, with inverse given by $f(dot, -1)$.
+
+  The set $f(p, S) = { f(p, i) | i in S }$ is called the _orbit_ or _trajectory_
+  of the point $p$. A _fixed point_ of $[X, S, f]$ is a point $p$ such that
+  $f(p, s) = p$ for all $s in S$. A set $A subset X$ is _invariant w.r.t._ $f$,
+  or simply _invariant_, if $f(A, s) subset A$ for all $s in S$.
+
+  The notions of equivalence and homomorphism are crucial. Two dynamical systems
+  $[X, S, f]$, $[Y, S, g]$ are said to be _isomorphic_ (also
+  _topologically equivalent_ or simply _equivalent_) if there exists a
+  homeomorphism $psi : X -> Y$ such that
+  $ psi(f(p, s)) = g(psi(p), s) $
+  for all $p in X$ and $s in S$. If the mapping $psi$ is only continuous, then
+  $[X, S, f]$ is said to be _homomorphic_ to $[Y, S, g]$. Homomorphisms preserve
+  trajectories, fixed points, and invariant sets.
+
+  In this paper, the continuous dynamical systems dealt with are defined by the
+  solutions of ordinary differential equations (ODEs):
+  $ x'(t) = f(x(t)) $
+  where $x(t) in X subset bb(R)^n$. The function $f : X -> bb(R)^n$ is called a
+  _vector field_ on $bb(R)^n$. The resulting dynamical system is then given by
+  $phi(x_0, t) = x(t)$ where $x(dot)$ is the solution to the above equation
+  starting at $x_0$ at $t = 0$. (We assume existence and uniqueness of
+  solutions; see [Hirsch, Smalle] for conditions). A system of ODEs is called
+  _autonomous_ or _time-invariant_ if its vector field does not depend
+  explicitly on time. Throughout, the shorthand _continuous_ (resp. _Lipschitz_)
+  _ODEs_ denotes ODEs with continuous (resp. Lipschitz) vector fields.
+
+  An ODE _with input and outputs_ is given by
+  $ x'(t) & = f(x(t), u(t)), \ y(t) & = h(x(t)) $
+  where $x(t) in X subset bb(R)^n$, $u(t) in U subset bb(R)^m$,
+  $y in Y subset bb(R)^p$, $f : bb(R)^n times bb(R)^m -> bb(R)^n$, and
+  $h : bb(R)^n -> bb(R)^p$. The functions $u(dot)$ and $y(dot)$ are the _inputs_
+  and _outputs_, respectively.
+])
+
 == Miscellaneous
 
 === Signals

src/bin/main.ml

@@ -27,7 +27,7 @@ let opts = [
   "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
   "-debug", Arg.Set Debug.debug, "\tPrint debug information";
   "-accelerate", Arg.Set accel, "\tConcatenate continuous functions";
-  "-sundials", Arg.Set sundials, "\tUse sundials (does not support acceleration)";
+  "-sundials", Arg.Set sundials, "\tUse sundials (doesn't support -accelerate)";
   "-inplace", Arg.Set inplace, "\tUse imperative solvers";
   "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
 ]
@@ -46,8 +46,6 @@ let m =
   | Some s        -> Format.eprintf "Unknown model: %s\n" s; exit 2
   with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
 
-let z = StatefulZ.(if !inplace then InPlace.zsolve else Functional.zsolve)
-
 let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
                      else (module State.FunctionalSimState : State.SimState)
 
@@ -55,12 +53,16 @@ let sim =
   if !sundials then
     let open StatefulSundials in
     let c = if !inplace then InPlace.csolve else Functional.csolve in
+    let open StatefulZ in
+    let z = if !inplace then InPlace.zsolve else Functional.zsolve in
     let s = Solver.solver c (d_of_dc z) in
     let open Sim.Sim(val st) in
     run_until_n (Output.print !sample (run m s))
   else
     let open StatefulRK45 in
     let c = if !inplace then InPlace.csolve else Functional.csolve in
+    let open StatefulZ in
+    let z = if !inplace then InPlace.zsolve else Functional.zsolve in
     let s = Solver.solver_c c z in
     let open Sim.Sim(val st) in
     let n = if !accel then accelerate m s else run m (d_of_dc s) in

src/lib/hsim/sim.ml

@@ -26,14 +26,16 @@ module Sim (S : SimState) =
           let mode, stop, now = Continuous, i.h, 0.0 in
           update ms ss (set_running ~mode ~input:i ~stop ~now s)
         end else set_running ~mode:Continuous s in
-      Some { h=0.0; c=Discontinuous; u=fun _ -> o }, s
+      Utils.dot o, s
 
     let step_continuous s step cset fout zset =
       let ms, ss = get_mstate s, get_sstate s in
       let i, now, stop = get_input s, get_now s, get_stop s in
       let (h, f, z), ss = step ss stop in
       let ms = cset ms (f h) in
-      let fout t = fout ms (i.u (now +. t)) (f (now +. t)) in
+      let fy t = f (now +. t) in
+      let fms t = cset ms (fy t) in
+      let fout t = fout ms (i.u (now +. t)) (fy t) in
       let s, c = match z with
       | None ->
           let s, c = if h >= stop
@@ -43,48 +45,102 @@ module Sim (S : SimState) =
       | Some z ->
           let s = set_running ~mode:Discrete ~now:h s in
           update (zset ms z) ss s, Discontinuous in
-      Some { h=h -. now; u=fout; c }, s
+      let h = h -. now in
+      { h; u=fout; c }, s, { h; c; u=fms }
 
     (** Simulation of a model with any solver. *)
     let run
-      (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
-      (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
+      (DNode s : ('y, 'yder, 'zin, 'zout) solver)
-      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
-    = let state = get_init model.state solver.state in
+    = let state = get_init m.state s.state in
-      let step_discrete s =
+      let step_discrete st =
-        step_discrete s model.step model.horizon model.fder model.fzer
+        let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget 
-          model.cget model.csize model.zsize model.jump solver.reset in
+          m.csize m.zsize m.jump s.reset in
+        Some o, s in
+      let step_continuous st =
+        let o, s, _ = step_continuous st s.step m.cset m.fout m.zset in
+        Some o, s in
 
-      let step_continuous s =
+      let step st = function
-        step_continuous s solver.step model.cset model.fout model.zset in
-
-      let step s = function
         | Some i ->
             let mode, now, stop = Discrete, 0.0, i.h in
-            step_discrete (set_running ~mode ~input:i ~now ~stop s)
+            step_discrete (set_running ~mode ~input:i ~now ~stop st)
         | None ->
-            if is_running s then match get_mode s with
+            if is_running st then match get_mode st with
-            | Discrete -> step_discrete s
+            | Discrete -> step_discrete st
-            | Continuous -> step_continuous s
+            | Continuous -> step_continuous st
-            else None, s in
+            else None, st in
 
-      let reset (pm, ps) s =
+      let reset (pm, ps) st =
-        let ms = model.reset pm (get_mstate s) in
+        let ms = m.reset pm (get_mstate st) in
-        let ss = solver.reset ps (get_sstate s) in
+        let ss = s.reset ps (get_sstate st) in
-        update ms ss (set_idle s) in
+        update ms ss (set_idle st) in
+
+      DNode { state; step; reset }
+
+    let rec run_assert :
+      'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a ->
+        (unit -> ('y, 'yder, 'zin, 'zout) solver) ->
+        ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim =
+    fun (HNodeA m) get_s ->
+      let DNode s = get_s () in
+      let al = List.map (fun a -> run_assert a get_s) m.assertions in
+      let state = get_init m.body.state s.state, al in
+
+      let step_discrete (st, al) =
+        let m=m.body in
+        let o, st =
+          step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize
+            m.jump s.reset in
+        let al = List.map (fun (DNode a) ->
+          let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
+          DNode { a with state }) al in
+        Some o, (st, al) in
+        
+      let step_continuous (st, al) =
+        let ({ h; _ } as o), st, u =
+          step_continuous st s.step m.body.cset m.body.fout m.body.zset in
+        let al = List.map (fun (DNode a) ->
+          (* Step assertions repeatedly until they reach the horizon. *)
+          let rec step s =
+            let o, s = a.step s None in
+            match o with None -> s | Some _ -> step s in
+          let state = step (snd @@ a.step a.state (Some u)) in
+          DNode { a with state }) al in
+        (* Reset the model's state to the reached horizon. *)
+        let st = set_mstate (u.u h) st in
+        Some o, (st, al) in
+
+      let step (st, al) = function
+        | Some i ->
+            let mode, now, stop = Discrete, 0.0, i.h in
+            step_discrete (set_running ~mode ~input:i ~now ~stop st, al)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete   -> step_discrete (st, al)
+            | Continuous -> step_continuous (st, al)
+            else None, (st, al) in
+
+      let reset (pm, ps) (st, al) =
+        let ms = m.body.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        update ms ss (set_idle st), al in
 
       DNode { state; step; reset }
 
     let accelerate
       (HNode  m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNodeC s : ('y, 'yder, 'zin, 'zout) solver_c)
-      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
     = let state = get_init m.state s.state in
       let step_discrete st =
-        step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize
+        let o, st = step_discrete st m.step m.horizon m.fder m.fzer m.cget
-          m.jump s.reset in
+          m.csize m.zsize m.jump s.reset in
+        Some o, st in
       let step_continuous st =
-        step_continuous st s.step m.cset m.fout m.zset in
+        let o, st, _ = step_continuous st s.step m.cset m.fout m.zset in
+        o, st in
 
       let rec step st = function
         | Some i ->
@@ -95,14 +151,11 @@ module Sim (S : SimState) =
             | Discrete -> step_discrete st
             | Continuous ->
               let o, st = step_continuous st in
-              match o with
+              match o.c with
-              | None   -> None, st
+              | Discontinuous -> Some o, st
-              | Some { c=Discontinuous; _ } -> o, st
+              | Continuous ->
-              | Some ({ c=Continuous; _ } as o) ->
                   let o', st = step st None in
-                  match o' with
+                  Some (Utils.concat [o; Option.get o']), st
-                  | None    -> assert false
-                  | Some o' -> Some (Utils.concat [o;o']), st
             else None, st in
 
       let reset (pm, ps) st =

src/lib/hsim/types.ml

@@ -6,7 +6,7 @@ type continuity = Continuous | Discontinuous
 type 'a value =
   { h : time;
     u : time -> 'a;                    (* Defined on [[0, h]].                *)
-    c : continuity }
+    c : continuity }                   (* Continuity w.r.t. the next value.   *)
 
 (** A time signal is a sequence of possibly absent α-values
     [{ h; u }] where:
@@ -14,34 +14,26 @@ type 'a value =
   - [u: [0, h] -> α] *)
 type 'a signal = 'a value option
 
+type ('s, 'p, 'a, 'b) drec =
+  { state : 's;
+    step  : 's -> 'a -> 'b * 's;
+    reset : 'p -> 's -> 's }
+
 (** A discrete node. *)
 type ('p, 'a, 'b) dnode =
-  DNode :
+  DNode : ('s, 'p, 'a, 'b) drec -> ('p, 'a, 'b) dnode
+
+type ('s, 'p, 'a, 'b) drec_c =
   { state : 's;
     step  : 's -> 'a -> 'b * 's;
     reset : 'p -> 's -> 's;
-    } -> ('p, 'a, 'b) dnode
+    copy  : 's -> 's }
 
 (** A discrete node which supports a state copy. *)
 type ('p, 'a, 'b) dnode_c =
-  DNodeC :
+  DNodeC : ('s, 'p, 'a, 'b) drec_c -> ('p, 'a, 'b) dnode_c
-    { state : 's;
-      step  : 's -> 'a -> 'b * 's;
-      reset : 'p -> 's -> 's;
-      copy  : 's -> 's;
-    } -> ('p, 'a, 'b) dnode_c
 
-(** A continuous node. *)
+type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
-type ('a, 'b, 'y, 'yder) cnode =
-  CNode :
-    { lsty : 'y;
-      fder : 'a -> 'y -> 'yder;
-      fout : 'a -> 'y -> 'b;
-    } -> ('a, 'b, 'y, 'yder) cnode
-
-(** A hybrid node. *)
-type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
-  HNode :
   { state: 's;
     step : 's -> 'a -> 'b * 's;        (** Discrete step function.            *)
     fder : 's -> 'a -> 'y -> 'yder;    (** Continuous derivative function.    *)
@@ -54,17 +46,32 @@ type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
     cset : 's -> 'y -> 's;             (** Set continuous state.              *)
     zset : 's -> 'zin -> 's;           (** Set zero-crossing state.           *)
     csize : int;
-      zsize : int;
+    zsize : int }
-    } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+
+(** A hybrid node. *)
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
+  HNode : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec ->
+          ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+
+(** A hybrid node with assertions. *)
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
+  HNodeA : {
+    body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list
+  } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
 
 (** The simulation of a hybrid system is a synchronous function on streams of
     functions. *)
-type ('p, 'a, 'b) lazy_sim = ('p, 'a signal, 'b signal) dnode
+type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
 
-(** Greedy simulation takes in an input and computes as many solver and
+(* Utils *)
-    subsystem steps as needed to reach the input's horizon. *)
-type ('p, 'a, 'b) greedy_sim = ('p, 'a value, 'b value list) dnode
-
-(** Utils *)
 
+(** Consider a node with state copying as a node without state copying. *)
 let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset }
+
+(** Consider a model without assertions as a model with an empty list of 
+    assertions. *)
+let a_of_h (HNode body) = HNodeA { body; assertions=[] }
+
+(** Consider a model without its assertions. *)
+let h_of_a (HNodeA { body; _ }) = HNode body

src/lib/hsim/utils.ml

@@ -1,18 +1,13 @@
 
 open Types
 
+let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v }
+
 (** Offset the [u] function by [now]. *)
 let offset (u : time -> 'a) (now : time) : time -> 'a =
   fun t -> u (t +. now)
 
-(**
+(** Concatenate functions. *)
-Concatenate functions. [
-^                    ^
-|     ---,           |     ---,
-|    ___  `---    =  |    _    `---
-| --'                | --'
-+-------------->     +-------------->]
-*)
 let rec concat = function
   | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
   | [f] -> f
@@ -22,6 +17,7 @@ let rec concat = function
       let { h=hr; u=ur; c } = concat l in
       { c; h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) }
 
+(** Sample a function at [n] equidistant points. *)
 let sample { h; u; _ } n =
   let hs = h /. float_of_int n in
   let rec step i =
@@ -44,9 +40,10 @@ let compose (DNode m) (DNode n) =
     (ms, ns) in
   DNode { state; step; reset }
 
-let compose_lazy
+(** Compose two simulations. *)
-  (DNode m : ('p, 'a, 'b) lazy_sim)
+let compose_sim
-  (DNode n : ('q, 'b, 'c) lazy_sim)
+  (DNode m : ('p, 'a, 'b) sim)
+  (DNode n : ('q, 'b, 'c) sim)
 = let state = m.state, n.state in
   let step (ms, ns) = function
     | Some i ->

src/lib/solvers/statefulZ.ml

@@ -58,7 +58,9 @@ module InPlace =
           (h, Some vec), s
         else (h, None), s in
 
-      let copy _ = raise Common.Errors.TODO in
+      let copy s =
+        let vec = zmake (length s.vec) in
+        blit s.vec vec; s.vec <- vec; s in
 
       DNodeC { state; step; reset; copy }
   end