commit a41e6b2faa19f05c5b70ce5ba42d5dd613eb2b17 Author: Henri Saudubray Date: Fri Mar 27 10:53:26 2026 +0100 chore: initial commit diff --git a/.envrc b/.envrc new file mode 100644 index 0000000..3550a30 --- /dev/null +++ b/.envrc @@ -0,0 +1 @@ +use flake diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8cb11e0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +out +_build diff --git a/dune-project b/dune-project new file mode 100644 index 0000000..f4f9ed8 --- /dev/null +++ b/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 ") + +(license Unlicense) + +(package + (name hsim_live) + (synopsis "Hybrid system runtime.") + (description "Backing project for a live-coding seminar.") + (depends ocaml zelus)) diff --git a/exm/ball_discrete.zls b/exm/ball_discrete.zls new file mode 100644 index 0000000..d8fca3b --- /dev/null +++ b/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 -> () diff --git a/exm/dune b/exm/dune new file mode 100644 index 0000000..c271fce --- /dev/null +++ b/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)) diff --git a/exm/main.zls b/exm/main.zls new file mode 100644 index 0000000..a82cd08 --- /dev/null +++ b/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) diff --git a/flake.lock b/flake.lock new file mode 100644 index 0000000..82d8ab7 --- /dev/null +++ b/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 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..5fae668 --- /dev/null +++ b/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 + ]); + }; + }; +} diff --git a/justfile b/justfile new file mode 100644 index 0000000..431ddbd --- /dev/null +++ b/justfile @@ -0,0 +1,3 @@ +exm1: + timeout 1s dune exec exm_ball_discrete \ + | feedgnuplot --stream --domain --lines diff --git a/lib/hsim/dune b/lib/hsim/dune new file mode 100644 index 0000000..50c2189 --- /dev/null +++ b/lib/hsim/dune @@ -0,0 +1,7 @@ +(env + (dev + (flags + (:standard -w -50)))) + +(library + (name hsim)) diff --git a/lib/hsim/fill.ml b/lib/hsim/fill.ml new file mode 100644 index 0000000..3513baf --- /dev/null +++ b/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 + + + + diff --git a/lib/hsim/full.ml b/lib/hsim/full.ml new file mode 100644 index 0000000..e2e01d3 --- /dev/null +++ b/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 }