chore: initial commit

.envrc

@@ -0,0 +1 @@
+use flake

.gitignore

@@ -0,0 +1,2 @@
+out
+_build

dune-project

@@ -0,0 +1,16 @@
+(lang dune 3.20)
+
+(name hsim_live)
+
+(source
+ (uri https://git.henri-saudubray.fr/hms/hsim-live))
+
+(authors "Henri Maïga Saudubray <hms@lmf.cnrs.fr>")
+
+(license Unlicense)
+
+(package
+ (name hsim_live)
+ (synopsis "Hybrid system runtime.")
+ (description "Backing project for a live-coding seminar.")
+ (depends ocaml zelus))

exm/ball_discrete.zls

@@ -0,0 +1,21 @@
+
+let dt = 0.001
+let g  = 9.81
+
+let node f_integr (x0, x') = x where
+  rec x = x0 -> pre (x +. x' *. dt)
+let node b_integr (x0, x') = x where
+  rec x = x0 -> (pre x) +. x' *. dt
+
+let node bouncing_ball (p0, v0) = p where
+  rec p = reset f_integr (q, v)    every z
+  and v = reset b_integr (w, -. g) every z
+  and q = p0 -> 0.0 and w = v0 -> -0.8 *. (pre v)
+  and z = false fby (p < 0.0)
+
+let node main () =
+  let rec t = 0.0 fby (dt +. t) in
+  let p = bouncing_ball (5.0, 0.0) in
+  match t <= 10.0 with
+  | true -> (print_float t; print_string "\t"; print_float p; print_newline ())
+  | false -> ()

exm/dune

@@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets exm_ball_discrete.ml ball_discrete.ml ball_discrete.zci)
+ (deps
+  (:zl ball_discrete.zls)
+  (package zelus))
+ (action
+  (run zeluc -s main -o exm_ball_discrete %{zl})))
+
+(executable
+ (name exm_ball_discrete)
+ (public_name exm_ball_discrete)
+ (libraries zelus))

exm/main.zls

@@ -0,0 +1,121 @@
+
+
+
+(** Zélus
+
+    Synchronous language kernel _à la_ Lustre:
+    - programs are Mealy machines (outputs on each transition)
+    - variables represent streams of values in time *)
+
+let node incr x = y where
+  y = x + 1
+
+(*   x │ 8  3  2  7  5  3  …
+    ───┼─────────────────────
+     y │ 9  4  3  8  6  4  …   *)
+
+
+
+
+
+
+(** - we can use values of the previous instants (using [pre]) and
+      initialize streams (using [->]) *)
+
+let node accumulate x = z where
+  rec w = pre x
+  and y = 0 -> pre x
+  and z = x -> (pre z) + x
+
+(*   x │ 1  2  5  2  5  3  …
+    ───┼─────────────────────
+     w │    1  2  5  2  5  …
+     y │ 0  1  2  5  2  5  …
+     z │ 1  3  8 10 15 18  …   *)
+
+
+
+(** - we can reset streams at will *)
+
+let node stay x = y where            (* output the first value forever        *)
+  rec y = x -> pre y
+
+let node from x = y where            (* count up from the first value         *)
+  rec y = x -> pre y + 1
+
+let node loop x = y where
+  rec y = reset from 0 every z
+  and z = false -> pre y >= w
+  and w = stay x
+
+(*        x │ 2  _  _  _  _  _  …
+    ────────┼─────────────────────
+     loop x │ 0  1  2  0  1  2  …   *)
+
+
+(** Already able to model physical behaviours! *)
+
+let dt = 0.001                        (* Integration step                     *)
+let g  = 9.81                         (* Gravitational constant               *)
+let node f_integr (x0, x') = x where  (* Forward Euler integrator             *)
+  rec x = x0 -> pre (x +. x' *. dt)
+let node b_integr (x0, x') = x where  (* Backward Euler integrator            *)
+  rec x = x0 -> (pre x) +. x' *. dt
+
+let node bouncing_ball (p0, v0) = p where
+  rec p = reset f_integr (q, v)    every z
+  and v = reset b_integr (w, -. g) every z
+  and q = p0 -> 0.0 and w = v0 -> -0.8 *. (pre v)
+  and z = false -> (pre p) < 0.0
+
+(** Quite cumbersome. *)
+
+
+
+
+
+
+(** Enter continuous-time constructs:
+    - express values with initial value problems with [der] and [init] *)
+
+let hybrid integr (x0, x') = x where
+  der x = x' init x0
+
+let hybrid position (p0, v0, a) = p where
+  rec der p = v init p0
+  and der v = a init v0
+
+
+
+
+
+
+
+
+
+
+
+
+
+(** We can intermingle discrete and continuous behaviours: *)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+(** We can now express physical systems much more precisely: *)
+
+let hybrid bouncing_ball (p0, v0) = p where
+  rec der p = v    init p0 reset z -> 0.0
+  and der v = -. g init v0 reset z -> -0.8 *. last v
+  and     z = up(-. p)

flake.lock

@@ -0,0 +1,27 @@
+{
+  "nodes": {
+    "nixpkgs": {
+      "locked": {
+        "lastModified": 1771848320,
+        "narHash": "sha256-0MAd+0mun3K/Ns8JATeHT1sX28faLII5hVLq0L3BdZU=",
+        "owner": "nixos",
+        "repo": "nixpkgs",
+        "rev": "2fc6539b481e1d2569f25f8799236694180c0993",
+        "type": "github"
+      },
+      "original": {
+        "owner": "nixos",
+        "ref": "nixos-unstable",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
+    "root": {
+      "inputs": {
+        "nixpkgs": "nixpkgs"
+      }
+    }
+  },
+  "root": "root",
+  "version": 7
+}

flake.nix

@@ -0,0 +1,29 @@
+{
+  description = "Live-coding of a hybrid language runtime.";
+
+  inputs.nixpkgs.url = "github:nixos/nixpkgs?ref=nixos-unstable";
+
+  outputs =
+    { nixpkgs, ... }:
+    let
+      system = "x86_64-linux";
+      pkgs = import nixpkgs {
+        inherit system;
+        config.allowUnfree = true;
+      };
+    in
+    {
+      devShells."${system}".default = pkgs.mkShell {
+        packages =
+          with pkgs;
+          [ feedgnuplot ]
+          ++ (with ocamlPackages; [
+            findlib
+            ocaml
+            ocaml-lsp
+            dune_3
+            zelus
+          ]);
+      };
+    };
+}

justfile

@@ -0,0 +1,3 @@
+exm1:
+    timeout 1s dune exec exm_ball_discrete \
+    | feedgnuplot --stream --domain --lines

lib/hsim/dune

@@ -0,0 +1,7 @@
+(env
+ (dev
+  (flags
+   (:standard -w -50))))
+
+(library
+ (name hsim))

lib/hsim/fill.ml

@@ -0,0 +1,270 @@
+[@@@warning "-27-50"]
+let todo = assert false
+
+
+
+
+(* Little OCaml reminder: *)
+type t = { a : int; b : int; c : int }
+
+let () =
+  let x = { a = 0; b = 1; c = 2 } in
+  let y = { x with a = 2 } in
+  assert (y = { a = 2; b = 1; c = 2 })
+
+
+
+
+
+
+
+
+(** Discrete-time node *)
+type ('i, 'o, 'r) dnode =
+  DNode : {
+    state : 's;                        (** current state                      *)
+    step  : 's -> 'i -> 's * 'o;       (** step function                      *)
+    reset : 's -> 'r -> 's;            (** reset function                     *)
+  } -> ('i, 'o, 'r) dnode
+
+
+(** Run a discrete node on a list of inputs *)
+let drun (DNode n : ('i, 'o, 'r) dnode) (i : 'i list) : 'o list =
+  todo
+
+
+
+
+
+type time =
+  float                                (** [≥ 0.0]                            *)
+
+
+(** Interval-defined functions *)
+type 'a dense =
+  { h : time;                          (** horizon                            *)
+    f : time -> 'a }                   (** [f : [0, h] -> α]                  *)
+
+
+(** Continuous-time signal *)
+type 'a signal =
+  'a dense option
+
+
+
+
+
+
+
+
+
+(** Initial value problem (IVP) *)
+type ('y, 'yder) ivp =
+  { y0   : 'y;                         (** initial position                   *)
+    fder : time -> 'y -> 'yder;        (** derivative function on [[0.0, h]]  *)
+    h    : time; }                     (** maximal horizon                    *) 
+
+
+
+
+
+
+
+
+
+
+
+
+
+(** ODE solver *)
+type ('y, 'yder) csolver =
+  (time,                               (** requested horizon                  *) 
+   'y dense,                           (** solution approximation             *) 
+   ('y, 'yder) ivp)                    (** initial value problem              *) 
+  dnode
+
+
+
+
+
+
+
+
+
+
+
+
+(** Zero-crossing problem (ZCP) *)
+type ('y, 'zin) zcp =
+  { y0   : 'y;                         (** initial position                   *)
+    fzer : time -> 'y -> 'zin;         (** zero-crossing function             *)
+    h    : time; }                     (** maximal horizon                    *)
+
+
+
+
+
+
+
+
+
+
+
+
+
+(** Zero-crossing solver *)
+type ('y, 'zin, 'zout) zsolver =
+  ('y dense,                           (** input value                        *)
+   time * 'zout,                       (** horizon and zero-crossing events   *)
+   ('y, 'zin) zcp)                     (** zero-crossing problem              *)
+  dnode
+
+
+
+
+
+
+
+
+
+
+
+
+(** Full solver (composition of an ODE and zero-crossing solver) *)
+type ('y, 'yder, 'zin, 'zout) solver =
+  (time,                               (** requested horizon                  *)
+   'y dense * 'zout,                   (** output and zero-crossing events    *)
+   ('y, 'yder) ivp * ('y, 'zin) zcp)   (** (re)initialization parameters      *)
+  dnode
+
+
+
+
+
+
+
+
+(** Compose an ODE solver and a zero-crossing solver *)
+let build_solver : ('y, 'yder) csolver ->
+                   ('y, 'zin, 'zout) zsolver ->
+                   ('y, 'yder, 'zin, 'zout) solver
+= fun (DNode cs) (DNode zs) ->
+  let state = (cs.state, zs.state) in
+  let step (cstate, zstate) (h : time) =
+    todo in
+
+
+  let reset (cstate, zstate) (ivp, zcp) =
+    (cs.reset cstate ivp, zs.reset zstate zcp) in
+  DNode { state; step; reset }
+
+
+
+
+
+(** Hybrid (discrete-time and continuous-time) node *)
+type ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode =
+  HNode : {
+    state : 's;                        (** current state                      *)
+    step  : 's -> 'i -> 's * 'o;       (** discrete step function             *)
+    reset : 's -> 'r -> 's;            (** reset function                     *)
+    fder  : 's -> 'i -> 'y -> 'yder;   (** derivative function                *)
+    fzer  : 's -> 'i -> 'y -> 'zin;    (** zero-crossing function             *)
+    fout  : 's -> 'i -> 'y -> 'o;      (** continuous output function         *)
+    cget  : 's -> 'y;                  (** continuous state getter            *)
+    cset  : 's -> 'y -> 's;            (** continuous state setter            *)
+    zset  : 's -> 'zout -> 's;         (** zero-crossing information setter   *)
+    jump  : 's -> bool;                (** discrete go-again indicator        *)
+  } -> ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode
+
+
+
+
+(** Simulation mode (either discrete or continuous) *)
+type mode = D | C
+
+(** Simulation state *)
+type ('i, 'o, 'r, 'y) state =
+  State : {
+    solver :                           (** solver state                       *)
+      ('y, 'yder, 'zin, 'zout) solver;
+    model :                            (** model state                        *)
+      ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode;
+    input : 'i signal;                 (** current input                      *)
+    time : time;                       (** current time                       *)
+    mode : mode;                       (** current step mode                  *)
+  } -> ('i, 'o, 'r, 'y) state
+
+
+
+(** Discrete simulation step *)
+let dstep (State ({ model = HNode m; solver = DNode s; _ } as state)) =
+  let i = Option.get state.input in
+  let (ms, o) = m.step m.state (todo (* current input? *)) in
+  let model = HNode { m with state = todo (* ? *) } in
+  let state =
+    if m.jump ms then State { state with model = todo (* ? *) }
+    else if state.time >= i.h then
+      State { state with input = todo (* ? *); model; time = todo (* ? *) }
+    else
+      let y0 = todo (* ? *) and h = i.h -. state.time in
+      let ivp = { h; y0; fder = fun t y -> m.fder ms (i.f todo (* ? *)) y } in
+      let zcp = { h; y0; fzer = fun t y -> m.fzer ms (i.f todo (* ? *)) y } in
+      let solver = DNode { s with state = s.reset s.state (ivp, zcp) } in
+      State { state with model; solver; mode = todo (* ? *) } in
+  (state, Some { h = 0.; f = fun _ -> o })
+
+
+  
+
+
+(** Continuous simulation step *)
+let cstep (State ({ model = HNode m; solver = DNode s; _ } as st)) =
+  let i = Option.get st.input in
+  let (ss, (y, z)) = s.step s.state todo (* ? *) in 
+  let ms = m.zset (m.cset m.state (y.f y.h)) z in
+  let out = Some { y with f = fun t -> m.fout ms todo (* ? *) (y.f t) } in
+  let mode = if m.jump ms || st.time +. y.h >= i.h then D else C in
+  let model = HNode { m with state = ms } in
+  let solver = DNode { s with state = ss } in
+  (State { st with model; solver; mode; time = todo (* ? *) }, out)
+
+
+
+
+
+(** Simulate a hybrid model with a solver *)
+let hsim : ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode ->
+           ('y, 'yder, 'zin, 'zout) solver ->
+           ('i signal, 'o signal, 'r) dnode
+= fun model solver ->
+  let state = State { model; solver; input = None; time = 0.; mode = D } in
+  let step (State s as st) input = match (input, s.input, s.mode) with
+    | Some _, None,   _ -> dstep (State { s with input; time = 0.; mode = D })
+    | None,   Some _, D -> dstep st
+    | None,   Some _, C -> cstep st
+    | None,   None,   _ -> (st, None)
+    | Some _, Some _, _ -> invalid_arg "Not done processing previous input" in
+  let reset (State ({ model = HNode m; _ } as s)) r =
+    let model = HNode { m with state = m.reset m.state r } in
+    State { s with model; input = None; time = 0.; mode = D } in
+  DNode { state; step; reset }
+
+
+
+
+(** Run a simulation on a list of inputs *)
+let hrun (model : ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode)
+         (solver : ('y, 'yder, 'zin, 'zout) solver)
+         (i : 'i dense list) : 'o dense list
+= let sim = hsim model solver and i = List.map Option.some i in
+  let rec step os (DNode sim) i =
+    let state, o = sim.step sim.state i in
+    let sim = DNode { sim with state } in
+    if o = None then (sim, List.rev_map Option.get os)
+    else step (o :: os) sim None in
+  List.fold_left_map (step []) sim i |> snd |> List.flatten
+
+
+
+

lib/hsim/full.ml

@@ -0,0 +1,280 @@
+[@@@warning "-27-50"]
+
+(** Discrete-time node
+
+    Low-level representation of a discrete node as a current [state]
+    and [step] and [reset] functions. The type parameters represent
+    the step function's inputs and outputs and the reset parameter. *)
+type ('i, 'o, 'r) dnode =
+  DNode : {
+    state : 's;                        (** current state                      *)
+    step  : 's -> 'i -> 's * 'o;       (** step function                      *)
+    reset : 's -> 'r -> 's;            (** reset function                     *)
+  } -> ('i, 'o, 'r) dnode
+
+
+(** Run a discrete node on a list of inputs *)
+let drun (DNode n : ('i, 'o, 'r) dnode) (i : 'i list) : 'o list =
+  List.fold_left_map n.step n.state i |> snd
+
+
+(** Time values
+
+    Must be positive. *)
+type time = float
+
+
+(** Interval-defined functions
+
+    A value [v] of type [α dense] represents a function from [time] to
+    [α] defined on the interval [[0, v.h]]. Calling [v.f] with a value
+    outside these bounds is undefined. For convenience, we assume all
+    functions are well-behaved w.r.t. the numerical methods we use.
+
+    In particular, if [v.h = 0], then the function is defined at a
+    single instant. *)
+type 'a dense =
+  { h : time;                          (** horizon                            *)
+    f : time -> 'a }                   (** [f : [0, h] -> α]                  *)
+
+
+(** Continuous-time signal
+
+    A value of type [α signal] is either empty or an interval-defined
+    function. A sequence of values of type [α signal] represents the
+    function obtained by the concatenation of the interval domains of
+    each [α value] in the sequence.
+
+    For instance, [[Some { h=3.0; f=f1 }; None; Some { h=2.0; f=f2 }]]
+    represents the function
+    [{ h=5.0; f=fun t -> if t <= 3.0 then f1 t else f2 (t - 3.0) }]. *)
+type 'a signal =
+  'a dense option
+
+
+(** Initial value problem (IVP)
+
+    Its solution is a function [f : [0, h] -> 'y] such that:
+    - [f 0.0 = y0] 
+    - [(df/dt) t = fder t (f t)] *)
+type ('y, 'yder) ivp =
+  { y0   : 'y;                         (** initial position                   *)
+    fder : time -> 'y -> 'yder;        (** derivative function                *)
+    h    : time; }                     (** maximal horizon                    *) 
+
+
+(** ODE solver
+
+    Given a requested horizon [t], the solver returns an approximation
+    of the solution to the IVP on [[0, t']] (where [t' ≤ t]).
+    Successive steps compute successive parts of the solution.
+
+    Its (re)initialization parameter is the IVP to solve. That is, the
+    solver must be initialized with an IVP before use. *)
+type ('y, 'yder) csolver =
+  (time, 'y dense, ('y, 'yder) ivp) dnode
+(*  ↑       ↑            ↑         *)
+(* input  output   reset parameter *)
+
+
+(** Zero-crossing problem (ZCP)
+
+    Paired with an approximation [f : [0, h] -> 'y], its solution is
+    the least instant [t ∈ [0, h]] such that [fzer t (y t)] crosses
+    zero at that instant, or [h] if no such crossing occurs. *)
+type ('y, 'zin) zcp =
+  { y0   : 'y;                         (** initial position                   *)
+    fzer : time -> 'y -> 'zin;         (** zero-crossing function             *)
+    h    : time; }                     (** maximal horizon                    *)
+
+
+(** Zero-crossing solver
+
+    Given an approximation [f : [0, h] -> 'y], the solver returns an
+    instant [t ∈ [0, h]] solving the ZCP, and an indicator ['zout] of
+    what zero-crossing event occured. *)
+type ('y, 'zin, 'zout) zsolver =
+  ('y dense, time * 'zout, ('y, 'zin) zcp) dnode
+(*    ↑           ↑              ↑         *)
+(*  input       output     reset parameter *)
+
+
+(** Full solver
+
+    Composes an ODE solver with a zero-crossing solver. *)
+type ('y, 'yder, 'zin, 'zout) solver =
+  (time, 'y dense * 'zout, ('y, 'yder) ivp * ('y, 'zin) zcp) dnode
+(*  ↑           ↑                         ↑        *)
+(* input      output               reset parameter *)
+
+
+(** Compose an ODE solver and a zero-crossing solver. *)
+let compose_solvers :
+    ('y, 'yder) csolver ->
+    ('y, 'zin, 'zout) zsolver ->
+    ('y, 'yder, 'zin, 'zout) solver =
+  fun (DNode csolver) (DNode zsolver) ->
+    let state = (csolver.state, zsolver.state) in
+    let step (cstate, zstate) h =
+      let cstate, f = csolver.step cstate h in
+      let zstate, (h, zout) = zsolver.step zstate f in
+      (cstate, zstate), ({ f with h }, zout) in
+    let reset (cstate, zstate) (ivp, zcp) =
+      (csolver.reset cstate ivp, zsolver.reset zstate zcp) in
+    DNode { state; step; reset }
+
+
+(** Hybrid (discrete-time and continuous-time) node
+
+    A hybrid node contains both a discrete [step] function and a
+    derivative function [fder], zero-crossing function [fzer], and
+    output function [fout], representing continuous behaviour.
+
+    Its state ['s] contains a continuous part ['y], which evolves
+    during continuous behaviour; functions [cget] and [cset] allow
+    reading of and writing to this continuous part.
+
+    The zero-crossing function [fzer] returns a vector ['zin]
+    of all the expressions to monitor for zero-crossings. When a
+    zero-crossing event is detected, the state may be updated using
+    the [zset] function.
+
+    After a discrete step, the model may require another discrete step
+    to be performed before resuming continuous behaviour, which it
+    notifies through the [jump] function. *)
+type ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode =
+  HNode : {
+    state : 's;                        (** current state                      *)
+    step  : 's -> 'i -> 's * 'o;       (** discrete step function             *)
+    reset : 's -> 'r -> 's;            (** reset function                     *)
+    fder  : 's -> 'i -> 'y -> 'yder;   (** derivative function                *)
+    fzer  : 's -> 'i -> 'y -> 'zin;    (** zero-crossing function             *)
+    fout  : 's -> 'i -> 'y -> 'o;      (** continuous output function         *)
+    cget  : 's -> 'y;                  (** continuous state getter            *)
+    cset  : 's -> 'y -> 's;            (** continuous state setter            *)
+    zset  : 's -> 'zout -> 's;         (** zero-crossing information setter   *)
+    jump  : 's -> bool;                (** discrete go-again function         *)
+  } -> ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode
+
+
+(* We consider the simulation of a hybrid node with a solver as a
+   particular case of a discrete node. That is, the simulation has an
+   internal state, step function and reset function. At each step of
+   the simulation, we operate according to one of two modes: discrete
+   and continuous.
+
+   In discrete mode, we perform a discrete step using the model's
+   [step] function. Physical time does not advance; and so the output
+   is a function defined at a single instant [0.0]. Additionally, we
+   may choose to reset the solver and switch to continuous mode in the
+   next step, according to the result of the model's [jump] function.
+
+   In continuous mode, we call the solver to obtain an approximation
+   of the solution to the model's IVP, obtained with its [fder]
+   function. Physical time advances up to the horizon reached by the
+   solver. If a zero-crossing event occurs, we switch to discrete mode
+   for the next step; otherwise we remain in continuous mode. *)
+
+
+(** Simulation mode
+
+    Either discrete ([D]) or continuous ([C]). *)
+type mode = D | C
+
+
+(** Simulation state
+
+    The simulation state must store the states of both the model and
+    solver. Additionally, it must store the current input, physical
+    time, and step mode. *)
+type ('i, 'o, 'r, 'y) state =
+  State : {
+    solver :                           (** solver state                       *)
+      ('y, 'yder, 'zin, 'zout) solver;
+    model :                            (** model state                        *)
+      ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode;
+    input : 'i signal;                 (** current input                      *)
+    time : time;                       (** current time                       *)
+    mode : mode;                       (** current step mode                  *)
+  } -> ('i, 'o, 'r, 'y) state
+
+
+(** Discrete simulation step
+
+    Performs a discrete step of the model, and resets the solver if
+    required by the model. *)
+let dstep (State ({ model = HNode m; solver = DNode s; _ } as st)) =
+  let i = Option.get st.input in
+  (* Step the model *)
+  let ms, o = m.step m.state (i.f st.time) in
+  let model = HNode { m with state = ms } in
+  let st =
+    if m.jump ms then
+      (* If the model asks for another discrete step, stay in discrete mode; *)
+      State { st with model }
+    else if st.time >= i.h then
+      (* if we reached the end of the input, remain idle; *)
+      State { st with input = None; model; time = 0. }
+    else
+      (* otherwise, reset the solver and switch to continuous mode. *)
+      let y0 = m.cget ms and h = i.h -. st.time in
+      let ivp = { h; y0; fder = fun t -> m.fder ms (i.f (t +. st.time)) } in
+      let zcp = { h; y0; fzer = fun t -> m.fzer ms (i.f (t +. st.time)) } in
+      let solver = DNode { s with state = s.reset s.state (ivp, zcp) } in
+      State { st with model; solver; mode = C } in
+  (* Return the state and the output function (defined only at [0.0]). *)
+  st, Some { h = 0.; f = fun _ -> o }
+  
+(** Continuous simulation step
+
+    Calls the solver to solve the IVP, and switch to discrete mode if
+    a zero-crossing event occurs or if the model asks for it. *)
+let cstep (State ({ model = HNode m; solver = DNode s; _ } as st)) =
+  let i = Option.get st.input in
+  (* Step the solver. *)
+  let ss, (y, z) = s.step s.state i.h in 
+  (* Update the model's state. *)
+  let ms = m.zset (m.cset m.state (y.f y.h)) z in
+  (* Create the output function. *)
+  let out = { y with f = fun t -> m.fout ms (i.f (t +. st.time)) (y.f t) } in
+  (* Compute the mode for the next step. *)
+  let mode = if m.jump ms || st.time +. y.h >= i.h then D else C in
+  let model = HNode { m with state = ms } in
+  let solver = DNode { s with state = ss } in
+  (* Return the state and the output function. *)
+  State { st with model; solver; mode; time = st.time +. y.h }, Some out
+
+(** Simulate a hybrid model with a solver
+
+    The [step] function calls [dstep] or [cstep] depending on the
+    current mode. If the current input is [None], nothing is done; the
+    simulation is awaiting input; and we are allowed to provide a new
+    function as input. If the current input is [Some f], the [step]
+    function expects [None] as input; providing a new input value
+    before the simulation is done with the previous one is an error. *)
+let hsim : ('i, 'o, 'r, 'y, 'yder, 'zin, 'zout) hnode ->
+           ('y, 'yder, 'zin, 'zout) solver ->
+           ('i signal, 'o signal, 'r) dnode =
+  fun model solver ->
+    let state = State { model; solver; input = None; time = 0.; mode = D } in
+    let step (State s as st) input = match (input, s.input, s.mode) with
+      | (Some _, None, _) ->
+          (* Accept the new input and reset the solver. *)
+          dstep (State { s with input; time = 0.; mode = D })
+      | (None, Some _, D) ->
+          (* Perform a discrete step on the current input. *)
+          dstep st
+      | None, Some _, C ->
+          (* Perform a continuous step on the current input. *)
+          cstep st
+      | (None, None, _) ->
+          (* Do nothing and wait for the next input. *)
+          (st, None)
+      | (Some _, Some _, _) ->
+          (* Got the next input too early! *)
+          invalid_arg "Not done processing previous input" in
+    let reset (State ({ model = HNode m; _ } as s)) r =
+      (* Reset the model; the solver will reset at the first discrete step. *)
+      let model = HNode { m with state = m.reset m.state r } in
+      State { s with model; input = None; time = 0.; mode = D } in
+    DNode { state; step; reset }