feat: lift runtime into language, start of zelus 2024 compatibility

doc/notes.typ

@@ -751,3 +751,20 @@ let hybrid ball () = y where
   and assert (y >= 0)
 ```
 Is this an issue?
+
+= Ideas
+
+== Continuous assertions with lifted runtime
+
+Given a synchronous node simulating a hybrid one :
+```zelus
+let hybrid sincos () = (sin, cos) where
+  rec der sin = cos init 0.0
+  and der cos = -. sin init 1.0
+
+let node sincos_sim = Solve.solve_ode45 sincos
+```
+We could create a primitive `Solve.assert_continuous` which would take as input
+an `'a value` and a function `'a -> Solve.cond`, where `Solve.cond` is a
+zero-crossing expression, for instance `fun v -> Solve.up(-. v)`. This could be
+passed along to a zero-crossing solver during continuous phases.

doc/rep.typ

@@ -11,10 +11,12 @@
 #let simulink = smallcaps[Simulink]
 #let sundials = smallcaps[Sundials CVODE]
 #let zelus = smallcaps[Zélus]
-#let TODO(..what) = {
+#let note(color, prefix, ..what) = {
   let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") }
-  block(fill: red, width: 100%, inset: 5pt, align(center, raw("TODO" + msg)))
+  block(fill: color, width: 100%, inset: 5pt, align(center, raw(prefix + msg)))
 }
+#let TODO(..what) = note(red, "TODO", ..what)
+#let MENTION(..what) = note(gray, "MENTION", ..what)
 #let adot(s) = $accent(#s, dot)$
 #let addot(s) = $accent(#s, dot.double)$
 
@@ -97,14 +99,14 @@ physical systems. Continuous phases are described using ordinary differential
 equations (ODEs), while discrete phases can be represented as a reactive
 program in a synchronous language such as #lustre or #esterel.
 
-As a first example, say we wish to model a bouncing ball. We could start by
+As an illustration, say we wished to model an extensively studied system: a
-describing its distance from the ground $y$ with a second-order differential
+bouncing ball. We could start by describing its distance from the ground $y$ as
-equation
+a function of time, with a second-order differential equation
-$ addot(y) = -9.81 $
+$ addot(y) = -9.81, $
 where $addot(y)$ denotes the second order derivative of $y$ with
-respect to time (the acceleration of the ball), and $9.81$ is the gravitational
+respect to time $(d^2y)/(d t^2)$ (the acceleration of the ball), and $9.81$ is
-constant $g$: the acceleration of the ball is its negation. We now give the
+the gravitational constant $g$: the acceleration of the ball is its negation. We
-initial position and speed of the ball:
+now give the initial position and speed of the ball:
 $ y(0) = 50 space space space adot(y)(0) = 0 $
 We have just described an initial value problem: given an ODE and an initial
 value for its dependent variable, its solution is a function $y(t)$ returning
@@ -113,17 +115,20 @@ this function using an ODE solver, such as #sundials.
 
 As of right now, our ball will fall until the end of time; we have not said
 anything about how it bounces when it hits the floor. To do so, we need a notion
-of discrete _events_. These are modelled by zero-crossings: we monitor a certain
+of _events_: we need to identify exactly when the ball hits the ground, so that
-value and stop when it goes from strictly negative to positive or null. For our
+we may take action to make it bounce. These events are modelled by
-purposes, we choose $-y$ as the monitored value, and call the zero-crossing
+zero-crossings: we monitor a certain value and stop when it goes from strictly
-event $z$. When $z$ occurs (i.e., when the ball touches the ground), we set the
+negative to positive or null. For our purposes, we choose $-y$ as the monitored
-speed $adot(y)$ to $-k dot "last"(adot(y))$, where $"last"(y)$ denotes the left
+value, and call the zero-crossing event $z$. When $z$ occurs (i.e., when the
-limit of $y$ (we cannot specify $adot(y)$ in terms of itself), and $k$ is a
+ball touches the ground), we set the speed $adot(y)$ to
-factor modelling the loss of inertia due to the collision (say, $0.8$). We can
+$-k dot #raw(lang: "zelus", "last")\(adot(y))$, where
-then resume the approximation of the solution.
+$#raw(lang: "zelus", "last")\(y)$ denotes the left limit of $y$ (we cannot
+specify $adot(y)$ in terms of itself), and $k$ is a factor modelling the loss of
+inertia due to the collision (say, $0.8$). We can then resume the approximation
+of the solution.
 
 @lst:ball.zls shows how such a model might be expressed in the concrete syntax
-of #zelus.
+of #zelus @cit:zelus_sync_lng_with_ode.
 
 #figure(placement: top, gap: 2em,
   ```zelus
@@ -135,22 +140,68 @@ of #zelus.
   caption: [The bouncing ball in #zelus]
 ) <lst:ball.zls>
 
-More formally, a hybrid system can be described as an automaton
+More formally, and as done in @cit:alg_ana_hyb_sys, a hybrid system $cal(H)$ can
+be defined as a graph whose nodes represent continuous behaviour, and whose
+edges represent discrete transitions:
+$ cal(H) = (L o c, V a r, E d g, A c t, I n v, I n i) $
+where:
+- $L o c$ is a finite set of _locations_;
+- $V a r$ is a finite set of _variables_;
+- $E d g subset.eq L o c times F times L o c$ is a finite set of _transitions_
+
 
 == Executing models
 
-Executing such a model is quite simple. There are two modes of execution:
+To execute such a model, we first compile it into a synchronous function, as
-discrete and continuous. In continuous mode, we call the ODE solver in order to
+described in @cit:sync_based_codegen_hyb_sys_lng. The details of this
-obtain an approximation of the variables defined through ODEs, and monitor for
+compilation step are not particularly relevant to our purposes, and can be
-zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in which
+ignored. What is more interesting is the output of this compilation step: a
-we perform computation steps as needed, until no other zero-crossing occurs, in
+single synchronous function. The simulation loop is then itself described as a
-which case we go back to the continuous mode, and repeat, as seen in @automaton.
+synchronous function operating on
 
-#figure(finite.automaton(
+#MENTION("Use of a single solver")
-  (D: (D: "cascade", C: "no cascade"),
+
-   C: (C: "no zero-crossing", D: "zero-crossing")),
+#pagebreak()
-  initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
+
-), caption: [High-level overview of the runtime], placement: top) <automaton>
+// The compilation of a hybrid model into a synchronous function is described in
+// detail in @cit:sync_based_codegen_hyb_sys_lng, but can be summarized quite
+// succintly as follows. By pairing this synchronous function with an
+// off-the-shelf ODE solver like #sundials, we can then simulate the dynamics of
+// the system. This is done by repeatedly performing execution steps according to
+// two different modes: discrete and continuous.
+
+// The continuous mode operates as follows: we first call the ODE solver in order
+// to approximate the dynamics of the model's continuous state.
+
+// Continuous steps first call the ODE solver to approximate the dynamics of the
+// model's continuous variables. The solver will return a function defined on a
+// time interval, which we then provide as input to the zero-crossing solver, which
+// will monitor the evolution of zero-crossing values along this interval. After
+// both solvers have been called, we choose what the next step's mode will be:
+// - if no zero-crossings have been detected, we output the entire solution
+//   provided by the ODE solver, and the next step remains continuous;
+// - if a zero-crossing occurs, we return the solution provided by the ODE solver
+//   up to the zero-crossing instant, and the next step becomes a discrete step.
+
+// Discrete steps perform state changes and side effects. We first call the model's
+// step function, which updates the state and outputs a value. We then decide what
+// the next step is. If a zero-crossing event occured due to the current step, the
+// next step is another discrete step. If no new event occured, we perform a
+// continuous step.
+
+// Executing such a model is quite simple. There are two modes of execution:
+// discrete and continuous. In continuous mode, we call the ODE solver in order
+// to obtain an approximation of the variables defined through ODEs, and monitor
+// for zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in
+// which we perform computation steps as needed, until no other zero-crossing
+// occurs, in which case we go back to the continuous mode, and repeat, as seen
+// in @automaton.
+
+// #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)
+// ), caption: [High-level overview of the runtime], placement: top) <automaton>
 
 = Runtime
 To solve this issue, we need to redefine what the runtime of our hybrid system
@@ -180,10 +231,10 @@ required by the assertion becomes a state variable.
 
 == Solvers as synchronous nodes
 == Simulations as synchronous nodes
-#TODO("talk about the new runtime")
+#MENTION("the new runtime")
 
 = Assertions
-#TODO("talk about how assertions are done")
+#MENTION("how assertions are done")
 
 #pagebreak()
 = Annex

doc/sources.bib

@@ -1,16 +1,137 @@
-@article{
+@article{cit:nonstd_sem_hyb_sys_mod,
-  ns_sem_benveniste_bourke_caillaud_pouzet,
+	title = {Non-standard semantics of hybrid systems modelers},
-  title={Non-Standard Semantics of Hybrid Systems Modelers},
+	journal = {Journal of Computer and System Sciences},
-  author={Benveniste, Albert and Bourke, Timothy and
+	volume = {78},
-          Caillaud, Benoıt and Pouzet, Marc},
+	number = {3},
-  year={2011},
+	pages = {877-910},
-  language={en}
+	year = {2012},
+	note = {In Commemoration of Amir Pnueli},
+	issn = {0022-0000},
+	doi = {https://doi.org/10.1016/j.jcss.2011.08.009},
+	url = {https://www.sciencedirect.com/science/article/pii/S0022000011001061},
+	author = {Albert Benveniste and Timothy Bourke and Benoît Caillaud and Marc
+	          Pouzet},
+	keywords = {Hybrid systems, Hybrid systems modelers, Non-standard analysis,
+	            Non-standard semantics, Constructive semantics, Kahn process
+	            networks, Compilation of hybrid systems},
+	abstract = {Hybrid system modelers have become a corner stone of complex
+	            embedded system development. Embedded systems include not only
+	            control components and software, but also physical devices. In
+	            this area, Simulink is a de facto standard design framework, and
+	            Modelica a new player. However, such tools raise several issues
+	            related to the lack of reproducibility of simulations
+	            (sensitivity to simulation parameters and to the choice of a
+	            simulation engine). In this paper we propose using techniques
+	            from non-standard analysis to define a semantic domain for hybrid
+	            systems. Non-standard analysis is an extension of classical
+	            analysis in which infinitesimal (the ε and η in the celebrated
+	            generic sentence ∀ε∃η… of college maths) can be manipulated as
+	            first class citizens. This approach allows us to define both a
+	            denotational semantics, a constructive semantics, and a Kahn
+	            Process Network semantics for hybrid systems, thus establishing
+	            simulation engines on a sound but flexible mathematical
+	            foundation. These semantics offer a clear distinction between the
+	            concerns of the numerical analyst (solving differential
+	            equations) and those of the computer scientist (generating
+	            execution schemes). We also discuss a number of practical and
+	            fundamental issues in hybrid system modelers that give rise to
+	            non-reproducibility of results, non-determinism, and undesirable
+	            side effects. Of particular importance are cascaded mode changes
+	            (also called “zero-crossings” in the context of hybrid systems
+	            modelers).},
 }
-@inbook{
+@inbook{cit:op_sem_hyb_sys,
-  opsem_lee_zheng,
+	address = {Berlin, Heidelberg},
+	series = {Lecture Notes in Computer Science},
 	title = {Operational Semantics of Hybrid Systems},
+	volume = {3414},
 	ISBN = {978-3-540-25108-8},
+	url = {http://link.springer.com/10.1007/978-3-540-31954-2_2},
+	DOI = {10.1007/978-3-540-31954-2_2},
+	abstractNote = {This paper discusses an interpretation of hybrid systems as
+	                executable models. A specification of a hybrid system for this
+	                purpose can be viewed as a program in a domain-specific
+	                programming language. We describe the semantics of HyVisual,
+	                which is such a domain-specific programming language. The
+	                semantic properties of such a language affect our ability to
+	                understand, execute, and analyze a model. We discuss several
+	                semantic issues that come in defining such a programming
+	                language, such as the interpretation of discontinuities in
+	                continuous-time signals, and the interpretation of
+	                discrete-event signals in hybrid systems, and the
+	                consequences of numerical ODE solver techniques. We describe
+	                the solution in HyVisual by giving its operational semantics.
+	                },
+	booktitle = {Hybrid Systems: Computation and Control},
+	publisher = {Springer Berlin Heidelberg},
 	author = {Lee, Edward A. and Zheng, Haiyang},
+	editor = {Morari, Manfred and Thiele, Lothar},
 	year = {2005},
-  language={en}
+	pages = {25–53},
+	collection = {Lecture Notes in Computer Science},
+	language = {en},
+}
+@inproceedings{cit:zelus_sync_lng_with_ode,
+	address = {Philadelphia Pennsylvania USA},
+	title = {Zélus: a synchronous language with ODEs},
+	ISBN = {978-1-4503-1567-8},
+	url = {https://dl.acm.org/doi/10.1145/2461328.2461348},
+	DOI = {10.1145/2461328.2461348},
+	abstractNote = { Z´elus is a new programming language for modeling systems
+	                that mix discrete logical time and continuous time behaviors.
+	                From a user’s perspective, its main originality is to extend
+	                an existing Lustre-like synchronous language with Ordinary
+	                Differential Equations (ODEs). The extension is conservative:
+	                any synchronous program expressed as dataflow equations and
+	                hierarchical automata can be composed arbitrarily with ODEs
+	                in the same source code. },
+	booktitle = { Proceedings of the 16th international conference on Hybrid
+	             systems: computation and control },
+	publisher = {ACM},
+	author = {Bourke, Timothy and Pouzet, Marc},
+	year = {2013},
+	month = apr,
+	pages = {113–118},
+	language = {en},
+}
+@inbook{cit:sync_based_codegen_hyb_sys_lng,
+	address = {Berlin, Heidelberg},
+	series = {Lecture Notes in Computer Science},
+	title = {A Synchronous-Based Code Generator for Explicit Hybrid Systems
+	         Languages},
+	volume = {9031},
+	rights = {http://www.springer.com/tdm},
+	ISBN = {978-3-662-46662-9},
+	url = {http://link.springer.com/10.1007/978-3-662-46663-6_4},
+	DOI = {10.1007/978-3-662-46663-6_4},
+	abstractNote = {Modeling languages for hybrid systems are cornerstones of
+	                embedded systems development in which software interacts with
+	                a physical environment. Sequential code generation from such
+	                languages is important for simulation efficiency and for
+	                producing code for embedded targets. Despite being routinely
+	                used in industrial compilers, code generation is rarely, if
+	                ever, described in full detail, much less formalized. Yet
+	                formalization is an essential step in building trustable
+	                compilers for critical embedded software development.},
+	booktitle = {Compiler Construction},
+	publisher = {Springer Berlin Heidelberg},
+	author = {Bourke, Timothy and Colaço, Jean-Louis and Pagano, Bruno and
+	          Pasteur, Cédric and Pouzet, Marc},
+	editor = {Franke, Björn},
+	year = {2015},
+	pages = {69–88},
+	collection = {Lecture Notes in Computer Science},
+	language = {en},
+}
+@article{cit:alg_ana_hyb_sys,
+	title = {The algorithmic analysis of hybrid systems},
+	author = {Alur, Rajeev and Courcoubetis, Costas and Halbwachs, Nicolas and
+	          Henzinger, Thomas A and Ho, P-H and Nicollin, Xavier and Olivero,
+	          Alfredo and Sifakis, Joseph and Yovine, Sergio},
+	journal = {Theoretical computer science},
+	volume = {138},
+	number = {1},
+	pages = {3--34},
+	year = {1995},
+	publisher = {Elsevier},
 }

exm/builtins/main.ml

@@ -73,21 +73,21 @@ let output =
 let sim =
   if !sundials then
     let open StatefulSundials in
-    let c = if !inplace then InPlace.csolve else Functional.csolve 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 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 !sample (run m s))
+    Hsim.Utils.run_until_n (output !sample (run m s))
   else
     let open StatefulRK45 in
-    let c = if !inplace then InPlace.csolve else Functional.csolve 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 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
-    run_until_n (output !sample n)
+    Hsim.Utils.run_until_n (output !sample n)
 
 let () = sim !stop !steps ignore
 

exm/zelus/ball/dune

@@ -4,9 +4,9 @@
    (:standard -w -a))))
 
 (rule
- (targets ballz.ml ballz.zci)
+ (targets ball.ml ball.zci)
  (deps
-  (:zl ballz.zls))
+  (:zl ball.zls))
  (action
   (run zeluc %{zl})))
 

exm/zelus/ball/main.ml

@@ -3,4 +3,4 @@ open Std
 
 let input _ = ()
 let output (now, (y, _, _)) = Format.printf "%.10e\t%.10e\n" now y
-let () = Runtime.go input Ballz.ball output
+let () = Runtime.go input Ball.ball output

exm/zelus/ball/ztypes.ml

@@ -0,0 +1,12 @@
+include Std
+include Ztypes
+include Solvers
+
+module type IGNORE = sig end
+module Defaultsolver : IGNORE = struct end
+
+module Zlsrun = struct
+  module Make (S : IGNORE) = struct
+    let go _ = ()
+  end
+end

exm/zelus/ballz/ztypes.ml

@@ -1,21 +0,0 @@
-include Std
-include Ztypes
-include Solvers
-
-module type IGNORE = sig end
-module Defaultsolver : IGNORE = struct end
-
-module Zlsrun = struct
-  module Make (S : IGNORE) = struct
-    let go s =
-      let s = Lift.lift_hsim s in
-      let open Hsim in
-      let state = (module State.InPlaceSimState : State.SimState) in
-      let solver =
-        Solver.solver (StatefulSundials.InPlace.csolve)
-        (Types.d_of_dc StatefulZ.InPlace.zsolve) in
-      let open Sim.Sim(val state) in
-      ()
-      (* run_until_n (Utils.ignore 0 (run s solver)) 30. 1 ignore *)
-  end
-end

exm/zelus/cradle/cradle.zls

@@ -0,0 +1,51 @@
+
+let mp6 = -. (3.1416 /. 6.)
+let g = 9.80665
+let l = 0.2
+
+let pi0 = mp6
+let pi1 = 0.
+let pi2 = 0.
+
+let acc x = -. g /. l *. (sin x)
+
+let hybrid cradle2() =
+  let rec der p0 = v0      init pi0 reset h01 -> last p1
+      and der v0 = acc(p0) init 0.0 reset h01 -> last v1
+      and der p1 = v1      init pi1 reset h01 -> last p0
+      and der v1 = acc(p1) init 0.0 reset h01 -> last v0
+      and h01 = up(last p0 -. last p1)
+      and init h = -0.1
+      and present h01 -> do h = -1.0 *. last h done
+             else do der h = 0.0 done
+  in (h, (p0, v0 /. 10.) , (p1, v1 /. 10.))
+
+let hybrid cradle3() =
+  let rec der p0 = v0      init pi0 reset h01 -> last p1
+      and der v0 = acc(p0) init 0.0 reset h01 -> last v1
+      and der p1 = v1      init pi1 reset h01 -> last p0 | h12 -> last p2
+      and der v1 = acc(p1) init 0.0 reset h01 -> last v0 | h12 -> last v2
+      and der p2 = v2      init pi2 reset h12 -> last p1
+      and der v2 = acc(p2) init 0.0 reset h12 -> last v1
+      and h01 = up(last p0 -. last p1)
+      and h12 = up(last p1 -. last p2)
+      and init h1 = -0.1
+      and present h01 -> do h1 = -1.0 *. last h1 done else do der h1 = 0.0 done
+      and init h2 = -0.1
+      and present h12 -> do h2 = -1.0 *. last h2 done else do der h2 = 0.0 done
+  in (p0, p1, p2, h1, h2)
+
+let node print(v, s) =
+  Format.printf "% .10e%s" v s
+
+let hybrid main() =
+  let der t = 1.0 init 0.0 in
+  let (p0, p1, p2, h1, h2) = cradle3() in
+  present (period(0.05)) -> (
+    print(t, "\t");
+	print(p0, "\t");
+	print(p1, "\t");
+	print(p2, "\t");
+	print(h1, "\t");
+	print(h2, "\n")
+  ); ()

exm/zelus/cradle/dune

@@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets cradle.ml cradle.zci format.zci)
+ (deps
+  (:zl cradle.zls)
+  (:zli format.zli))
+ (action
+  (run zeluc %{zli} %{zl})))
+
+(executable
+ (public_name cradle.exe)
+ (name main)
+ (libraries std))

exm/zelus/cradle/format.zli

@@ -0,0 +1,2 @@
+
+val printf : string -> float -> string -> unit

exm/zelus/cradle/main.ml

@@ -0,0 +1,30 @@
+
+open Std
+
+let input2 _ = ()
+let output2 (now, (h, (p0, v0), (p1, v1))) =
+  Format.printf "%.10e\t%.10e\t%.10e\n" now p0 p1
+
+let input3 _ = ()
+let output3 (now, (p0, p1, p2, h1, h2)) =
+  Format.printf "%.10e\t%.10e\t%.10e\t%.10e\t%.10e\t%.10e\n"
+  now p0 (p1 +. 1.0) (p2 +. 2.0) (h1 +. 3.0) (h2 +. 4.0)
+
+let input_main _ = ()
+let output_main (now, ()) = ()
+
+let three = ref false
+let main  = ref false
+
+let toggle y n () =
+  y := true;
+  List.iter (fun n -> n := false) n
+
+let () =
+  Runtime.register_args [
+    "-three", Arg.Unit (toggle three [main]), "\tUse the third model";
+    "-main", Arg.Unit (toggle main [three]), "\tUse the main model";
+  ];
+  if !main then Runtime.go input_main Cradle.main output_main
+  else if !three then Runtime.go input3 Cradle.cradle3 output3
+  else Runtime.go input2 Cradle.cradle2 output2

exm/zelus/cradle/ztypes.ml

@@ -0,0 +1,12 @@
+include Std
+include Ztypes
+include Solvers
+
+module type IGNORE = sig end
+module Defaultsolver : IGNORE = struct end
+
+module Zlsrun = struct
+  module Make (S : IGNORE) = struct
+    let go _ = ()
+  end
+end

exm/zelus/sincos/ztypes.ml

@@ -7,15 +7,6 @@ module Defaultsolver : IGNORE = struct end
 
 module Zlsrun = struct
   module Make (S : IGNORE) = struct
-    let go s =
+    let go _ = ()
-      let s = Lift.lift_hsim s in
-      let open Hsim in
-      let state = (module State.InPlaceSimState : State.SimState) in
-      let solver =
-        Solver.solver (StatefulSundials.InPlace.csolve)
-        (Types.d_of_dc StatefulZ.InPlace.zsolve) in
-      let open Sim.Sim(val state) in
-      ()
-      (* run_until_n (Utils.ignore 0 (run s solver)) 30. 1 ignore *)
   end
 end

exm/zelus/solve/dune

@@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets time.ml time.zci)
+ (deps
+  (:zl time.zls)
+  (:zli solve.zli))
+ (action
+  (run zeluc %{zli} %{zl})))
+
+(executable
+ (public_name time.exe)
+ (name main)
+ (libraries std))

exm/zelus/solve/main.ml

@@ -0,0 +1,10 @@
+
+open Std
+
+let input () = ()
+let output () = flush stdout
+
+let () =
+  Runtime.parse_args ();
+  Runtime.go_discrete input Time.main output
+

exm/zelus/solve/solve.zli

@@ -0,0 +1,23 @@
+
+type time = float
+type 'a value
+type 'a signal = 'a value option
+type 'a signal_t = ('a value * time) option
+
+val horizon : 'a value -> time
+val make    : time * (time -> 'a) -> 'a value
+val apply   : 'a value * time -> 'a
+
+val solve_ode45    : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t
+val solve_sundials : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t
+
+val synchr :
+  ('a signal -D-> 'b signal_t) -S->
+  ('a signal -D-> 'c signal_t) -S->
+  'a signal -D-> ('b * 'c) signal_t
+
+val iter    : int -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val iter_t  : int -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
+
+val check   : int -S-> ('a        -D-> bool) -S-> 'a signal_t -D-> unit
+val check_t : int -S-> (time * 'a -D-> bool) -S-> 'a signal_t -D-> unit

exm/zelus/solve/time.zls

@@ -0,0 +1,59 @@
+
+let epsilon = 0.0001
+
+let input _ = ()
+
+let hybrid sincos() =
+  let rec der sin = cos    init 0.0
+      and der cos = -. sin init 1.0
+  in (sin, cos)
+
+let sincos_ode45    = Solve.solve_ode45(sincos)
+let sincos_sundials = Solve.solve_sundials(sincos)
+let sincos_both     = Solve.synchr(sincos_ode45)(sincos_sundials)
+
+let hybrid ball () =
+  let rec der y  = y'    init 50.0 reset z -> 0.0
+      and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y')
+      and z = up(-. y)
+  in y
+
+let ball_ode45    = Solve.solve_ode45(ball)
+let ball_sundials = Solve.solve_sundials(ball)
+let ball_both     = Solve.synchr(ball_ode45)(ball_sundials)
+
+let node print_ball_both (now, (y1, y2)) =
+  print_float(now); print_string("\t");
+  print_float(y1);  print_string("\t");
+  print_float(y2);  print_string("\n");
+  ()
+
+let node print_sincos (now, (sin, cos)) =
+  print_float now; print_string "\t";
+  print_float sin; print_string "\t";
+  print_float cos; print_string "\n"
+
+let node print_sincos2 (now, ((sin1, cos1), (sin2, cos2))) =
+  print_float now;  print_string "\t";
+  print_float sin1; print_string "\t";
+  print_float sin2; print_string "\t";
+  print_float cos1; print_string "\t";
+  print_float cos2; print_string "\n"
+
+let node check_sincos (now, (sin, cos)) =
+  print_sincos (now, (sin, cos));
+  sin <= 1.0 +. epsilon && sin >= -1.0 -. epsilon &&
+  cos <= 1.0 +. epsilon && cos >= -1.0 -. epsilon
+
+let node check_sincos2 (now, ((sin1, cos1), (sin2, cos2))) =
+  print_sincos2 (now, ((sin1, cos1), (sin2, cos2)));
+  sin1 <= 1.0 +. epsilon && sin1 >= -1.0 -. epsilon &&
+  cos1 <= 1.0 +. epsilon && cos1 >= -1.0 -. epsilon &&
+  sin2 <= 1.0 +. epsilon && sin2 >= -1.0 -. epsilon &&
+  cos2 <= 1.0 +. epsilon && cos2 >= -1.0 -. epsilon
+
+let node main() =
+  let input = Some (Solve.make (30.0, input)) fby None in
+  let o = run sincos_sundials input in
+  Solve.check_t 100 check_sincos o
+

exm/zelus/solve/ztypes.ml

@@ -0,0 +1,16 @@
+include Std
+include Ztypes
+include Solvers
+
+module type IGNORE = sig end
+module Defaultsolver : IGNORE = struct end
+
+module Zlsrun = struct
+  module Make (S : IGNORE) = struct
+    let go _ = ()
+  end
+end
+
+module Stdlib = struct
+  type nonrec 'a option = 'a option
+end

exm/zelus_2024/ball/ball.ml

@@ -0,0 +1,58 @@
+(* The Zelus compiler, version 2024-dev
+  (2025-06-4-15:49) *)
+open Ztypes
+
+type ('e, 'd, 'c, 'b, 'a) ball =
+  { mutable time: 'e; mutable major: 'd; mutable up: 'c;
+    mutable y': 'b; mutable y: 'a }
+
+let ball =
+  let machine cstate = 
+    let alloc _ =
+      cstate.cmax <- cstate.cmax + 1;
+      cstate.zmax <- cstate.zmax + 1;
+      { time = -1.;
+        major = false;
+        up = { zin = false; zout = 1. };
+        y' = -1.;
+        y = { pos = -1.; der = 0. };
+      } in
+    let step self _ =
+      let cindex = cstate.cindex in
+      let cpos = ref cindex in
+      let zindex = cstate.zindex in
+      let zpos = ref zindex in
+      cstate.cindex <- cstate.cindex + 1;
+      cstate.zindex <- cstate.zindex + 1;
+      self.major <- cstate.major;
+      self.time <- cstate.time;
+      if cstate.major then
+        for i = cindex to 0 do Zls.set cstate.dvec i 0. done
+      else begin
+        self.y.pos <- Zls.get cstate.cvec !cpos;
+        cpos := !cpos + 1
+      end;
+      let result =
+        self.up.zout <- -. self.y.pos;
+        if self.up.zin then self.y' <- -0.8 *. self.y';
+        self.y.der <- self.y';
+        self.y.pos, self.y', self.up.zin in
+      cpos := cindex;
+      if cstate.major then begin
+         Zls.set cstate.cvec !cpos self.y.pos;
+         cpos := !cpos + 1;
+         self.up.zin <- false
+      end else begin
+        self.up.zin <- Zls.get_zin cstate.zinvec !zpos;
+        zpos := !zpos + 1
+      end;
+      zpos := zindex;
+      Zls.set cstate.zoutvec !zpos self.up.zout;
+      zpos := !zpos + 1;
+      Zls.set cstate.dvec !cpos self.y.der;
+      cpos := !cpos + 1;
+      result in
+    let reset self =
+      self.y.pos <- 50.; self.y' <- 0. in
+    Node { alloc; step; reset } in
+  machine

exm/zelus_2024/ball/dune

@@ -0,0 +1,9 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(executable
+ (public_name newball.exe)
+ (name main)
+ (libraries std))

exm/zelus_2024/ball/main.ml

@@ -0,0 +1,7 @@
+
+open Std
+
+let input _ = ()
+let output (now, (y, _, _)) = Format.printf "%.10e\t%.10e\n" now y
+let () = Runtime.go_2024 input Ball.ball output
+

exm/zelus_2024/ball/ztypes.ml

@@ -0,0 +1,12 @@
+include Std
+include Ztypes
+include Solvers
+
+module type IGNORE = sig end
+module Defaultsolver : IGNORE = struct end
+
+module Zlsrun = struct
+  module Make (S : IGNORE) = struct
+    let go _ = ()
+  end
+end

src/lib/hsim/sim.ml

@@ -171,34 +171,4 @@ module Sim (S : SimState) =
         update ms ss (set_idle st) in
 
       DNode { state; step; reset }
-
-    (** Run the model on the given input until the end of the input or until the
-        model stops answering. *)
-    let run_on (DNode n) input use =
-      let out = n.step n.state (Some input) in
-      let state = match out with None, s -> s | Some o, s -> use o; s in
-      let rec loop state =
-        let o, state = n.step state None in
-        match o with None -> () | Some o -> use o; loop state in
-      loop state
-
-    (** Run the model on multiple inputs. *)
-    let run_on_n (DNode n) inputs use =
-      ignore @@ List.fold_left (fun state i ->
-        let o, state = n.step state (Some i) in
-        begin match o with None -> () | Some o -> use o end;
-        let rec loop state =
-          let o, state = n.step state None in
-          match o with None -> state | Some o -> use o; loop state in
-        loop state) n.state inputs
-
-    (** Run the model autonomously until [h], or until the model stops
-        answering. *)
-    let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
-
-    (** Run the model autonomously until [h], split in [k] steps. *)
-    let run_until_n n h k =
-      let h = h /. float_of_int k in
-      run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () }))
-
   end

src/lib/hsim/types.ml

@@ -14,6 +14,10 @@ type 'a value =
   - [u: [0, h] -> α] *)
 type 'a signal = 'a value option
 
+(** A time signal with absolute timestamps added.
+    These represent the starting date for the value. *)
+type 'a signal_t = ('a value * time) option
+
 type ('s, 'p, 'a, 'b) drec =
   { state : 's;
     step  : 's -> 'a -> 'b * 's;

src/lib/hsim/utils.ml

@@ -7,6 +7,34 @@ let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v }
 let offset (u : time -> 'a) (now : time) : time -> 'a =
   fun t -> u (t +. now)
 
+(** Cut a value into two at a specified timestamp. *)
+let cutoff { h; u; c } t =
+  if t < 0.0 then
+    raise (Invalid_argument "Cutoff point cannot be negative");
+  if t > h then
+    raise (Invalid_argument "Cutoff point cannot be greater than horizon");
+  { h=t; c=Continuous; u }, { h=h -. t; c; u=fun n -> u (t +. n) }
+
+(** Join two values. *)
+let join { h=hl; u=ul; c=cl } { h=hr; u=ur; c=cr } =
+  let h = min hl hr in
+  let u = fun t -> ul t, ur t in
+  let c = match cl, cr with
+    | Continuous, Continuous -> Continuous
+    | _ -> Discontinuous in
+  { h; u; c }
+
+(** Map a function. *)
+let map_value f ({ u; _ } as v) =
+  { v with u=fun t -> f (u t) }
+
+(** Swap a pair. *)
+let swap v = map_value (fun (a, b) -> b, a) v
+
+let map_signal f v = Option.map (map_value f) v
+
+let swap_signal v = Option.map swap v
+
 (** Concatenate functions. *)
 let rec concat = function
   | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
@@ -67,7 +95,7 @@ let compose_sim
   DNode { state; step; reset }
 
 (** Track the evolution of a signal in time. *)
-let track : (unit, 'a signal, ('a value * time) option) dnode =
+let track : (unit, 'a signal, 'a signal_t) dnode =
   let state = 0.0 in
   let step now = function
     | None -> None, now
@@ -101,3 +129,31 @@ let do_and_reset (DNode m) (DNode n) f =
     m.reset ms mp, n.reset ns np in
   DNode { state; step; reset }
 
+(** Run a model on the given input until the end of the input or until the model
+    stops answering. *)
+let run_on (DNode n) input use =
+  let out = n.step n.state (Some input) in
+  let state = match out with None, s -> s | Some o, s -> use o; s in
+  let rec loop state =
+    let o, state = n.step state None in
+    match o with None -> () | Some o -> use o; loop state in
+  loop state
+
+(** Run the model on multiple inputs. *)
+let run_on_n (DNode n) inputs use =
+  Stdlib.ignore @@ List.fold_left (fun state i ->
+    let o, state = n.step state (Some i) in
+    begin match o with None -> () | Some o -> use o end;
+    let rec loop state =
+      let o, state = n.step state None in
+      match o with None -> state | Some o -> use o; loop state in
+    loop state) n.state inputs
+
+(** Run the model autonomously until [h], or until the model stops answering. *)
+let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
+
+(** Run the model autonomously until [h], split in [k] steps. *)
+let run_until_n n h k =
+  let h = h /. float_of_int k in
+  run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () }))
+

src/lib/solvers/statefulRK45.ml

@@ -7,7 +7,8 @@ module Functional =
   struct
     type ('state, 'vec) state = { state: 'state; vec: 'vec }
 
-    let csolve : (carray, carray) csolver_c =
+    let csolve () : (carray, carray) csolver_c =
+      Common.Debug.print "Instantiating RK45";
       let open Odexx.Ode45 in
 
       let state =
@@ -37,7 +38,8 @@ module InPlace =
   struct
     type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
 
-    let csolve : (carray, carray) csolver_c =
+    let csolve () : (carray, carray) csolver_c =
+      Common.Debug.print "Instantiating RK45";
       let open Odexx.Ode45 in
 
       let state =

src/lib/solvers/statefulSundials.ml

@@ -7,7 +7,8 @@ module Functional =
   struct
     type ('state, 'vec) state = { state : 'state; vec : 'vec }
 
-    let csolve : (carray, carray) csolver =
+    let csolve () : (carray, carray) csolver =
+      Format.printf "Instantiating Sundials";
       let open Cvode in
 
       let state =
@@ -37,7 +38,8 @@ module InPlace =
   struct
     type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
 
-    let csolve : (carray, carray) csolver =
+    let csolve () : (carray, carray) csolver =
+      Common.Debug.print "Instantiating Sundials";
       let open Cvode in
 
       let state =

src/lib/solvers/statefulZ.ml

@@ -7,7 +7,7 @@ module Functional =
   struct
     type ('state, 'vec) state = { state: 'state; vec: 'vec }
 
-    let zsolve : (carray, zarray, carray) zsolver_c =
+    let zsolve () : (carray, zarray, carray) zsolver_c =
       let open Illinois in
 
       let state =
@@ -38,7 +38,7 @@ module InPlace =
   struct
     type ('state, 'vec) state = { mutable state : 'state; mutable vec : 'vec }
 
-    let zsolve : (carray, zarray, carray) zsolver_c =
+    let zsolve () : (carray, zarray, carray) zsolver_c =
       let open Illinois in
 
       let state =

src/lib/std/lift.ml

@@ -6,8 +6,10 @@ open Ztypes
 type ('s, 'a) state =
   { mutable state : 's; mutable input : 'a option; mutable time : time }
 
-let lift f =
+let lift
-  let cstate =
+  (f : cstate -> (time * 'a, 'b) node)
+: (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode
+= let cstate =
     { cvec = cmake 0; dvec = cmake 0; cindex = 0; zindex = 0;
       cend = 0; zend = 0; cmax = 0; zmax = 0;
       zinvec = zmake 0; zoutvec = cmake 0;
@@ -61,7 +63,7 @@ let lift f =
     let o = f_step state (st.time, input) in
     o, st in
 
-  let reset _ ({ state; _ } as st) = f_reset state; st in
+  let reset () ({ state; _ } as st) = f_reset state; st in
 
   (* horizon *)
   let horizon { time; _ } =
@@ -140,3 +142,107 @@ let lift_hsim n =
     derivative state cstates ignore_der no_roots_in no_roots_out no_time; cstates in
 
   HNode { state; fder; fzer; fout; step; reset; horizon; jump; cget; cset; zset; csize; zsize }
+
+let lift_2024
+  (f : Ztypes.cstate_new -> (time * 'a, 'b) node)
+: (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode
+= let cstate =
+    { cvec = cmake 0; dvec = cmake 0; cindex = 0; zindex = 0;
+      cend = 0; zend = 0; cmax = 0; zmax = 0;
+      zinvec = zmake 0; zoutvec = cmake 0;
+      major = false; horizon = max_float; time=0.0 } in
+  let Node { alloc=f_alloc; step=f_step; reset=f_reset } = f cstate in
+  let state = { state = f_alloc (); input = None; time = 0.0 } in
+  let csize, zsize = cstate.cmax, cstate.zmax in
+  let no_roots_in  = zmake zsize in
+  let no_roots_out = cmake zsize in
+  let ignore_der   = cmake csize in
+  let cstates      = cmake csize in
+  cstate.cvec <- cstates;
+  f_reset state.state;
+
+  let no_time = -1.0 in
+
+  (* the function that compute the derivatives *)
+  let fder { state; time; _ } offset input y =
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time;
+    ignore (f_step state (time +. offset, input));
+    cstate.dvec in
+
+  (* the function that compute the zero-crossings *)
+  let fzer { state; time; _ } offset input y =
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time;
+    ignore (f_step state (time +. offset, input));
+    cstate.zoutvec in
+
+  (* the function which compute the output during integration *)
+  let fout { state; time; _ } offset input y =
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time;
+    f_step state (time +. offset, input) in
+
+  (* the function which compute a discrete step *)
+  let step ({ state; time; _ } as st) offset input =
+    st.input <- Some input;
+    st.time <- time +. offset;
+    cstate.time <- time;
+    cstate.major <- true;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    let o = f_step state (st.time, input) in
+    o, st in
+
+  let reset () ({ state; _ } as st) = f_reset state; st in
+
+  (* horizon *)
+  let horizon { time; _ } =
+    cstate.horizon -. time in
+
+  let jump _ = true in
+
+  (* the function which sets the zinvector into the *)
+  (* internal zero-crossing variables *)
+  let zset ({ state; input; _ } as st) zinvec =
+    cstate.major <- false;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.zinvec <- zinvec;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    st in
+
+  let cset ({ state; input; _ } as st) _ =
+    cstate.major <- false;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    st in
+
+  let cget { state; input; _ } =
+    cstate.major <- false;
+    cstate.horizon <- infinity;
+    cstate.zinvec <- no_roots_in;
+    cstate.zoutvec <- no_roots_out;
+    cstate.dvec <- ignore_der;
+    cstate.cindex <- 0;
+    cstate.zindex <- 0;
+    ignore (f_step state (no_time, Option.get input));
+    cstate.cvec in
+
+  HNode
+    { state; fder; fzer; step; fout; reset;
+      horizon; cset; cget; zset; zsize; csize; jump }

src/lib/std/runtime.ml

@@ -3,29 +3,53 @@ open Hsim.Types
 
 let sample   = ref 0
 let stop     = ref 10.0
+let sundials = ref false
 
-let options = [
+let opts = ref [
   "-sample", Arg.Set_int sample, "\tSampling frequency (default=0)";
   "-stop", Arg.Set_float stop, "\tStop time (default=10.0)";
   "-debug", Arg.Set Common.Debug.debug, "\tShow debug information";
+  "-sundials", Arg.Set sundials, "\tUse sundials cvode";
 ]
 
-let arg s =
+let anon = ref (fun s -> Format.eprintf "Unexpected argument: %s\n" s; exit 1)
-  Format.eprintf "Unexpected argument: %s\n" s; exit 1
 
 let usage = ""
 
+let register_args l = opts := !opts @ l
+let register_anon f = anon := f
+let parse_args () = Arg.parse (Arg.align !opts) !anon usage
+
 let go
   (input : time -> 'a)
   (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node)
   (output : (time * 'b) -> unit)
-= Arg.parse options arg usage;
+: unit
+= parse_args ();
   let input = { h=(!stop); c=Discontinuous; u=input } in
   let output o = List.iter output @@ Hsim.Utils.sample_tracked o !sample in
-  let model = Lift.lift model in
+  let solver = Solve.(if !sundials then Sundials else RK45) in
-  let open Hsim in
+  Hsim.Utils.run_on (Solve.build_sim solver model) input output
-  let solver = Solver.solver_c Solvers.StatefulRK45.InPlace.csolve
+
-                               Solvers.StatefulZ.InPlace.zsolve in
+let go_discrete
-  let open Sim.Sim(State.InPlaceSimState) in
+  (input : unit -> 'a)
-  let sim = Hsim.Utils.(compose (run model (d_of_dc solver)) track) in
+  (Ztypes.Node { alloc; step; reset } : ('a, 'b) Ztypes.node)
-  run_on sim input output
+  (output : 'b -> unit)
+: unit
+= parse_args ();
+  let mem = alloc () in
+  reset mem;
+  while true do
+    input () |> step mem |> output
+  done; ()
+
+let go_2024
+  (input : time -> 'a)
+  (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node)
+  (output : (time * 'b) -> unit)
+: unit
+= parse_args ();
+  let input = { h=(!stop); c=Discontinuous; u=input } in
+  let output o = List.iter output @@ Hsim.Utils.sample_tracked o !sample in
+  let solver = Solve.(if !sundials then Sundials else RK45) in
+  Hsim.Utils.run_on (Solve.build_sim_2024 solver model) input output

src/lib/std/runtime.mli

@@ -0,0 +1,24 @@
+
+open Hsim.Types
+
+val register_args : (string * Arg.spec * string) list -> unit
+val register_anon : (string -> unit) -> unit
+val parse_args : unit -> unit
+
+val go :
+  (time -> 'a) ->
+  (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) ->
+  ((time * 'b) -> unit) ->
+  unit
+
+val go_2024 :
+  (time -> 'a) ->
+  (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) ->
+  ((time * 'b) -> unit) ->
+  unit
+
+val go_discrete :
+  (unit -> 'a) ->
+  ('a, 'b) Ztypes.node ->
+  ('b -> unit) ->
+  unit

src/lib/std/solve.ml

@@ -0,0 +1,228 @@
+
+open Hsim
+open Types
+
+type nonrec 'a value = 'a value
+type nonrec 'a signal = 'a signal
+type nonrec 'a signal_t = 'a signal_t
+
+type time = float
+
+type solver = RK45 | Sundials
+
+(** Get a value's horizon [h] (reminder: a value is defined on [[0,h]]). *)
+let horizon { h; _ } = h
+
+(** Create a value from a horizon and function. *)
+let make (h, u) = { h; u; c=Discontinuous }
+
+(** Apply a value at a time t. *)
+let apply ({ u; h; _ }, t) =
+  if t > h then raise (Invalid_argument (Format.sprintf
+    "Requested time t=%.10e is greater than the horizon h=%.10e" t h)); 
+  u t
+
+let build_sim
+  (solver : solver)
+  (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node)
+: (unit *
+   ((Ztypes.cvec, Ztypes.dvec) Solver.ivp *
+    (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode
+= let model = Lift.lift model in
+  let solver = Hsim.Solver.solver
+    (match solver with
+     | RK45 -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve ()
+     | Sundials -> Solvers.StatefulSundials.InPlace.csolve ())
+    (d_of_dc @@ Solvers.StatefulZ.InPlace.zsolve ()) in
+  let open Hsim.Sim.Sim(Hsim.State.InPlaceSimState) in
+  let DNode s = Hsim.Utils.(compose (run model solver) track) in
+  DNode { s with reset=fun p -> s.reset (p, ())}
+
+let build_sim_2024
+  (solver : solver)
+  (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node)
+: (unit *
+   ((Ztypes.cvec, Ztypes.dvec) Solver.ivp *
+    (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode
+= let model = Lift.lift_2024 model in
+  let solver = Hsim.Solver.solver
+    (match solver with
+     | RK45 -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve ()
+     | Sundials -> Solvers.StatefulSundials.InPlace.csolve ())
+    (d_of_dc @@ Solvers.StatefulZ.InPlace.zsolve ()) in
+  let open Hsim.Sim.Sim(Hsim.State.InPlaceSimState) in
+  let DNode s = Hsim.Utils.(compose (run model solver) track) in
+  DNode { s with reset=fun p -> s.reset (p, ())}
+
+(** Lift a hybrid node into a discrete node on streams of functions. *)
+let solve
+  (solver : solver)
+  (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node)
+: ('a signal, 'b signal_t) Ztypes.node
+= let DNode sim = build_sim solver model in
+  let alloc () = ref sim.state in
+  let step s a = let b, s' = sim.step !s a in s := s'; b in
+  let reset _ = () in
+  Ztypes.Node { alloc; step; reset }
+
+let solve_2024
+  (solver : solver)
+  (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node)
+: ('a signal, 'b signal_t) Ztypes.node
+= let DNode sim = build_sim_2024 solver model in
+  let alloc () = ref sim.state in
+  let step s a = let b, s' = sim.step !s a in s := s'; b in
+  let reset _ = () in
+  Ztypes.Node { alloc; step; reset }
+
+let solve_ode45         m = solve      RK45     m
+let solve_ode45_2024    m = solve_2024 RK45     m
+let solve_sundials      m = solve      Sundials m
+let solve_sundials_2024 m = solve_2024 Sundials m
+
+(** Utility function for [synchr].
+
+    During synchronization, step the simulation that is lagging behind ([m]) and
+    join it with the stored value for the other ([n]).
+    Takes as arguments:
+    - The step method for [m];
+    - The input;
+    - The last stop times for [m] and [n];
+    - The state of [m];
+    - The stored value for [n].
+
+    Returns:
+    - The common output value up to the common reached date;
+    - The new reached date of [m];
+    - The stored value for [m];
+    - The stored value for [n]. *)
+let synchr_neq
+  (m_step : 'ms -> 'a signal -> 'b signal_t)
+  (input : 'a signal)
+  (m_stop : time) (n_stop : time) (m_state : 'ms) (n_value : 'c value)
+: ('b * 'c) signal_t * time * 'b signal * 'c signal
+= match m_step m_state input with
+  | None -> None, m_stop, None, Some n_value
+  | Some (m_value, m_start) ->
+      let m_stop = m_start +. m_value.h in
+      let m_value, n_value, m_rest, n_rest =
+        (* Three possible scenarios: *)
+        if m_stop < n_stop then begin
+          (* [m] is still behind [n]: cut off [n_value] at [m_stop'] *)
+          let n_value, n_rest = Utils.cutoff n_value m_value.h in
+          m_value, n_value, None, Some n_rest
+        end else if n_stop < m_stop then begin
+          (* [m] overtakes [n]: cut off [m_value] at [n_stop] *)
+          let m_value, m_rest = Utils.cutoff m_value (n_stop -. m_start) in
+          m_value, n_value, Some m_rest, None
+        end else
+          (* [m] reaches [n] exactly: *)
+          m_value, n_value, None, None in
+      let mn_value = Utils.join m_value n_value in
+      Some (mn_value, m_start), m_stop, m_rest, n_rest
+
+(** Utility function for [synchr].
+
+    During synchronization, step both simulations at the same time.
+    Takes as arguments:
+    - The step functions for both simulations;
+    - The input;
+    - The states of both simulations;
+    - The last stop times of both simulations.
+
+    Returns:
+    - The common output value up to the common reached date;
+    - The new stop times for both simulations;
+    - The new stored values for both simulations. *)
+let synchr_eq
+  (m_step : 'ms -> 'a signal -> 'b signal_t)
+  (n_step : 'ns -> 'a signal -> 'c signal_t)
+  (input : 'a signal) (m_state : 'ms) (n_state : 'ns)
+  (m_stop : time) (n_stop : time)
+: ('b * 'c) signal_t * time * time * 'b signal * 'c signal
+= match m_step m_state input, n_step n_state input with
+  | Some (m_value, m_start), Some (n_value, n_start) ->
+      let m_stop, n_stop = m_start +. m_value.h, n_start +. n_value.h in
+      let m_value, n_value, m_rest, n_rest =
+        if m_stop < n_stop then
+          let n_value, n_rest = Utils.cutoff n_value m_value.h in
+          m_value, n_value, None, Some n_rest
+        else if m_stop > n_stop then
+          let m_value, m_rest = Utils.cutoff m_value n_value.h in
+          m_value, n_value, Some m_rest, None
+        else m_value, n_value, None, None in
+      let mn_value = Utils.join m_value n_value in
+      Some (mn_value, min m_stop n_stop), m_stop, n_stop, m_rest, n_rest
+  | None, None -> None, m_stop, n_stop, None, None
+  | _ -> assert false
+
+(** Synchronize two simulations as much as possible. *)
+let synchr
+  (m : ('a signal, 'b signal_t) Ztypes.node)
+  (n : ('a signal, 'c signal_t) Ztypes.node)
+: ('a signal, ('b * 'c) signal_t) Ztypes.node
+= let Ztypes.Node { alloc=m_alloc; step=m_step; reset=m_reset } = m in
+  let Ztypes.Node { alloc=n_alloc; step=n_step; reset=n_reset } = n in
+  let alloc () =
+    ref ((0.0, None, m_alloc ()), (0.0, None, n_alloc ())) in
+  let step state input =
+    let (m_stop, m_value, m_state), (n_stop, n_value, n_state) = !state in
+    let m_stop, m_rest, m_state, n_stop, n_rest, n_state, output =
+      if m_stop < n_stop then
+        let n_value = Option.get n_value in
+        let output, m_stop, m_rest, n_rest =
+          synchr_neq m_step input m_stop n_stop m_state n_value in
+        m_stop, m_rest, m_state, n_stop, n_rest, n_state, output
+      else if m_stop > n_stop then
+        let m_value = Option.get m_value in
+        let output, n_stop, n_rest, m_rest =
+          synchr_neq n_step input n_stop m_stop n_state m_value in
+        let output = Option.map (fun (u, t) -> Utils.swap u, t) output in
+        m_stop, m_rest, m_state, n_stop, n_rest, n_state, output
+      else
+        let output, m_stop, n_stop, m_rest, n_rest =
+          synchr_eq m_step n_step input m_state n_state m_stop n_stop in
+        m_stop, m_rest, m_state, n_stop, n_rest, n_state, output in
+    state := (m_stop, m_rest, m_state), (n_stop, n_rest, n_state);
+    output in
+  let reset ({ contents=((_, _, ms), (_, _, ns)) } as s) =
+    n_reset ns; m_reset ms; s := (0.0, None, ms), (0.0, None, ns) in
+  Ztypes.Node { alloc; step; reset }
+
+(** Sample a value [n] times and iterate [f] on the samples. *)
+let iter n f =
+  let Ztypes.Node { alloc; step; reset } = f in
+  let step s =
+    Option.iter @@ fun (v, _) ->
+      List.iter (fun (_, v) -> step s v) @@ Utils.sample v n in
+  Ztypes.Node { alloc; step; reset }
+
+(** Sample a value [n] times and iterate [f] on the dated samples. *)
+let iter_t n f =
+  let Ztypes.Node { alloc; step; reset } = f in
+  let step s =
+    Option.iter @@ fun (v, h) ->
+      List.iter (fun (t, v) -> step s (t +. h, v)) @@ Utils.sample v n in
+  Ztypes.Node { alloc; step; reset }
+
+(** Sample a value [n] times and assert [f] on the samples. *)
+let check
+  (n : int)
+  (Ztypes.Node { alloc; step; reset } : ('a, bool) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let step s (now, v) =
+    try assert (step s v)
+    with Assert_failure _ ->
+      (Format.eprintf "Assertion failed at time %.10e\n" now; exit 1) in
+  iter_t n (Ztypes.Node { alloc; reset; step })
+
+(** Sample a value [n] times and assert [f] on the dated samples. *)
+let check_t
+  (n : int)
+  (Ztypes.Node { alloc; step; reset } : (time * 'a, bool) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let step s (now, v) =
+    try assert (step s (now, v))
+    with Assert_failure _ ->
+      (Format.eprintf "Assertion failed at time %.10e\n" now; exit 1) in
+  iter_t n (Ztypes.Node { alloc; reset; step })

src/lib/std/solve.mli

@@ -0,0 +1,61 @@
+
+type time = float
+type 'a value = 'a Hsim.Types.value
+type 'a signal = 'a value option
+type 'a signal_t = ('a value * time) option
+
+type solver = RK45 | Sundials
+
+val horizon : 'a value -> time
+val make    : time * (time -> 'a) -> 'a value
+val apply   : 'a value * time -> 'a
+
+val build_sim :
+  solver ->
+  (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) ->
+  (unit *
+   ((Ztypes.cvec, Ztypes.dvec) Hsim.Solver.ivp *
+    (Ztypes.cvec, Ztypes.zoutvec) Hsim.Solver.zc),
+   'a signal, 'b signal_t) Hsim.Types.dnode
+
+val build_sim_2024 :
+  solver ->
+  (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) ->
+  (unit *
+   ((Ztypes.cvec, Ztypes.dvec) Hsim.Solver.ivp *
+    (Ztypes.cvec, Ztypes.zoutvec) Hsim.Solver.zc),
+   'a signal, 'b signal_t) Hsim.Types.dnode
+
+val solve :
+  solver ->
+  (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+
+val solve_2024 :
+  solver ->
+  (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+
+val solve_ode45 :
+  (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+val solve_ode45_2024 :
+  (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+val solve_sundials :
+  (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+val solve_sundials_2024 :
+  (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> 
+  ('a signal, 'b signal_t) Ztypes.node
+
+val synchr :
+  ('a signal, 'b signal_t) Ztypes.node ->
+  ('a signal, 'c signal_t) Ztypes.node ->
+  ('a signal, ('b * 'c) signal_t) Ztypes.node
+
+val iter   : int -> ('a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node
+val iter_t : int -> (time * 'a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node
+
+val check   : int -> ('a, bool) Ztypes.node -> ('a signal_t, unit) Ztypes.node
+val check_t : int -> (time * 'a, bool) Ztypes.node -> ('a signal_t, unit) Ztypes.node

src/lib/std/solve.zli

@@ -0,0 +1,23 @@
+
+type time = float
+type 'a value
+type 'a signal = 'a value option
+type 'a signal_t = ('a value * time) option 
+
+val horizon : 'a value -> time
+val make    : time -> (time -> 'a) -> 'a value
+val apply   : 'a value -> time -> 'a
+
+val solve_ode45    : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t
+val solve_sundials : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t
+
+val synchr :
+  ('a signal -D-> 'b signal_t) -S->
+  ('a signal -D-> 'c signal_t) -S->
+  'a signal -D-> ('b * 'c) signal_t
+
+val iter    : int -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val iter_t  : int -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
+
+val check   : int -S-> ('a        -D-> bool) -S-> 'a signal_t -D-> unit
+val check_t : int -S-> (time * 'a -D-> bool) -S-> 'a signal_t -D-> unit

src/lib/std/ztypes.ml

@@ -67,6 +67,23 @@ type cstate =
     mutable major : bool; (* integration iff [major = false] *)
   }
 
+(* The interface with the ODE solver (new Zélus version). *)
+type cstate_new =
+  { mutable dvec    : dvec;          (* Derivative vector.                    *)
+    mutable cvec    : cvec;          (* Position vector.                      *)
+    mutable zinvec  : zinvec;        (* Zero-crossing result vector.          *)
+    mutable zoutvec : zoutvec;       (* Zero-crossing value vector.           *)
+    mutable cindex  : int;           (* Position in position vector.          *)
+    mutable zindex  : int;           (* Position in zero-crossing vector.     *)
+    mutable cend    : int;           (* End of position vector.               *)
+    mutable zend    : int;           (* End of zero-crossing vector.          *)
+    mutable cmax    : int;           (* Maximum size of position vector.      *)
+    mutable zmax    : int;           (* Maximum size of zero-crossing vector. *)
+    mutable horizon : float;         (* Next horizon.                         *)
+    mutable major   : bool;          (* Step mode: major <-> discrete.        *)
+    mutable time    : float;         (* Simulation time.                      *)
+  }
+
 (* A hybrid node is a node that is parameterised by a continuous state *)
 (* all instances points to this global parameter and read/write on it *)
 type ('a, 'b) hnode = cstate -> (time * 'a, 'b) node