feat: a LOT of stuff (final report, examples, simulation of a single assert, move from node instances to node definitions, etc.)

2025-08-20 18:20:46 +02:00 · 2025-08-20 18:20:46 +02:00 · f2c545ce2c
commit f2c545ce2c
parent ba5db5bd99
49 changed files with 12377 additions and 1898 deletions
--- a/doc/data/middle.csv
+++ b/doc/data/middle.csv
@ -0,0 +1,101 @@
+0.1,1.,1.
+0.2,1.09,0.9
+0.3,1.1606355,0.706355
+0.4,1.20740674526,0.467712452587
+0.5,1.23139725894,0.23990513678
+0.6,1.23688017859,0.0548291965681
+0.7,1.22854167208,-0.0833850651095
+0.8,1.21004121057,-0.185004615107
+0.9,1.18373429159,-0.263069189855
+1.,1.15086755111,-0.328667404789
+1.1,1.1118247,-0.390428511056
+1.2,1.06627366782,-0.455510321832
+1.3,1.01317876753,-0.530949002857
+1.4,0.950656415679,-0.625223518541
+1.5,0.875618537104,-0.750378785745
+1.6,0.783071598116,-0.925469389886
+1.7,0.664795417332,-1.18276180784
+1.8,0.506869445305,-1.57925972026
+1.9,0.285198697436,-2.21670747869
+2.,-0.0411442507765,-3.26342948213
+2.1,-0.529971005861,-4.88826755084
+2.2,-1.18926322222,-6.59292216359
+2.3,-1.70007492651,-5.10811704286
+2.4,-1.71110375487,-0.110288283653
+2.5,-1.6943904386,0.167133162707
+2.6,-1.67636818234,0.180222562559
+2.7,-1.65789428106,0.184739012827
+2.8,-1.63899329307,0.189009879909
+2.9,-1.61963873086,0.193545622084
+3.,-1.59979623066,0.19842500202
+3.1,-1.57942644945,0.203697812082
+3.2,-1.5584846181,0.20941831349
+3.3,-1.53691956293,0.215650551743
+3.4,-1.5146724274,0.222471355336
+3.5,-1.4916750512,0.229973761947
+3.6,-1.46784790356,0.2382714764
+3.7,-1.44309742078,0.247504827798
+3.8,-1.41731253591,0.257848848732
+3.9,-1.39036009703,0.269524388767
+4.,-1.36207873371,0.282813633255
+4.1,-1.3322705209,0.298082128026
+4.2,-1.30068946176,0.315810591386
+4.3,-1.26702528703,0.336641747347
+4.4,-1.23088021532,0.361450717085
+4.5,-1.1917348921,0.391453232213
+4.6,-1.14889727973,0.428376123701
+4.7,-1.10142396201,0.474733177166
+4.8,-1.04799551139,0.534284506253
+4.9,-0.986712969524,0.612825418644
+5.,-0.914754444531,0.719585249924
+5.1,-0.827775684935,0.869787595966
+5.2,-0.718829259023,1.08946425911
+5.3,-0.576368481233,1.42460777791
+5.4,-0.380576409527,1.95792071706
+5.5,-0.0972616434138,2.83314766113
+5.6,0.327343067813,4.24604711227
+5.7,0.938227766901,6.10884699088
+5.8,1.57630039075,6.38072623847
+5.9,1.72492948547,1.48629094717
+6.,1.70950968622,-0.1541979925
+6.1,1.6918164813,-0.176932049138
+6.2,1.67367963513,-0.18136846171
+6.3,1.6551400682,-0.185395669303
+6.4,1.63617378367,-0.189662845297
+6.5,1.61674960397,-0.194241796996
+6.6,1.59683206882,-0.199175351464
+6.7,1.57638103454,-0.204510342852
+6.8,1.55535084938,-0.210301851627
+6.9,1.53368929463,-0.216615547499
+7.,1.51133625499,-0.223530396355
+7.1,1.48822203852,-0.231142164754
+7.2,1.46426523263,-0.239568058841
+7.3,1.43936993788,-0.248952947563
+7.4,1.41342215567,-0.259477822069
+7.5,1.38628500991,-0.271371457617
+7.6,1.35779233631,-0.28492673593
+7.7,1.32773994949,-0.300523868279
+7.8,1.29587354655,-0.318664029389
+7.9,1.26187164466,-0.34001901883
+8.,1.22532103041,-0.365506142586
+8.1,1.18568065733,-0.396403730709
+8.2,1.14222727487,-0.434533824692
+8.3,1.09397137632,-0.482558985489
+8.4,1.03952350195,-0.544478743691
+8.5,0.976874889936,-0.626486120115
+8.6,0.903025524298,-0.738493656375
+8.7,0.813331641428,-0.896938828698
+8.8,0.700324127342,-1.13007514087
+8.9,0.551522102003,-1.48802025338
+9.,0.345434892204,-2.060872098
+9.1,0.0451454340179,-3.00289458186
+9.2,-0.405434196095,-4.50579630113
+9.3,-1.04021686075,-6.34782664655
+9.4,-1.63855504333,-5.9833818258
+9.5,-1.7164488535,-0.778938101736
+9.6,-1.70137983829,0.150690152195
+9.7,-1.68357260618,0.178072321076
+9.8,-1.66526259017,0.183100160094
+9.9,-1.64653268832,0.187299018484
+10.,-1.62736154517,0.191711431492
+
--- a/doc/data/sincos_discrete.csv
+++ b/doc/data/sincos_discrete.csv
--- a/doc/data/xsmall.csv
+++ b/doc/data/xsmall.csv
--- a/doc/remarques.typ
+++ b/doc/remarques.typ
@ -1,3 +1,7 @@
+#set par(justify: true)
+
+= Version 1.
+
 La phrase "it is less computationally intensive... the compiler". est mal dite;
 tu peux dire que c'est la methode classique, par exemple, implementee pour
 simulink et aussi Zelus.
@ -224,5 +228,121 @@ phases d'integration (le temps ronronne) et des pas discrets (reactions
 instantanees).

 Super. Continue ! --Marc
+#pagebreak()
+
+= Version 2.
+
+*L.230 - 235*. Le discours peut donner l'impression de confondre simuation et
+  modèle  mathématique. Tu écris un modèle mathématique idéalisé, c.à.d. avec un
+  choc élastique où la vitesse change instantanément de direction. Pour cela, la
+  manière de simuler doit changer. La, tu peux garder le discours.
+
+"_Since time is logical in discrete nodes_". Tu veux dire plutôt que l'on
+choisit de décrire des réactions en temps zéro (ou instantanées), c'est bien
+ça ? "_Since time is logical in discrete nodes, nothing tells us when, in
+continuous time, we should perform discrete steps._" Je ne comprends pas cette
+phrase !
+
+*L. 242-245*. Si tu veux définir le zéro-crossing, ne faudrait-il pas plutôt le
+faire avec une définition d'analyse d'abord et décrire ensuite, si besoin
+l'algorithme (mais informellement) ? L'intuition, c'est qu'il y a un
+zéro-crossing de $z$ en $t_0$, entre $t_"left"$ et $t_"right"$ lorsque il existe
+une boule d'épaisseur non nulle autour de $t_0$ (d'épaisseur $epsilon$) telle
+que pour tous les points a gauche de t0 (notons
+$t_0^- = { t | t <= t_0 and |t_0 - t| <= epsilon }$), $z(t_0^-) <= 0$ et tous
+les points a droite, $z(t_0^+) > 0$. Tim t'as donné l'article qui décrit l'algo.
+Illinois. Je ne l'ai pas sur moi. Le signal $z: [t_"left", t_"right"]$ a un
+zéro-crossing en t0 lorsque il existe un $epsilon > 0$ tq pour tout
+$alpha < epsilon$: ... . Qu'en penses-tu ?
+
+*Listing 5*. Tu n'expliques pas `last y'`.


+*P.10*. Tu devrais d'abord expliquer le modele mathématique, ce dont tu as
+besoin, puis l'exemple; ou bien l'inverse, expliquer d'abord intuitivement
+l'exemple, les ingrédients dont tu as besoin et que tu introduits, plus le
+modèle mathématique. Et ensuite, les choses dont tu as besoin pour simuler.
+E.g., connaitre $f_"der"$, $f_"zero"$, l'état initial continu, l'état discret,
+etc.
+
+*L.302*.
+Je trouve l'écriture $h in [0, v.h]$ un peu troublante. Est-ce qu'une notation
+$v\#h$ ou autre ne serait pas mieux ? Ou alors, utiliser plutôt un champ
+"horizon" que h, dans la notation en point. Ou "right" ?
+
+"_Multiple methods exist (Zélus uses the Illinois method [Sny53])_". La méthode
+Illinois est la plus connue et utilisée. C'est mal formulé. Multiple methods
+exist; one of the oldest and most-used method is the Illinois method. It is
+used, for example, in the Sundials CVODE suite (donner la référence; regarde
+le manuel). Simulink also used this method by default. Zelus also implements
+this method.
+
+*L.305-309*.
+C'est un peu confus. Dis simplement ce que fait une méthode de zéro-crossing et
+des ingrédients dont elle a besoin, avant d'expliquer, plus tard, comment on
+s'en sert dans la simulation.
+Tu as besoin de $g: T i m e -> X -> Z_o$; de $x_0: X$; de $t_"left"$, de
+$t_"right"$ et de $d e n s e: T i m e -> X$, défini sur cet intervalle. La
+fonction de zéro-crossing indique qui, sur le vecteur $Z_o$, traverse zéro. Tu
+peux signaler les difficultés éventuelles (on peut rater un événement, c'est
+sensible a la largeur de l'intervalle de détection, au fait que le nombre de
+traversées peut être paire et on rate l'événement, etc.). Et dire que on ne fait
+rien là dessus (pas plus que Simulink, Modelica, et les autres d'ailleurs).
+
+*L. 322-324*. Je ne comprends pas la définition entre la ligne 322 et 324 et pas
+bien le paragraphe précédent.
+
+*L. 418*. Tu veux plutôt dire que tu voudrais pouvoir agréger les solutions
+denses successives ($d k y$). Les solveurs classiques sont impératifs,
+c'est-à-dire qu'ils ont un état interne et que chaque appel à "step" le modifie
+physiquement. Pour pouvoir agréger plusieurs solutions successives, il faut
+qu'il fournisse un moyen de le faire. (Rmq: ce n'est pas forcément une "copie
+d'état" dont on a besoin. Pour RK, je l'ai fait en fournissant un moyen d'avoir
+une copie de $d k y$. Ça suffit.)
+
+*L. 459*. Si tu parles d'assertions dans les programmes, cite/lis les articles
+classiques car c'est une construction très ancienne des langages de
+programmation. Je ne suis pas spécialiste mais j'en ai lu deux vieux (de
+mémoire, un de Hoare; un de Dijsktra). Je suis sûr qu'il doit y avoir un article
+de survey que tout le monde cite la dessus. Je regarderai en rentrant. En somme,
+dis aussi pourquoi c'est intéressant/utile de pouvoir écrire des assertions et
+les difficultés que cela pose dans le cas présent.
+
+"_An important property of assertions is that they are transparent: their
+presence does not affect the result of the computation._" An expected feature of
+run-time assertions is that they should not affect the rest of the computation
+(except stopping execution when they are not fulfilled). We call them
+"transparent", in the sense that, running the program with or without, if no
+error is raised, should produce the same result.
+
+*Figure 9*. Je ne comprends pas la figure 9. Tu veux dire que, en Lustre, cela
+correspond à avoir un seul noeud qui calcule les deux en parallèle ? On n'écrira
+jamais le code de la partie droite. Relis l'article sur les assertions de Lustre
+(de memoire, AMAST 93). En Lustre, les assertions ont un role qui consiste à
+contraindre l'environnement, avant tout. C'est utile quand on veut vérifier des
+propriétes. En fait, une assertion se decompose en deux parties, une hypothèse
+(qui parle des entrées, incontrolâbles), et une conclusion (assume/guarantee).
+Dans notre cas, comme on veut s'en servir d'abord comme on le fait dans un
+langage généraliste et, pour le moment, pas pour faire de la vérification, on ne
+distingue pas les deux. Le programme de gauche est equivalent à
+```
+let node f (x) =
+ let v = ... in
+ let assertion = (let p = integr(0.0, v) in p >= 0.0) in
+ (v, assertion)
+```
+c.à.d. que assertion est un flot comme les autres. Est-ce autre chose que cela
+que tu veux dire ?
+
+*L. 477*.
+Pas vraiment. Un noeud Lustre avec une assertion a vérifier dynamiquement,
+s'implémente en calculant un flot supplementaire et en vérifiant qu'il est vrai
+à chaque instant. On le calcule donc avec le reste du code, comme on le fait
+habituellement pour les assertions dans les langages classiques. Attention: je
+ne parle pas ici de vérification formelle. On traite de manière particulière les
+assertions quand on cherche a vérifier une propriété. En Lustre, on considère
+que les assertions sont vraies et on vérifie qu'elles ne dépendent pas des
+sorties; sinon, elle ne sont pas causales. C'est pour cela, qu'il faut
+distinguer les hypothèses (assume) des résultats attendus (guarantee).
+
+Il n'y a aucun exemple ?
--- a/doc/rep.typ
+++ b/doc/rep.typ
--- a/doc/sources.bib
+++ b/doc/sources.bib
@ -212,3 +212,91 @@
 	volume = {4},
 	year = {1953},
 }
+@article{cit:assertion_hist,
+	title = {A historical perspective on runtime assertion checking in software
+	         development},
+	volume = {31},
+	ISSN = {0163-5948},
+	DOI = {10.1145/1127878.1127900},
+	abstractNote = {This report presents initial results in the area of software
+	                testing and analysis produced as part of the Software
+	                Engineering Impact Project. The report describes the
+	                historical development of runtime assertion checking,
+	                including a description of the origins of and significant
+	                features associated with assertion checking mechanisms, and
+	                initial findings about current industrial use. A future
+	                report will provide a more comprehensive assessment of
+	                development practice, for which we invite readers of this
+	                report to contribute information.},
+	number = {3},
+	journal = {ACM SIGSOFT Software Engineering Notes},
+	author = {Clarke, Lori A. and Rosenblum, David S.},
+	year = {2006},
+	month = may,
+	pages = {25–37},
+	language = {en},
+}
+@article{cit:assertion_axiom,
+	title = {An axiomatic basis for computer programming},
+	volume = {12},
+	ISSN = {0001-0782},
+	DOI = {10.1145/363235.363259},
+	abstractNote = {In this paper an attempt is made to explore the logical
+	                foundations of computer programming by use of techniques
+	                which were first applied in the study of geometry and have
+	                later been extended to other branches of mathematics. This
+	                involves the elucidation of sets of axioms and rules of
+	                inference which can be used in proofs of the properties of
+	                computer programs. Examples are given of such axioms and
+	                rules, and a formal proof of a simple theorem is displayed.
+	                Finally, it is argued that important advantage, both
+	                theoretical and practical, may follow from a pursuance of
+	                these topics.},
+	number = {10},
+	journal = {Commun. ACM},
+	author = {Hoare, C. A. R.},
+	year = {1969},
+	month = oct,
+	pages = {576–580},
+}
+@article{cit:assertion_lustre,
+	title = {Programming and verifying real-time systems by means of the
+	         synchronous data-flow language LUSTRE},
+	volume = {18},
+	rights = {
+	          https://ieeexplore.ieee.org/Xplorehelp/downloads/license-information/IEEE.html
+	          },
+	ISSN = {00985589},
+	DOI = {10.1109/32.159839},
+	number = {9},
+	journal = {IEEE Transactions on Software Engineering},
+	author = {Halbwachs, N. and Lagnier, F. and Ratel, C.},
+	year = {1992},
+	month = sep,
+	pages = {785–793},
+}
+@inbook{cit:hyb_auto,
+	address = {Berlin, Heidelberg},
+	title = {The Theory of Hybrid Automata},
+	ISBN = {978-3-642-59615-5},
+	url = {https://doi.org/10.1007/978-3-642-59615-5_13},
+	DOI = {10.1007/978-3-642-59615-5_13},
+	abstractNote = {A hybrid automaton is a formal model for a mixed
+	                discrete-continuous System. W e classify hybrid automata
+	                acoording to what questions about their behavior can be
+	                answered algorithmically. The Classification reveals
+	                structure on mixed discrete-continuous State Spaces that was
+	                previously studied on purely discrete state Spaces only. In
+	                particular, various classes of hybrid automata induce
+	                finitary trace equivalence (or similarity, or bisimilarity)
+	                relations on an uncountable State space, thus permitting the
+	                application of various model-checking techniques that were
+	                originally developed for finitestate Systems.},
+	booktitle = {Verification of Digital and Hybrid Systems},
+	publisher = {Springer},
+	author = {Henzinger, Thomas A.},
+	editor = {Inan, M. Kemal and Kurshan, Robert P.},
+	year = {2000},
+	pages = {265–292},
+	language = {en},
+}
--- a/exm/builtins/ball_assert.ml
+++ b/exm/builtins/ball_assert.ml
@ -0,0 +1,83 @@
+
+open Hsim.Types
+open Solvers.Zls
+
+(* let hybrid bouncing () = (x, y) where
+     rec der y  = y'   init y0
+     and der y' = -. g init y'0 reset z -> -0.8 *. last y'
+     and z = up (-. y)
+     and assert (up (-. (y +. epsilon))) *)
+
+let of_array (a : time array) : carray =
+  Bigarray.(Array1.of_array Float64 c_layout) a
+
+type state =
+  { zin : zarray;
+    lx  : carray; (* [| y'; y |] *)
+    i   : bool; }
+
+let g = -9.81
+let y0 = 50.0
+let y'0 = 0.0
+let x0 = 0.0
+let x'0 = 1.0
+
+let ball ()
+: (state, unit, unit, carray, carray, carray, zarray, carray) hrec
+= let zsize = 1 in
+  let csize = 2 in
+  let yd     = cmake csize in
+  let zout   = cmake zsize in
+  let zfalse = zmake 1 in
+  let state = { zin=zfalse; lx=of_array [| y'0; y0 |]; i=true } in
+  let fder _ _ () y = yd.{0} <- g; yd.{1} <- y.{0}; yd in
+  let fzer _ _ () y = zout.{0} <- -. y.{1}; zout in
+  let fout _ _ _ y = y in
+  let step s _ () =
+    let lx =
+      if s.zin.{0} = 1l then of_array [| -. 0.8 *. s.lx.{0}; s.lx.{1} |]
+      else s.lx in
+    s.lx, { zin=zfalse; lx; i=false } in
+  let reset () _ = state in
+  let horizon _ = max_float in
+  let jump _ = true in
+  let cset s lx = { s with lx } in
+  let cget s = s.lx in
+  let zset s zin = { s with zin } in
+  { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset;
+    csize; zsize }
+
+type astate = { zin_a: zarray }
+
+let aball epsilon
+: (astate, unit, state, bool, carray, carray, zarray, carray) hrec
+= let zsize = 1 in
+  let csize = 0 in
+  let zin_a = zmake zsize in
+  let yd   = cmake csize in
+  let zout = cmake zsize in
+  let state = { zin_a } in
+  let fder _ _ _ _ = yd in
+  let fzer _ _ st _ = zout.{0} <- -. (st.lx.{1} +. epsilon); zout in
+  let fout _ _ st _ = st.lx.{1} +. epsilon >= 0.0 in
+  let step { zin_a } _ st =
+    zin_a.{0} <> 1l && st.lx.{1} +. epsilon >= 0.0, { zin_a } in
+  let reset _ _ = state in
+  let horizon _ = max_float in
+  let jump _ = true in
+  let cset s _ = s in
+  let cget s = yd in
+  let zset _ zin_a = { zin_a } in
+  { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; csize; zsize }
+
+let errmsg_invalid = "Invalid arguments to model (needed: [float])"
+let errmsg_few     = "Too few arguments to model (needed: [float])"
+let errmsg_many    = "Too many arguments to model (needed: [float])"
+let init = function
+  | [eps] ->
+      let eps = try float_of_string eps
+      with Failure _ -> raise (Invalid_argument errmsg_invalid) in
+      let a = HNodeA { body=aball eps; assertions=[] } in
+      HNodeA { body=ball (); assertions=[a] }
+  | [] -> raise (Invalid_argument errmsg_few)
+  | _  -> raise (Invalid_argument errmsg_many)
--- a/exm/builtins/main.ml
+++ b/exm/builtins/main.ml
@ -4,19 +4,21 @@ open Solvers
 open Common
 open Types

-let sample   = ref 1
-let stop     = ref 10.0
-let accel    = ref false
-let inplace  = ref false
-let sundials = ref false
-let speed    = ref false
-let steps    = ref 1
-let model    = ref None
-let minstep  = ref None
-let maxstep  = ref None
-let mintol   = ref None
-let maxtol   = ref None
-let no_print = ref false
+let sample    = ref 1
+let stop      = ref 10.0
+let accel     = ref false
+let inplace   = ref false
+let sundials  = ref false
+let speed     = ref false
+let steps     = ref 1
+let model     = ref None
+let minstep   = ref None
+let maxstep   = ref None
+let mintol    = ref None
+let maxtol    = ref None
+let no_print  = ref false
+let no_assert = ref false
+let c_assert  = ref false

 let gt0i v i  = v := if i <= 0   then 1   else i
 let gt0f v f  = v := if f <= 0.0 then 1.0 else f
@ -29,7 +31,7 @@ let set_model s =
  | Some _ -> modelargs := s :: !modelargs

 let opts = [
-  "-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)";
+  "-sample", Arg.Int (gt0i sample), "n \tSample count (default=1)";
  "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
  "-debug", Arg.Set Debug.debug, "\tPrint debug information";
  "-accelerate", Arg.Set accel, "\tConcatenate continuous functions";
@ -41,7 +43,9 @@ let opts = [
  "-maxstep", Arg.String (opt maxstep), "\tSet maximum solver step length";
  "-mintol",  Arg.String (opt mintol),  "\tSet minimum solver tolerance";
  "-maxtol",  Arg.String (opt maxtol),  "\tSet maximum solver tolerance";
-  "-no-print", Arg.Set no_print, "\tDo not print output values";
+  "-noprint", Arg.Set no_print, "\tDo not print output values";
+  "-noassert", Arg.Set no_assert, "\tDo not check assertions";
+  "-cassert", Arg.Set c_assert, "\tCheck assertions continuously";
 ]

 let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
@ -53,12 +57,13 @@ let args = List.rev !modelargs
 let m =
  try match !model with
  | None          -> Format.eprintf "Missing model\n"; exit 2
-  | Some "ball"   -> Ball.init args
-  | Some "vdp"    -> Vdp.init args
-  | Some "sincos" -> Sincos.init args
-  | Some "sqrt"   -> Sqrt.init args
-  | Some "sin1x"  -> Sin1x.init args
-  | Some "sin1xd" -> Sin1x_der.init args
+  | Some "ball"   -> a_of_h @@ Ball.init args
+  | Some "vdp"    -> a_of_h @@ Vdp.init args
+  | Some "sincos" -> a_of_h @@ Sincos.init args
+  | Some "sqrt"   -> a_of_h @@ Sqrt.init args
+  | Some "sin1x"  -> a_of_h @@ Sin1x.init args
+  | Some "sin1xd" -> a_of_h @@ Sin1x_der.init args
+  | Some "aball"  -> Ball_assert.init args
  | Some s        -> Format.eprintf "Unknown model: %s\n" s; exit 2
  with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2

@ -73,21 +78,30 @@ 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 s = Solver.solver c (d_of_dc z) in
+    let z = if !inplace then InPlace.zsolve else Functional.zsolve in
+    let s = Solver.solver c (fun () -> d_of_dc (z ())) in
    let open Sim.Sim(val st) in
-    Hsim.Utils.run_until_n (output !sample (run m s))
+    let sim = if !no_assert then run (fun () -> h_of_a m) s
+              else if !c_assert then run_assert_continuous (fun () -> m) s
+              else run_assert_sample !sample (fun () -> m) s in 
+    Hsim.Utils.run_until_n (output !sample sim ())
  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
-    Hsim.Utils.run_until_n (output !sample n)
+    let sim =
+      if !no_assert then
+        if !accel then accelerate (fun () -> h_of_a m) s
+        else run (fun () -> h_of_a m) (fun () -> d_of_dc (s ()))
+      else
+        if !c_assert then run_assert_continuous (fun () -> m) (fun () -> d_of_dc (s ()))
+        else run_assert_sample !sample (fun () -> m) (fun () -> d_of_dc (s ())) in
+    Hsim.Utils.run_until_n (output !sample sim ())

-let () = sim !stop !steps ignore
+let () = ignore @@ sim !stop !steps ignore

--- a/exm/zelus/ballcos/ball.zci
+++ b/exm/zelus/ballcos/ball.zci
--- a/exm/zelus/ballcos/ball.zls
+++ b/exm/zelus/ballcos/ball.zls
@ -0,0 +1,41 @@
+(* Ball rolling on a cosine curve.                          *)
+(* Illustrates the impact of an observer on the simulation. *)
+
+let g  = 9.81
+let mu = 0.5  (* Friction coefficient. *)
+
+let hybrid ball(v0) = (x, v) where
+  rec der x = v init 0.0
+  and der v = a *. (cos x) init v0
+  and     a = g *. (sin x) -. mu *. v /. (cos x)
+
+let hybrid vdp_c(mu) = (x, y) where
+  rec der x = y init 1.0
+  and der y = (mu *. (1.0 -. (x *. x)) *. y) -. x init 1.0
+
+let hybrid print(p)(t, x, v, x', y) = () where
+  present(period(p)) -> do
+    () = print_endline(String.concat ",\t\t" (List.map string_of_float [t;x;v;x';y]))
+  done
+
+(* Changing the period for [print] changes the result. *)
+let hybrid main () = () where
+  rec der t  = 1.0 init 0.0
+  and (x, v) = ball(2.953)
+  and (x', y) = vdp_c(0.5)
+  and ()     = print(0.5)(t, x, v, x', y)
+
+(*
+let input _ = 2.953
+
+let node print_discrete (now, (x, v)) =
+  print_endline (String.concat ",\t\t" (List.map string_of_float [now;x;v]))
+
+let ball_discrete = Solve.solve_sundials(ball)
+
+let node main_discrete () =
+  let input = Some (Solve.make(30.0, input)) fby None in
+  let o = run ball_discrete input in
+  Solve.period'_t 1.0 print_discrete o
+*)
+
--- a/exm/zelus/ballcos/dune
+++ b/exm/zelus/ballcos/dune
@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets ball.ml ball.zci)
+ (deps
+  (:zl ball.zls)
+  (:zli solve.zli))
+ (action
+  (run zeluc %{zli} %{zl})))
+
+(executable
+ (public_name ballcos.exe)
+ (name main)
+ (libraries std))
--- a/exm/zelus/ballcos/main.ml
+++ b/exm/zelus/ballcos/main.ml
@ -0,0 +1,4 @@
+
+open Std
+
+let () = Runtime.go_discrete ignore Ball.main_discrete ignore
--- a/exm/zelus/ballcos/solve.zli
+++ b/exm/zelus/ballcos/solve.zli
@ -0,0 +1,27 @@
+
+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 sustain : 'a -> 'a value
+
+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
+
+val period'   : float -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val period'_t : float -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
--- a/exm/zelus/ballcos/tmp/ball.ml
+++ b/exm/zelus/ballcos/tmp/ball.ml
@ -0,0 +1,172 @@
+(* The Zelus compiler, version 2.2-dev
+  (2025-08-14-22:1) *)
+open Ztypes
+let g = 9.81
+
+let mu = 0.5
+
+type ('d , 'c , 'b , 'a) _ball =
+  { mutable major_68 : 'd ;
+    mutable i_72 : 'c ; mutable x_71 : 'b ; mutable v_70 : 'a }
+
+let ball (cstate_97:Ztypes.cstate) = 
+  
+  let ball_alloc _ =
+    cstate_97.cmax <- (+) cstate_97.cmax  2;
+    { major_68 = false ;
+      i_72 = (false:bool) ;
+      x_71 = { pos = 42.; der = 0. } ; v_70 = { pos = 42.; der = 0. } } in
+  let ball_step self ((time_67:float) , (v0_66:float)) =
+    ((let (cindex_98:int) = cstate_97.cindex in
+      let cpos_100 = ref (cindex_98:int) in
+      cstate_97.cindex <- (+) cstate_97.cindex  2 ;
+      self.major_68 <- cstate_97.major ;
+      (if cstate_97.major then
+       for i_1 = cindex_98 to 1 do Zls.set cstate_97.dvec  i_1  0. done
+       else ((self.x_71.pos <- Zls.get cstate_97.cvec  !cpos_100 ;
+              cpos_100 := (+) !cpos_100  1) ;
+             (self.v_70.pos <- Zls.get cstate_97.cvec  !cpos_100 ;
+              cpos_100 := (+) !cpos_100  1))) ;
+      (let (result_102:(float  * float)) =
+           (if self.i_72 then self.v_70.pos <- v0_66) ;
+           self.i_72 <- false ;
+           (let (a_69:float) =
+                (-.) (( *. ) g  (sin self.x_71.pos)) 
+                     ((/.) (( *. ) mu  self.v_70.pos)  (cos self.x_71.pos)) in
+            self.v_70.der <- ( *. ) a_69  (cos self.x_71.pos) ;
+            self.x_71.der <- self.v_70.pos ; (self.x_71.pos , self.v_70.pos)) in
+       cpos_100 := cindex_98 ;
+       (if cstate_97.major then
+        (((Zls.set cstate_97.cvec  !cpos_100  self.x_71.pos ;
+           cpos_100 := (+) !cpos_100  1) ;
+          (Zls.set cstate_97.cvec  !cpos_100  self.v_70.pos ;
+           cpos_100 := (+) !cpos_100  1)))
+        else (((Zls.set cstate_97.dvec  !cpos_100  self.x_71.der ;
+                cpos_100 := (+) !cpos_100  1) ;
+               (Zls.set cstate_97.dvec  !cpos_100  self.v_70.der ;
+                cpos_100 := (+) !cpos_100  1)))) ; result_102)):float * float) in
+   let ball_reset self  =
+     ((self.i_72 <- true ; self.x_71.pos <- 0.):unit) in
+  Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
+type ('c , 'b , 'a) _vdp_c =
+  { mutable major_75 : 'c ; mutable y_77 : 'b ; mutable x_76 : 'a }
+
+let vdp_c (cstate_103:Ztypes.cstate) = 
+  
+  let vdp_c_alloc _ =
+    cstate_103.cmax <- (+) cstate_103.cmax  2;
+    { major_75 = false ;
+      y_77 = { pos = 42.; der = 0. } ; x_76 = { pos = 42.; der = 0. } } in
+  let vdp_c_step self ((time_74:float) , (mu_73:float)) =
+    ((let (cindex_104:int) = cstate_103.cindex in
+      let cpos_106 = ref (cindex_104:int) in
+      cstate_103.cindex <- (+) cstate_103.cindex  2 ;
+      self.major_75 <- cstate_103.major ;
+      (if cstate_103.major then
+       for i_1 = cindex_104 to 1 do Zls.set cstate_103.dvec  i_1  0. done
+       else ((self.y_77.pos <- Zls.get cstate_103.cvec  !cpos_106 ;
+              cpos_106 := (+) !cpos_106  1) ;
+             (self.x_76.pos <- Zls.get cstate_103.cvec  !cpos_106 ;
+              cpos_106 := (+) !cpos_106  1))) ;
+      (let (result_108:(float  * float)) =
+           self.y_77.der <- (-.) (( *. ) (( *. ) mu_73 
+                                                 ((-.) 1. 
+                                                       (( *. ) self.x_76.pos 
+                                                               self.x_76.pos)))
+                                          self.y_77.pos)  self.x_76.pos ;
+           self.x_76.der <- self.y_77.pos ; (self.x_76.pos , self.y_77.pos) in
+       cpos_106 := cindex_104 ;
+       (if cstate_103.major then
+        (((Zls.set cstate_103.cvec  !cpos_106  self.y_77.pos ;
+           cpos_106 := (+) !cpos_106  1) ;
+          (Zls.set cstate_103.cvec  !cpos_106  self.x_76.pos ;
+           cpos_106 := (+) !cpos_106  1)))
+        else (((Zls.set cstate_103.dvec  !cpos_106  self.y_77.der ;
+                cpos_106 := (+) !cpos_106  1) ;
+               (Zls.set cstate_103.dvec  !cpos_106  self.x_76.der ;
+                cpos_106 := (+) !cpos_106  1)))) ; result_108)):float * float) in
+  
+  let vdp_c_reset self  =
+    ((self.y_77.pos <- 1. ; self.x_76.pos <- 1.):unit) in
+  Node { alloc = vdp_c_alloc; step = vdp_c_step ; reset = vdp_c_reset }
+type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _main =
+  { mutable i_96 : 'g ;
+    mutable i_95 : 'f ;
+    mutable major_79 : 'e ;
+    mutable h_94 : 'd ;
+    mutable i_92 : 'c ; mutable h_90 : 'b ; mutable t_80 : 'a }
+
+let main (cstate_109:Ztypes.cstate) = 
+  let Node { alloc = i_96_alloc; step = i_96_step ; reset = i_96_reset } = ball 
+  cstate_109 in 
+  let Node { alloc = i_95_alloc; step = i_95_step ; reset = i_95_reset } = vdp_c 
+  cstate_109 in
+  let main_alloc _ =
+    cstate_109.cmax <- (+) cstate_109.cmax  1;
+    { major_79 = false ;
+      h_94 = 42. ;
+      i_92 = (false:bool) ;
+      h_90 = (42.:float) ; t_80 = { pos = 42.; der = 0. };
+      i_96 = i_96_alloc () (* continuous *)  ;
+      i_95 = i_95_alloc () (* continuous *)  } in
+  let main_step self ((time_78:float) , ()) =
+    ((let (cindex_110:int) = cstate_109.cindex in
+      let cpos_112 = ref (cindex_110:int) in
+      cstate_109.cindex <- (+) cstate_109.cindex  1 ;
+      self.major_79 <- cstate_109.major ;
+      (if cstate_109.major then
+       for i_1 = cindex_110 to 0 do Zls.set cstate_109.dvec  i_1  0. done
+       else ((self.t_80.pos <- Zls.get cstate_109.cvec  !cpos_112 ;
+              cpos_112 := (+) !cpos_112  1))) ;
+      (let (result_114:unit) =
+           let h_93 = ref (infinity:float) in
+           (if self.i_92 then self.h_90 <- (+.) time_78  0.) ;
+           (let (z_91:bool) = (&&) self.major_79  ((>=) time_78  self.h_90) in
+            self.h_90 <- (if z_91 then (+.) self.h_90  0.01 else self.h_90) ;
+            h_93 := min !h_93  self.h_90 ;
+            self.h_94 <- !h_93 ;
+            self.i_92 <- false ;
+            (let ((x_82:float) , (v_81:float)) =
+                 i_96_step self.i_96 (time_78 , 2.953) in
+             let ((x'_83:float) , (y_84:float)) =
+                 i_95_step self.i_95 (time_78 , 0.0000000000000001) in
+             (begin match z_91 with
+                    | true ->
+                        let () =
+                            print_endline (String.concat ",\t\t" 
+                                                         (List.map string_of_float
+                                                                   
+                                                                   ((::) 
+                                                                    (
+                                                                    self.t_80.pos
+                                                                    ,
+                                                                    (
+                                                                    (::) 
+                                                                    (
+                                                                    x_82 ,
+                                                                    (
+                                                                    (::) 
+                                                                    (
+                                                                    v_81 ,
+                                                                    (
+                                                                    (::) 
+                                                                    (
+                                                                    x'_83 ,
+                                                                    (
+                                                                    (::) 
+                                                                    (
+                                                                    y_84 ,
+                                                                    ([]))))))))))))) in
+                        () | _ -> ()  end) ; self.t_80.der <- 1. ; ())) in
+       cstate_109.horizon <- min cstate_109.horizon  self.h_94 ;
+       cpos_112 := cindex_110 ;
+       (if cstate_109.major then
+        (((Zls.set cstate_109.cvec  !cpos_112  self.t_80.pos ;
+           cpos_112 := (+) !cpos_112  1)))
+        else (((Zls.set cstate_109.dvec  !cpos_112  self.t_80.der ;
+                cpos_112 := (+) !cpos_112  1)))) ; result_114)):unit) in 
+  let main_reset self  =
+    ((self.i_92 <- true ;
+      self.t_80.pos <- 0. ; i_96_reset self.i_96  ; i_95_reset self.i_95 ):
+    unit) in
+  Node { alloc = main_alloc; step = main_step ; reset = main_reset }
--- a/exm/zelus/ballcos/tmp/dune
+++ b/exm/zelus/ballcos/tmp/dune
@ -0,0 +1,8 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(executable
+ (name main_b)
+ (libraries zelus))
--- a/exm/zelus/ballcos/tmp/main_b.ml
+++ b/exm/zelus/ballcos/tmp/main_b.ml
@ -0,0 +1,31 @@
+open Ztypes
+open Zls
+
+(* simulation (continuous) function *)
+let main = 
+  let cstate = 
+    { dvec = cmake 0; cvec = cmake 0; zinvec = zmake 0; zoutvec = cmake 0; 
+      cindex = 0; zindex = 0; cend = 0; zend = 0; cmax = 0; zmax = 0; 
+      major = false; horizon = 0.0 } in 
+  let Node { alloc = alloc; step = hstep; reset = reset } = Ball.main cstate in 
+  let step mem cvec dvec zin t = 
+    cstate.major <- true; cstate.cvec <- cvec; cstate.dvec <- dvec; 
+    cstate.cindex <- 0; cstate.zindex <- 0; cstate.horizon <- infinity;  
+    hstep mem (t, ()) in
+  let derivative mem cvec dvec zin zout t = 
+    cstate.major <- false;  cstate.cvec <- cvec; cstate.dvec <- dvec; 
+    cstate.zinvec <- zin; cstate.zoutvec <- zout; cstate.cindex <- 0; 
+    cstate.zindex <- 0; ignore (hstep mem (t, ())) in
+  let crossings mem cvec zin zout t = 
+    cstate.major <- false;  cstate.cvec <- cvec; cstate.zinvec <- zin; 
+    cstate.zoutvec <- zout; cstate.cindex <- 0; cstate.zindex <- 0; 
+    ignore (hstep mem (t, ())) in
+  let maxsize mem = cstate.cmax, cstate.zmax in
+  let csize mem = cstate.cend in
+  let zsize mem = cstate.zend in
+  let horizon mem = cstate.horizon in
+  Hsim { alloc; step; reset; derivative; crossings; maxsize;  csize;  zsize; 
+         horizon };;
+(* instantiate a numeric solver *)
+module Runtime = Zlsrun.Make (Defaultsolver)
+let _ = Runtime.go main
--- a/exm/zelus/ballcos/ztypes.ml
+++ b/exm/zelus/ballcos/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
--- a/exm/zelus/brusselator/brusselator.zls
+++ b/exm/zelus/brusselator/brusselator.zls
@ -0,0 +1,28 @@
+(* The Brusselator. *)
+let hybrid brusselator(a, b) = (x, y) where
+	rec der x = a +. x *. x *. y -. b *. x -. x init 1.0
+	and der y = b *. x -. x *. x *. y init 1.0
+
+let pi = 3.141592653589793
+
+(* Add another oscillator. *)
+let hybrid harmonic(p) = x where
+	rec der x = v init 1.0
+	and der v = -2.0 *. pi *. x /. p init 0.0
+
+(* Putting the harmonic besides the brusselator changes the output of the first.
+   To visualize:
+
+    dune exec ./run.exe -- -speedup 1000 -maxstep 1.0 | feedgnuplot --stream --domain --lines
+*)
+let hybrid print(t, x) =
+	present (period (100.0)) ->
+		(print_endline (String.concat " " (List.map string_of_float [t; x])))
+	else ()
+
+let hybrid simu() =
+	let der t = 1.0 init 0.0 in
+	let (x, y) = brusselator(1.0, 2.001) in
+	let z = harmonic(1e-5) in
+	print(t, x)
+
--- a/exm/zelus/brusselator/dune.bak
+++ b/exm/zelus/brusselator/dune.bak
@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets brusselator.ml)
+ (deps
+  (:zl brusselator.zls))
+ (action
+  (run zeluc %{zl})))
+
+(executable
+ (name main)
+ (public_name brusselator.exe)
+ (libraries std)
+ (promote (until-clean)))
--- a/exm/zelus/brusselator/main.ml
+++ b/exm/zelus/brusselator/main.ml
@ -0,0 +1,6 @@
+
+open Std
+let input _ = ()
+let output (_, ()) = ()
+
+let () = Runtime.go input Brusselator.simu output
--- a/exm/zelus/odes/dune
+++ b/exm/zelus/odes/dune
@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets odes.ml odes.zci)
+ (deps
+  (:zl odes.zls)
+  (:zli solve.zli))
+ (action
+  (run zeluc %{zli} %{zl})))
+
+(executable
+ (public_name odes.exe)
+ (name main)
+ (libraries std))
--- a/exm/zelus/odes/main.ml
+++ b/exm/zelus/odes/main.ml
@ -0,0 +1,9 @@
+
+open Std
+
+let input _ = ()
+let output () = ()
+
+let () = Runtime.go_discrete input Odes.main output
+
+
--- a/exm/zelus/odes/odes.zls
+++ b/exm/zelus/odes/odes.zls
@ -0,0 +1,37 @@
+
+let hybrid vdp mu () = (x, y) where
+  rec der x = y init 1.0
+  and der y = (mu *. (1.0 -. (x *. x)) *. y) -. x init 1.0
+
+let hybrid sincos () = (sin, cos) where
+  rec der sin = cos    init 0.0
+  and der cos = -. sin init 1.0
+
+let hybrid both () = (x, y, s, c) where
+      (x, y) = vdp 5.0 ()
+  and (s, c) = sincos ()
+
+let vdp_d    = Solve.solve_sundials (vdp 5.0)
+let sincos_d = Solve.solve_sundials sincos
+let both_d   = Solve.solve_sundials both
+
+let main_d = Solve.synchr sincos_d both_d
+
+let node print_vdp (now, (x, y)) =
+  print_endline (String.concat ",\t\t" (List.map string_of_float [now;x;y]))
+
+let node print_sincos (now, (s, c)) =
+  print_endline (String.concat ",\t\t" (List.map string_of_float [now;s;c]))
+
+let node print_both (now, (x, y, s, c)) =
+  print_endline (String.concat ",\t\t" (List.map string_of_float [now;x;y;s;c]))
+
+let node print_main (now, ((s1, c1), (x, y, s2, c2))) =
+  print_endline (String.concat ",\t\t" (List.map string_of_float [now;x;y]))
+
+let input _ = ()
+
+let node main () =
+  let i = Some (Solve.make(1000.0, input)) fby None in
+  let o = run main_d i in
+  Solve.period'_t 0.01 print_main o
--- a/exm/zelus/odes/solve.zli
+++ b/exm/zelus/odes/solve.zli
@ -0,0 +1,26 @@
+
+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
+
+val period'   : float -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val period'_t : float -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
--- a/exm/zelus/odes/ztypes.ml
+++ b/exm/zelus/odes/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
--- a/exm/zelus/parallel/dune
+++ b/exm/zelus/parallel/dune
@ -0,0 +1,12 @@
+(env (dev (flags (:standard -w -a))))
+
+(rule
+ (targets main.ml parallel.ml parallel.zci)
+ (deps (:zl parallel.zls) (:zli solve.zli))
+ (action
+  (run zeluc -deps -s main %{zli} %{zl})))
+
+(executable
+ (public_name parallel.exe)
+ (name main)
+ (libraries std))
--- a/exm/zelus/parallel/parallel.zls
+++ b/exm/zelus/parallel/parallel.zls
@ -0,0 +1,34 @@
+(* Parallel simulation of harmonic oscillators.             *)
+(* Illustrates the impact of unrelated parallel simulation. *)
+
+let pi = 3.141592653589793
+
+let hybrid harmonic(p) = x where
+  rec der x = v init 1.0
+  and der v = -2.0 *. pi *. x /. p init 0.0
+
+let hybrid f () = (t, x, y) where
+  rec der t = 1.0 init 0.0
+  and x = harmonic(100.0)
+  and y = harmonic(1000.0)
+
+let hybrid main' () =
+  let t, x, y = f () in
+  present (period (0.001)) ->
+    print_endline (String.concat ",\t" (List.map string_of_float [t;x;y]))
+  else ()
+
+let hybrid f' () = harmonic(100.0)
+let hybrid g' () = (harmonic(100.0), harmonic(1e-3))
+
+let f_d = Solve.solve_sundials(f')
+let g_d = Solve.solve_sundials(g')
+let m   = Solve.synchr f_d g_d
+
+let node print (now, (xf, (xg, _))) =
+  print_endline (String.concat ",\t" (List.map string_of_float [now;xf;xg]))
+
+let input _ = ()
+let node main () =
+  let input = Some (Solve.make (100.0, input)) fby None in
+  Solve.period'_t 0.01 print (run m input)
--- a/exm/zelus/parallel/solve.zli
+++ b/exm/zelus/parallel/solve.zli
@ -0,0 +1,27 @@
+
+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 sustain : 'a -> 'a value
+
+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
+
+val period'   : float -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val period'_t : float -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
--- a/exm/zelus/parallel/ztypes.ml
+++ b/exm/zelus/parallel/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
--- a/exm/zelus/sincos/main.ml
+++ b/exm/zelus/sincos/main.ml
@ -1,6 +1,10 @@

 open Std

+(* let input _ = () *)
+(* let output (now, (sin, cos)) = Format.printf "%.10e\t%.10e\t%.10e\n" now sin cos *)
+(* let () = Runtime.go input Sincosz.g output *)
+
 let input _ = ()
-let output (now, (sin, cos)) = Format.printf "%.10e\t%.10e\t%.10e\n" now sin cos
-let () = Runtime.go input Sincosz.g output
+let output (now, sin, cos) = Format.printf "%.10e,%.10e,%.10e\n" now sin cos
+let () = Runtime.go_discrete input Sincosz.sincos output
--- a/exm/zelus/sincos/sincosz.zls
+++ b/exm/zelus/sincos/sincosz.zls
@ -1,4 +1,4 @@
-
+(*
 let hybrid g () = (sin, cos) where
  rec der sin = cos init 0.0
  and der cos = -. sin init 1.0
@ -15,3 +15,13 @@ let hybrid f () =
    print_float cos;
    print_newline ()
  ); ()
+*)
+let h = 0.01
+
+let node integr(x0, x') = (x) where
+	rec x = x0 -> pre(x +. x' *. h)
+
+let node sincos() = (now, sin, cos) where
+	rec sin = integr(0.0, cos)
+	and cos = integr(1.0, -. sin)
+	and now = integr(0.0, 1.0)
--- a/exm/zelus/solve/time.zls
+++ b/exm/zelus/solve/time.zls
@ -3,6 +3,16 @@ let epsilon = 0.0001

 let input _ = ()

+let time t = t
+
+let hybrid fsin  t =  s     where der s = cos t init 0.0
+let hybrid fcos  t =     c  where der c = -. (sin t) init 1.0
+let hybrid fboth t = (s, c) where (s, c) = (fsin t, fcos t)
+
+let fsind  = Solve.solve_sundials(fsin)
+let fcosd  = Solve.solve_sundials(fcos)
+let fbothd = Solve.solve_sundials(fboth)
+
 let hybrid sincos() =
  let rec der sin = cos    init 0.0
      and der cos = -. sin init 1.0
@ -22,38 +32,36 @@ 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 print1 (now, v) =
+  print_float(now); print_string "\t";
+  print_float v; print_string "\n"

-let node print_sincos (now, (sin, cos)) =
+let node print2 (now, (l, r)) =
  print_float now; print_string "\t";
-  print_float sin; print_string "\t";
-  print_float cos; print_string "\n"
+  print_float l; print_string "\t";
+  print_float r; print_string "\n"

-let node print_sincos2 (now, ((sin1, cos1), (sin2, cos2))) =
+let node print22 (now, ((ll, rl), (lr, rr))) =
  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"
+  print_float ll; print_string "\t";
+  print_float lr; print_string "\t";
+  print_float rl; print_string "\t";
+  print_float rr; print_string "\n"

 let node check_sincos (now, (sin, cos)) =
-  print_sincos (now, (sin, cos));
+  print2 (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)));
+  print22 (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
+  let input = Some (Solve.make (100.0, time)) fby None in
+  let o = run fbothd input in
+  Solve.iter_t 100 print2 o

--- a/exm/zelus/vdp/dune
+++ b/exm/zelus/vdp/dune
@ -0,0 +1,17 @@
+(env
+ (dev
+  (flags
+   (:standard -w -a))))
+
+(rule
+ (targets main.ml vdp.ml vdp.zci)
+ (deps
+  (:zl vdp.zls)
+  (:zli solve.zli))
+ (action
+  (run zeluc -deps -s main_d -o main %{zli} %{zl})))
+
+(executable
+ (public_name vdp.exe)
+ (name main)
+ (libraries std))
--- a/exm/zelus/vdp/solve.zli
+++ b/exm/zelus/vdp/solve.zli
@ -0,0 +1,27 @@
+
+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 sustain : 'a -> 'a value
+
+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
+
+val period'   : float -S-> ('a        -D-> unit) -S-> 'a signal_t -D-> unit
+val period'_t : float -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit
--- a/exm/zelus/vdp/vdp.zls
+++ b/exm/zelus/vdp/vdp.zls
@ -0,0 +1,54 @@
+
+let mu = 5.0
+
+let hybrid vdp_c() = (x, y) where
+  rec der x = y init 1.0
+  and der y = (mu *. (1.0 -. (x *. x)) *. y) -. x init 1.0
+
+let node forward(h)(x0, x') = x where
+  rec x = x0 fby (x +. h *. x')
+
+let node backward(h)(x0, x') = x where
+  rec x = x0 -> pre x +. h *. x'
+
+let node vdp_d(h)() = (x, y) where
+  rec x = backward(h)(1.0, y)
+  and y = forward(h)(1.0, (mu *. (1.0 -. (x *. x)) *. y) -. x)
+
+let stop_time = 50.0
+
+let node print (t, (x, y)) =
+  print_endline (String.concat ",\t" (List.map string_of_float [t;x;y]))
+
+let node main_d() =
+  let rec t = 0.0 -> pre t +. 0.001 in
+  print(t, vdp_d(0.001)())
+
+let node main_dc() =
+  let rec (t0, (x0, y0)) = ((0.0 -> pre t0 +. 0.1), vdp_d(0.1)()) in
+  let rec (t1, (x1, y1)) = ((0.0 -> pre t1 +. 0.2), vdp_d(0.2)()) in
+  let rec (t2, (x2, y2)) = ((0.0 -> pre t2 +. 0.3), vdp_d(0.3)()) in
+  let rec (t3, (x3, y3)) = ((0.0 -> pre t3 +. 0.4), vdp_d(0.4)()) in
+  let rec (t4, (x4, y4)) = ((0.0 -> pre t4 +. 0.5), vdp_d(0.5)()) in
+  let rec (t5, (x5, y5)) = ((0.0 -> pre t5 +. 0.6), vdp_d(0.6)()) in
+  let rec (t6, (x6, y6)) = ((0.0 -> pre t6 +. 0.7), vdp_d(0.7)()) in
+  let rec (t7, (x7, y7)) = ((0.0 -> pre t7 +. 0.8), vdp_d(0.8)()) in
+  let rec (t8, (x8, y8)) = ((0.0 -> pre t8 +. 0.9), vdp_d(0.9)()) in
+  let rec (t9, (x9, y9)) = ((0.0 -> pre t9 +. 1.0), vdp_d(1.0)()) in
+  print_endline (String.concat "\t" [string_of_float t0; "x0"; string_of_float x0; "y0"; string_of_float y0]);
+  print_endline (String.concat "\t" [string_of_float t1; "x1"; string_of_float x1; "y1"; string_of_float y1]);
+  print_endline (String.concat "\t" [string_of_float t2; "x2"; string_of_float x2; "y2"; string_of_float y2]);
+  print_endline (String.concat "\t" [string_of_float t3; "x3"; string_of_float x3; "y3"; string_of_float y3]);
+  print_endline (String.concat "\t" [string_of_float t4; "x4"; string_of_float x4; "y4"; string_of_float y4]);
+  print_endline (String.concat "\t" [string_of_float t5; "x5"; string_of_float x5; "y5"; string_of_float y5]);
+  print_endline (String.concat "\t" [string_of_float t6; "x6"; string_of_float x6; "y6"; string_of_float y6]);
+  print_endline (String.concat "\t" [string_of_float t7; "x7"; string_of_float x7; "y7"; string_of_float y7]);
+  print_endline (String.concat "\t" [string_of_float t8; "x8"; string_of_float x8; "y8"; string_of_float y8]);
+  print_endline (String.concat "\t" [string_of_float t9; "x9"; string_of_float x9; "y9"; string_of_float y9])
+
+let input _ = ()
+let vdp_s = Solve.solve_sundials vdp_c
+
+let node main_c() =
+  let o = run vdp_s (Some (Solve.make(stop_time, input)) fby None) in
+  Solve.period'_t 1.0 print o
--- a/exm/zelus/vdp/ztypes.ml
+++ b/exm/zelus/vdp/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
--- a/exm/zelus_2024/ball/ball.ml
+++ b/exm/zelus_2024/ball/ball.ml
@ -1,58 +1,72 @@
 (* The Zelus compiler, version 2024-dev
  (2025-06-4-15:49) *)
 open Ztypes
+type ('c, 'b, 'a) machine_17 =
+ { mutable _up_16: 'c;
+   mutable y'_12: 'b;
+   mutable y_11: 'a }

-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
+let (ball) =
+      let ball_10 =
+          let machine_17 cstate_18 = 
+            
+            let machine_17_alloc _ =
+              cstate_18.cmax <- (+) cstate_18.cmax 2;
+              cstate_18.zmax <- (+) cstate_18.zmax 1;
+              { _up_16 = { zin = false; zout = 1. };
+                y'_12 = { pos = (-1.); der = 0. };
+                y_11 = { pos = (-1.); der = 0. } } in
+            let machine_17_step self _ =
+              ((let cindex_19 = cstate_18.cindex in
+                let cpos_21 = ref (cindex_19:int) in
+                let zindex_20 = cstate_18.zindex in
+                let zpos_22 = ref (zindex_20:int) in
+                cstate_18.cindex <- (+) cstate_18.cindex 2;
+                cstate_18.zindex <- (+) cstate_18.zindex 1;
+                (if cstate_18.major
+                 then
+                   for i_1 = cindex_19 to 1
+                   do Zls.set cstate_18.dvec i_1 0. done
+                 else
+                   ((self.y'_12.pos <- Zls.get cstate_18.cvec !cpos_21;
+                     cpos_21 := (+) !cpos_21 1);
+                    (self.y_11.pos <- Zls.get cstate_18.cvec !cpos_21;
+                     cpos_21 := (+) !cpos_21 1)));
+                (let result_23 =
+                     self._up_16.zout <- (~-.) self.y_11.pos;
+                     self.y'_12.der <- (-9.81);
+                     (let z_13 = self._up_16.zin in
+                      let lx_15 = self.y'_12.pos in
+                      (match z_13 with
+                         | true ->
+                             let v_14 = lx_15 in
+                             self.y'_12.pos <- ( *. ) (-0.8) v_14 | _ -> () );
+                      self.y_11.der <- self.y'_12.pos;
+                      (self.y_11.pos, self.y'_12.pos, z_13)) in
+                 cpos_21 := cindex_19;
+                 (if cstate_18.major
+                  then
+                    (((Zls.set cstate_18.cvec !cpos_21 self.y'_12.pos;
+                       cpos_21 := (+) !cpos_21 1);
+                      (Zls.set cstate_18.cvec !cpos_21 self.y_11.pos;
+                       cpos_21 := (+) !cpos_21 1));
+                     ((self._up_16.zin <- false)))
+                  else
+                    (((self._up_16.zin <- Zls.get_zin cstate_18.zinvec
+                                                      !zpos_22;
+                       zpos_22 := (+) !zpos_22 1));
+                     zpos_22 := zindex_20;
+                     ((Zls.set cstate_18.zoutvec !zpos_22 self._up_16.zout;
+                       zpos_22 := (+) !zpos_22 1));
+                     ((Zls.set cstate_18.dvec !cpos_21 self.y'_12.der;
+                       cpos_21 := (+) !cpos_21 1);
+                      (Zls.set cstate_18.dvec !cpos_21 self.y_11.der;
+                       cpos_21 := (+) !cpos_21 1)))); result_23)):(float *
+                                                                   float *
+                                                                   bool)) in
+            let machine_17_reset self  =
+              ((self.y_11.pos <- 50.; self.y'_12.pos <- 0.):unit) in
+            Node { alloc = machine_17_alloc; step = machine_17_step;
+                                             reset = machine_17_reset } in
+            machine_17 in
+      ball_10
--- a/src/lib/hsim/sim.ml
+++ b/src/lib/hsim/sim.ml
@ -8,7 +8,7 @@ module Sim (S : SimState) =
    include S

    (** Discrete step. *)
-    let step_discrete
+    let dstep
      (s : ('a, 'b, 'ms, 'ss, 'zin) state)
      (step : 'ms -> time -> 'a -> 'b * 'ms)
      (hor : 'ms -> time)
@ -53,7 +53,7 @@ module Sim (S : SimState) =
      o, (set_last (Some o) (set_zin None s))

    (** Continuous step. *)
-    let step_continuous
+    let cstep
      (s : ('a, 'b, 'ms, 'ss, 'zin) state)
      (step : 'ss -> time -> (time * (time -> 'y) * 'zin option) * 'ss)
      (cset : 'ms -> 'y -> 'ms)
@ -85,18 +85,20 @@ module Sim (S : SimState) =

    (** Simulation of a model with any solver. *)
    let run
-      (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
-      (DNode s : ('y, 'yder, 'zin, 'zout) solver)
-    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
-    = let state = get_init m.state s.state in
+      (m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode def)
+      (s : ('y, 'yder, 'zin, 'zout) solver def)
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim def
+    = fun () ->
+      let HNode m = m () in
+      let DNode s = s () in
+      let state = get_init m.state s.state in
      let dstep ?(reinit=false) st =
-        let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset
+        let o, s = dstep st m.step m.horizon m.fder m.fzer m.cget m.zset
          m.csize m.zsize m.jump s.reset reinit in
        Some o, s in
      let cstep st =
-        let o, s, _ = step_continuous st s.step m.cset m.fout m.horizon in
+        let o, s, _ = cstep st s.step m.cset m.fout m.horizon in
        Some o, s in
-
      let step st = function
        | Some i ->
            let mode, now, stop = Discrete, 0.0, i.h in
@ -107,76 +109,26 @@ module Sim (S : SimState) =
            | Discrete   -> dstep st
            | Continuous -> cstep st
            else None, st in
-
      let reset (pm, ps) st =
        let ms = m.reset pm (get_mstate st) in
        let ss = s.reset ps (get_sstate st) in
        update ms ss (set_idle st) in
-
-      DNode { state; step; reset }
-
-    let rec run_assert :
-      'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a ->
-        (unit -> ('y, 'yder, 'zin, 'zout) solver) ->
-        ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim =
-    fun (HNodeA { body=m; assertions }) get_s ->
-      let DNode s = get_s () in
-      let al = List.map (fun a -> run_assert a get_s) assertions in
-      let state = get_init m.state s.state, al in
-
-      let dstep ?(reinit=false) (st, al) =
-        let o, st =
-          step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset m.csize
-          m.zsize m.jump s.reset reinit in
-        let al = List.map (fun (DNode a) ->
-          let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
-          DNode { a with state }) al in
-        Some o, (st, al) in
-
-      let cstep (st, al) =
-        let ({ h; _ } as o), st, u =
-          step_continuous st s.step m.cset m.fout m.horizon in
-        let al = List.map (fun (DNode a) ->
-          (* Step assertions repeatedly until they reach the horizon. *)
-          let rec step s =
-            let o, s = a.step s None in
-            match o with None -> s | Some _ -> step s in
-          let state = step (snd @@ a.step a.state (Some u)) in
-          DNode { a with state }) al in
-        (* Reset the model's state to the reached horizon. *)
-        let st = set_mstate (u.u h) st in
-        Some o, (st, al) in
-
-      let step (st, al) = function
-        | Some input ->
-            let mode, now, stop = Discrete, 0.0, input.h in
-            dstep (set_running ~mode ~input ~now ~stop st, al)
-        | None ->
-            if is_running st then match get_mode st with
-            | Discrete   -> dstep (st, al)
-            | Continuous -> cstep (st, al)
-            else None, (st, al) in
-
-      let reset (pm, ps) (st, al) =
-        let ms = m.reset pm (get_mstate st) in
-        let ss = s.reset ps (get_sstate st) in
-        let al = List.map (fun (DNode a) ->
-          DNode { a with state = a.reset (pm, ps) a.state }) al in
-        update ms ss (set_idle st), al in
-
      DNode { state; step; reset }

    let accelerate
-      (HNode  m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
-      (DNodeC s : ('y, 'yder, 'zin, 'zout) solver_c)
-    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
-    = let state = get_init m.state s.state in
+      (m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode def)
+      (s : ('y, 'yder, 'zin, 'zout) solver_c def)
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim def
+    = fun () ->
+      let HNode m = m () in
+      let DNodeC s = s () in
+      let state = get_init m.state s.state in
      let step_discrete ?(reinit=false) st =
-        let o, st = step_discrete st m.step m.horizon m.fder m.fzer m.cget
+        let o, st = dstep st m.step m.horizon m.fder m.fzer m.cget
          m.zset m.csize m.zsize m.jump s.reset reinit in
        Some o, st in
      let step_continuous st =
-        let o, st, _ = step_continuous st s.step m.cset m.fout m.horizon in
+        let o, st, _ = cstep st s.step m.cset m.fout m.horizon in
        o, st in

      let rec step st = function
@ -201,4 +153,121 @@ module Sim (S : SimState) =
        update ms ss (set_idle st) in

      DNode { state; step; reset }
+
+    let rec run_assert_continuous :
+      'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a def ->
+        (('y, 'yder, 'zin, 'zout) solver def) ->
+        ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim def =
+    fun m mk_s () ->
+      let HNodeA { body=m; assertions } = m () in
+      let DNode s = mk_s () in
+      let al = List.map (fun a -> run_assert_continuous (fun () -> a) mk_s ())
+        assertions in
+      let state = get_init m.state s.state, al in
+
+      let dstep ?(reinit=false) (st, al) =
+        let o, st =
+          dstep st m.step m.horizon m.fder m.fzer m.cget m.zset m.csize
+          m.zsize m.jump s.reset reinit in
+        let al = List.map (fun (DNode a) ->
+          let _, state = a.step a.state (Some (Utils.dot (get_mstate st))) in
+          DNode { a with state }) al in
+        Some o, (st, al) in
+
+      let cstep (st, al) =
+        let ({ h; _ } as o), st, u = cstep st s.step m.cset m.fout m.horizon in
+        let al = List.map (fun a -> Utils.run_on a u ignore) al in
+        let st = set_mstate (u.u h) st in
+        Some o, (st, al) in
+
+      let step (st, al) = function
+        | Some input ->
+            let mode, now, stop = Discrete, 0.0, input.h in
+            dstep (set_running ~mode ~input ~now ~stop st, al)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete   -> dstep (st, al)
+            | Continuous -> cstep (st, al)
+            else None, (st, al) in
+      let reset (pm, ps) (st, al) =
+        let ms = m.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        let al = List.map (fun (DNode a) ->
+          DNode { a with state = a.reset (pm, ps) a.state }) al in
+        update ms ss (set_idle st), al in
+      DNode { state; step; reset }
+
+    let rec run_assert_sample :
+      'a 'b. int ->
+        (('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a def) ->
+        (('y, 'yder, 'zin, 'zout) solver def) ->
+        ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim def =
+    fun n m get_s () ->
+      let HNodeA { body=m; assertions } = m () in
+      let DNode s = get_s () in
+      let al = List.map (fun a -> run_assert_sample n (fun () -> a) get_s ()) assertions in
+      let state = get_init m.state s.state, al in
+      let check v = List.iter (fun (_, v) -> assert v) @@ Utils.sample v n in
+      let dstep ?(reinit=false) (st, al) =
+        let o, st =
+          dstep st m.step m.horizon m.fder m.fzer m.cget m.zset m.csize m.zsize
+            m.jump s.reset reinit in
+        let al = List.map (fun (DNode a) ->
+          let o, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
+          Option.iter check o; DNode { a with state }) al in
+        Some o, (st, al) in
+      let cstep (st, al) =
+        let ({ h; _ } as o), st, u = cstep st s.step m.cset m.fout m.horizon in
+        let al = List.map (fun a -> Utils.run_on a u check) al in
+        let st = set_mstate (u.u h) st in
+        Some o, (st, al) in
+      let step (st, al) = function  
+        | Some input ->
+            let mode, now, stop = Discrete, 0.0, input.h in
+            dstep (set_running ~mode ~input ~now ~stop st, al)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete   -> dstep (st, al)
+            | Continuous -> cstep (st, al)
+            else None, (st, al) in
+      let reset (pm, ps) (st, al) =
+        let ms = m.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        let al = List.map (fun (DNode a) ->
+          DNode { a with state=a.reset (pm, ps) a.state}) al in
+        update ms ss (set_idle st), al in
+      DNode { state; step; reset }
+
+    let run_single_assert
+      (m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_sa def)
+      (s : ('y, 'yder, 'zin, 'zout) solver def)
+    : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim def
+    = fun () ->
+      let (HNodeSA { model=m; check=(DNode c) }) = m () in
+      let (DNode s) = s () in
+      let state = get_init m.state s.state, (DNode c) in
+      let dstep ?(reinit=false) (st, DNode c) =
+        let o, st = dstep st m.step m.horizon m.fder m.fzer m.cget m.zset
+          m.csize m.zsize m.jump s.reset reinit in
+        let b, state = c.step c.state @@ Some (Utils.dot @@ get_mstate st) in
+        assert b; Some o, (st, DNode { c with state }) in
+      let cstep (st, (DNode c)) =
+        let o, st, u = cstep st s.step m.cset m.fout m.horizon in
+        let b, state = c.step c.state (Some u) in 
+        assert b; Some o, (st, DNode { c with state }) in
+      let step (st, c) = function
+        | Some input ->
+            let mode, now, stop = Discrete, 0.0, input.h in
+            dstep (set_running ~mode ~input ~now ~stop st, c)
+        | None ->
+            if is_running st then match get_mode st with
+            | Discrete   -> dstep (st, c)
+            | Continuous -> cstep (st, c)
+            else None, (st, c) in
+      let reset (pm, ps) (st, DNode c) =
+        let ms = m.reset pm (get_mstate st) in
+        let ss = s.reset ps (get_sstate st) in
+        let c  = DNode { c with state=c.reset pm c.state } in
+        update ms ss (set_idle st), c in
+      DNode { state; step; reset }
  end
--- a/src/lib/hsim/solver.ml
+++ b/src/lib/hsim/solver.ml
@ -3,15 +3,15 @@ open Types

 (** An Initial Value Problem. *)
 type ('y, 'yder) ivp =
-  { init : 'y;                         (** [y₀]: initial value of y.          *)
-    fder : time -> 'y -> 'yder;        (** [dy/dt]: derivative of y.          *)
-    stop : time;                       (** Stop time.                         *)
+  { init : 'y;                      (** [y₀]: initial value of y.             *)
+    fder : time -> 'y -> 'yder;     (** [dy/dt]: derivative of y on [0,stop]. *)
+    stop : time;                    (** Stop time.                            *)
    size : int }

 (** A zero-crossing expression. *)
 type ('y, 'zout) zc =
-  { init : 'y;                         (** Value to watch for zero-crossings. *)
-    fzer : time -> 'y -> 'zout;        (** Zero-crossing function.            *)
+  { init : 'y;                      (** Value to watch for zero-crossings.    *)
+    fzer : time -> 'y -> 'zout;     (** Zero-crossing function.               *)
    size : int }

 (** An ODE solver is a synchronous function with:
@ -55,9 +55,12 @@ type ('y, 'yder, 'zin, 'zout) solver_c =
   time * (time -> 'y) * 'zin option) dnode_c

 (** Build a full solver from an ODE solver and a zero-crossing solver. *)
-let solver (DNode csolver : ('y, 'yder) csolver)
-           (DNode zsolver : ('y, 'zin, 'zout) zsolver)
-    : ('y, 'yder, 'zin, 'zout) solver =
+let solver (csolver : ('y, 'yder) csolver def)
+           (zsolver : ('y, 'zin, 'zout) zsolver def)
+    : ('y, 'yder, 'zin, 'zout) solver def =
+fun () ->
+  let DNode csolver = csolver () in
+  let DNode zsolver = zsolver () in
  let state = csolver.state, zsolver.state in
  let step (cstate, zstate) h =
    let (h', f), cstate = csolver.step cstate h in
@ -70,9 +73,12 @@ let solver (DNode csolver : ('y, 'yder) csolver)
  DNode { state; step; reset }

 (** Build a full solver supporting state copies. *)
-let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
-             (DNodeC zsolver : ('y, 'zin, 'zout) zsolver_c)
-    : ('y, 'yder, 'zin, 'zout) solver_c =
+let solver_c (csolver : ('y, 'yder) csolver_c def)
+             (zsolver : ('y, 'zin, 'zout) zsolver_c def)
+    : ('y, 'yder, 'zin, 'zout) solver_c def =
+fun () ->
+  let DNodeC csolver = csolver () in
+  let DNodeC zsolver = zsolver () in
  let state = csolver.state, zsolver.state in
  let step (cstate, zstate) h =
    let (h', f), cstate = csolver.step cstate h in
--- a/src/lib/hsim/types.ml
+++ b/src/lib/hsim/types.ml
@ -1,4 +1,6 @@

+type 'a def = unit -> 'a
+
 type time = float
 type continuity = Continuous | Discontinuous

@ -37,6 +39,10 @@ type ('s, 'p, 'a, 'b) drec_c =
 type ('p, 'a, 'b) dnode_c =
  DNodeC : ('s, 'p, 'a, 'b) drec_c -> ('p, 'a, 'b) dnode_c

+(** The simulation of a hybrid system is a synchronous function on streams of
+    functions. *)
+type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
+
 type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
  { state: 's;
    step : 's -> time -> 'a -> 'b * 's;     (** Step function.                *)
@ -52,7 +58,7 @@ type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
    csize : int;
    zsize : int }

-(** A hybrid node. *)
+(** A hybrid node instance. *)
 type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
  HNode : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec ->
          ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
@ -61,12 +67,15 @@ type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
 type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
  HNodeA : {
    body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
-    assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list
+    assertions : ('p, 's, bool, 'y, 'yder, 'zin, 'zout) hnode_a list
  } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a

-(** The simulation of a hybrid system is a synchronous function on streams of
-    functions. *)
-type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
+(** A hybrid node and a simulation of its assertions. *)
+type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_sa =
+  HNodeSA : {
+    model : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
+    check : ('p, 's signal, bool) dnode;
+  } -> ('p, 'a, 'b, 'y,' yder, 'zin,' zout) hnode_sa

 (* Utils *)

--- a/src/lib/hsim/utils.ml
+++ b/src/lib/hsim/utils.ml
@ -58,8 +58,24 @@ let sample { h; u; _ } n =
 let sample_tracked (o, t) n =
  List.map (fun (h, v) -> h +. t, v) @@ sample o n

+let period ?(offset=0.0) { h; u; _ } p =
+  let rec step now =
+    if now >= h then [], now -. h
+    else let l, o = step (now +. p) in
+    (now, u now) :: l, o in
+  if h <= 0.0 then [(0.0, u 0.0)], offset
+  else step offset
+
+let period_tracked ?(offset=0.0) (v, t) p =
+  let l, o = period ~offset v p in
+  List.map (fun (h, v) -> h +. t, v) l, o
+
 (** Compose two nodes together. *)
-let compose (DNode m) (DNode n) =
+let compose (m : ('p, 'a, 'b) dnode def) (n : ('q, 'b, 'c) dnode def)
+  : ('p * 'q, 'a, 'c) dnode def =
+fun () ->
+  let DNode m = m () in
+  let DNode n = n () in
  let state = m.state, n.state in
  let step (ms, ns) i =
    let v, ms = m.step ms i in
@ -73,9 +89,13 @@ let compose (DNode m) (DNode n) =

 (** Compose two simulations. *)
 let compose_sim
-  (DNode m : ('p, 'a, 'b) sim)
-  (DNode n : ('q, 'b, 'c) sim)
-= let state = m.state, n.state in
+  (m : ('p, 'a, 'b) sim def)
+  (n : ('q, 'b, 'c) sim def)
+  : ('p * 'q, 'a, 'c) sim def =
+fun () ->
+  let DNode m = m () in
+  let DNode n = n () in
+  let state = m.state, n.state in
  let step (ms, ns) = function
    | Some i ->
        let v, ms = m.step ms (Some i) in
@ -95,7 +115,9 @@ let compose_sim
  DNode { state; step; reset }

 (** Track the evolution of a signal in time. *)
-let track : (unit, 'a signal, 'a signal_t) dnode =
+let track
+  : (unit, 'a signal, 'a signal_t) dnode def =
+fun () ->
  let state = 0.0 in
  let step now = function
    | None -> None, now
@ -104,17 +126,29 @@ let track : (unit, 'a signal, 'a signal_t) dnode =
  DNode { state; step; reset }

 (** Apply a function to a signal. *)
-let map f =
+let sigmap (f : 'a -> 'b)
+  : (unit, 'a option, 'b option) dnode def =
+fun () ->
  let state = () in
  let step () = function None -> None, () | Some v -> Some (f v), () in
  let reset () () = () in
  DNode { state; step; reset }

-let ignore _ n =
+let dimap (i : 'a -> 'b) (o : 'c -> 'd) (n : ('p, 'b, 'c) dnode def)
+  : ('p, 'a, 'd) dnode def
+= fun () ->
+  let DNode { state; step; reset } = n () in
+  DNode { state; reset; step=fun s a -> let c, s = step s (i a) in o c, s }
+
+(* let preprocess (f : 'a -> 'b) (n : ('p, 'b, 'c) dnode def) *)
+  (* : ('p, 'a, 'c) dnode def *)
+(* = *) 
+
+let ignore _ n () =
  let state = () in
  let step () = function None -> None, () | Some _ -> Some (), () in
  let reset () () = () in
-  let DNode n = compose n @@ DNode { state; step; reset } in
+  let DNode n = compose n (fun () -> DNode { state; step; reset }) () in
  DNode { n with reset=fun p -> n.reset (p, ()) }

 let do_and_reset (DNode m) (DNode n) f =
@ -136,7 +170,7 @@ let run_on (DNode n) input use =
  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
+    match o with None -> DNode { n with state } | Some o -> use o; loop state in
  loop state

 (** Run the model on multiple inputs. *)
@ -156,4 +190,3 @@ let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
 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 _ -> () }))
-
--- a/src/lib/solvers/statefulZ.ml
+++ b/src/lib/solvers/statefulZ.ml
@ -7,7 +7,8 @@ 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 def =
+      fun () ->
      let open Illinois in

      let state =
@ -38,7 +39,8 @@ 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 def =
+      fun () ->
      let open Illinois in

      let state =
--- a/src/lib/std/lift.ml
+++ b/src/lib/std/lift.ml
@ -6,243 +6,202 @@ open Ztypes
 type ('s, 'a) state =
  { mutable state : 's; mutable input : 'a option; mutable time : time }

-let lift
-  (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;
-      major = false; horizon = max_float } in
-  let Node { alloc=f_alloc; step=f_step; reset=f_reset } = f cstate in
+(* Wrappers around the [step] function. *)
+
+let mkfder (c : cstate) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time +. o, i)); c.dvec
+
+let mkfzer (c : cstate) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time +. o, i)); c.zoutvec
+
+let mkfout (c : cstate) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0; f state (time +. o, i)
+
+let mkstep (c : cstate) f der zin zout ({ state; time; _ } as st) o i =
+  st.input <- Some i; st.time <- time +. o; c.major <- true;
+  c.horizon <- infinity; c.zinvec <- zin; c.zoutvec <- zout; c.dvec <- der;
+  c.cindex <- 0; c.zindex <- 0; let o = f state (st.time, i) in o, st
+
+let mkzset (c : cstate) f der zout time ({ state; input; _ } as s) zin =
+  c.major <- false; c.zoutvec <- zout; c.dvec <- der; c.zinvec <- zin;
+  c.cindex <- 0; c.zindex <- 0; ignore (f state (time, Option.get input)); s
+
+let mkcset (c : cstate) f der zin zout time ({ state; input; _ } as st) _ =
+  c.major <- false; c.horizon <- infinity; c.zinvec <- zin; c.zoutvec <- zout;
+  c.dvec <- der; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time, Option.get input)); st
+
+let mkcget (c : cstate) f der zin zout time { state; input; _ } =
+  c.major <- false; c.horizon <- infinity; c.zinvec <- zin; c.zoutvec <- zout;
+  c.dvec <- der; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time, Option.get input)); c.cvec
+
+(** Lift the inner record of a node. Needed to keep the typechecker aware of the
+    existential type. *)
+let lift_inner (cs : cstate) (f : ('s, time * 'a, 'b) node_rec)
+  : (('s, 'a) state, unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) hrec
+= let { alloc=f_alloc; step=f_step; reset=f_reset } = f 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 csize, zsize = cs.cmax, cs.zmax in
+  let no_zin, no_zout, no_der = zmake zsize, cmake zsize, cmake csize in
+  cs.cvec <- cmake csize; 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;
-    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;
-    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;
-    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.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 fder = mkfder cs f_step no_der no_zin no_zout in
+  let fzer = mkfzer cs f_step no_der no_zin no_zout in
+  let fout = mkfout cs f_step no_der no_zin no_zout in
+  let step = mkstep cs f_step no_der no_zin no_zout in
  let reset () ({ state; _ } as st) = f_reset state; st in
-
-  (* horizon *)
-  let horizon { time; _ } =
-    cstate.horizon -. time in
-
+  let horizon { time; _ } = cs.horizon -. time in
  let jump _ = true in
+  let zset = mkzset cs f_step no_der no_zout no_time in
+  let cset = mkcset cs f_step no_der no_zin no_zout no_time in
+  let cget = mkcget cs f_step no_der no_zin no_zout no_time in
+  { state; fder; fzer; step; fout; reset;
+    horizon; cset; cget; zset; zsize; csize; jump }

-  (* 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
+(** Main lifting function. *)
+let lift (f : cstate -> (time * 'a, 'b) node)
+  : (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode def
+= fun () ->
+  let cs = { cvec=cmake 0; dvec=cmake 0; zinvec=zmake 0; zoutvec=cmake 0;
+             cend=0; zend=0; zindex=0; cindex=0; cmax=0; zmax=0; major=false;
+             horizon=max_float } in
+  let Node m = f cs in
+  HNode (lift_inner cs m)

-  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 }
-
-let lift_hsim n =
+(** Lift a simulation (obtained from zeluc with the [-s] flag). *)
+let lift_hsim (n : unit hsimu)
+  : (unit, unit, unit, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode def
+= fun () ->
  let Hsim {
    alloc; step; reset; derivative; crossings; maxsize; horizon; _
  } = n in
-  let s = alloc () in
-  let state = { state = s; input = None; time = 0.0 } in
-  let csize, zsize = maxsize s 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
-  let no_time = -1.0 in
-  reset s;
+  let state = { state=alloc (); input=None; time=0.0 } in
+  let csize, zsize = maxsize state.state in
+  let no_zin, no_zout = zmake zsize, cmake zsize in
+  let no_der, pos = cmake csize, cmake csize in
+  let no_time = -1.0 in reset state.state;
  let fder { state; time; _ } offset () y =
-    derivative state y ignore_der no_roots_in no_roots_out (time +. offset);
-    ignore_der in
+    derivative state y no_der no_zin no_zout (time +. offset);
+    no_der in
  let fzer { state; time; _ } offset () y =
-    crossings state y no_roots_in no_roots_out (time +. offset); no_roots_out in
+    crossings state y no_zin no_zout (time +. offset); no_zout in
  let fout _ _ () _ = () in
  let step { state; time; _ } offset () =
-    step state cstates ignore_der no_roots_in (time +. offset),
+    step state pos no_der no_zin (time +. offset),
    { state; time=time +. offset; input=Some () } in
-  let reset _ ({ state; _ } as st) = reset state; st in
+  let reset () ({ state; _ } as st) = reset state; st in
  let horizon { state; time; _ } = horizon state -. time in
  let jump _ = true in
  let cset ({ state; _ } as st) _ =
-    derivative state cstates ignore_der no_roots_in no_roots_out no_time; st in
+    derivative state pos no_der no_zin no_zout no_time; st in
  let zset ({ state; _ } as st) zinvec =
-    derivative state cstates ignore_der zinvec no_roots_out no_time; st in
+    derivative state pos no_der zinvec no_zout no_time; st in
  let cget { state; _ } =
-    derivative state cstates ignore_der no_roots_in no_roots_out no_time; cstates in
-
+    derivative state pos no_der no_zin no_zout no_time; pos 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
+
+(* Wrappers around the [step] function (for Zelus 2024). *)
+
+let nmkfder (c : cstate_new) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0; c.time <- time;
+  ignore (f state (time +. o, i)); c.dvec
+
+let nmkfzer (c : cstate_new) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0; c.time <- time;
+  ignore (f state (time +. o, i)); c.zoutvec
+
+let nmkfout (c : cstate_new) f der zin zout { state; time; _ } o i y =
+  c.major <- false; c.cvec <- y; c.dvec <- der; c.zinvec <- zin; c.time <- time;
+  c.zoutvec <- zout; c.cindex <- 0; c.zindex <- 0; f state (time +. o, i)
+
+let nmkstep (c : cstate_new) f der zin zout ({ state; time; _ } as s) o i =
+  s.input <- Some i; s.time <- time +. o; c.time <- time;
+  c.major <- true; c.horizon <- infinity; c.zinvec <- zin;
+  c.zoutvec <- zout; c.dvec <- der; c.cindex <- 0; c.zindex <- 0;
+  let o = f state (s.time, i) in o, s
+
+let nmkzset (c : cstate_new) f der zout time ({ state; input; _ } as s) zin =
+  c.major <- false; c.zoutvec <- zout; c.dvec <- der; c.zinvec <- zin;
+  c.cindex <- 0; c.zindex <- 0; ignore (f state (time, Option.get input)); s
+
+let nmkcset (c : cstate_new) f der zin zout time ({ state; input; _ } as st) _ =
+  c.major <- false; c.horizon <- infinity; c.zinvec <- zin; c.zoutvec <- zout;
+  c.dvec <- der; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time, Option.get input)); st
+
+let nmkcget (c : cstate_new) f der zin zout time { state; input; _ } =
+  c.major <- false; c.horizon <- infinity; c.zinvec <- zin; c.zoutvec <- zout;
+  c.dvec <- der; c.cindex <- 0; c.zindex <- 0;
+  ignore (f state (time, Option.get input)); c.cvec
+
+(** Lift the inner record of a node. *)
+let nlift_inner cs (f : ('s, time * 'a, 'b) node_rec)
+  : (('s, 'a) state, unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) hrec
+= let { alloc=f_alloc; step=f_step; reset=f_reset } = f 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 csize, zsize = cs.cmax, cs.zmax in
+  let no_roots_in, no_roots_out = zmake zsize, cmake zsize in
+  let ignore_der = cmake csize in cs.cvec <- cmake csize; 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 fder = nmkfder cs f_step ignore_der no_roots_in no_roots_out in
+  let fzer = nmkfzer cs f_step ignore_der no_roots_in no_roots_out in
+  let fout = nmkfout cs f_step ignore_der no_roots_in no_roots_out in
+  let step = nmkstep cs f_step ignore_der no_roots_in no_roots_out in
  let reset () ({ state; _ } as st) = f_reset state; st in
-
-  (* horizon *)
-  let horizon { time; _ } =
-    cstate.horizon -. time in
-
+  let horizon { time; _ } = cs.horizon -. time in
  let jump _ = true in
+  let zset = nmkzset cs f_step ignore_der no_roots_out no_time in
+  let cset = nmkcset cs f_step ignore_der no_roots_in no_roots_out no_time in
+  let cget = nmkcget cs f_step ignore_der no_roots_in no_roots_out no_time in
+  { state; fder; fzer; step; fout; reset;
+    horizon; cset; cget; zset; zsize; csize; jump }

-  (* 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 nlift (f : Ztypes.cstate_new -> (time * 'a, 'b) node)
+  : (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode def
+= fun () ->
+  let cs = { 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 m = f cs in
+  HNode (nlift_inner cs m)

-  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
+(** Wrap a discrete node into the format expected by Zélus.
+    Resets allocate a fresh node. *)
+let wrap (n : ('p, 'a, 'b) dnode def) : ('a, 'b) node =
+  let alloc () = ref (n ()) in
+  let step s a =
+    let DNode n = !s in
+    let b, state = n.step n.state a in
+    s := DNode { n with state }; b in
+  let reset s = s := n () in
+  Node { alloc; step; reset }

-  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
+(** Unwrap a Zélus node into a discrete node. *)
+let unwrap (n : ('a, 'b) node) : (unit, 'a, 'b) dnode def = fun () ->
+  let Node { alloc; step; reset } = n in
+  let state = alloc () in
+  let step s a = let b = step s a in b, s in 
+  let reset () s = reset s; s in
+  DNode { state; step; reset }

-  HNode
-    { state; fder; fzer; step; fout; reset;
-      horizon; cset; cget; zset; zsize; csize; jump }
+(* let lift_assert (n : ('a, 'b) hnodea) *)
+  (* : (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) hnode_sa def *)
+(* = fun () -> *)
+  (* let NodeA { body; check } = n in *)
+  (* let cs = { 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 } in *)
+  (* let model = lift_inner cs (body cs) in *)
+  (* let proj ({ u; _ } as v) = { v with u=fun t -> (u t).state } in *)
+  (* let proj = Hsim.Utils.sigmap proj in *)
+  (* let check = Hsim.Utils.compose proj (unwrap check) in *)
+  (* HNodeSA { model; check } *)
--- a/src/lib/std/output.ml
+++ b/src/lib/std/output.ml
@ -44,10 +44,10 @@ let print_limits { h; _ } =
  if h <= 0.0 then Format.printf "D: % .10e\n" 0.0
  else Format.printf "C: % .10e\t% .10e\n" 0.0 h

-let print samples n =
-  let DNode m = compose n (compose track (map (print_sample samples))) in
+let print samples n () =
+  let DNode m = compose n (compose track (sigmap (print_sample samples))) () in
  DNode { m with reset=fun p -> m.reset (p, ((), ())) }

-let print_h samples n =
-  let DNode m = compose n (compose track (map (print_sample_h samples))) in
+let print_h samples n () =
+  let DNode m = compose n (compose track (sigmap (print_sample_h samples))) () in
  DNode { m with reset=fun p -> m.reset (p, ((), ())) }
--- a/src/lib/std/runtime.ml
+++ b/src/lib/std/runtime.ml
@ -10,6 +10,7 @@ let opts = ref [
  "-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";
+  "-check", Arg.Set_int Solve.assertion_samples, "\tAssertion checking frequency (default: 0.0)";
 ]

 let anon = ref (fun s -> Format.eprintf "Unexpected argument: %s\n" s; exit 1)
@ -29,7 +30,7 @@ let go
  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 solver model) input output
+  ignore @@ Hsim.Utils.run_on (Solve.build_sim solver model ()) input output

 let go_discrete
  (input : unit -> 'a)
@ -52,4 +53,4 @@ let go_2024
  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
+  ignore @@ Hsim.Utils.run_on (Solve.build_sim_2024 solver model ()) input output
--- a/src/lib/std/solve.ml
+++ b/src/lib/std/solve.ml
@ -27,15 +27,16 @@ let build_sim
  (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
+    (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode def
+= fun () ->
+  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
+     | RK45 -> (fun () -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve ())
+     | Sundials -> Solvers.StatefulSundials.InPlace.csolve)
+    (fun () -> 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
+  let DNode s = Hsim.Utils.(compose (run model solver) track) () in
  DNode { s with reset=fun p -> s.reset (p, ())}

 let build_sim_2024
@ -43,15 +44,16 @@ let build_sim_2024
  (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
+    (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode def
+= fun () ->
+  let model = Lift.nlift 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
+     | RK45 -> (fun () -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve ())
+     | Sundials -> Solvers.StatefulSundials.InPlace.csolve)
+    (fun () -> 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
+  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. *)
@ -59,21 +61,13 @@ 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 }
+= Lift.wrap @@ build_sim solver model

 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 }
+= Lift.wrap @@ build_sim_2024 solver model

 let solve_ode45         m = solve      RK45     m
 let solve_ode45_2024    m = solve_2024 RK45     m
@ -190,20 +184,26 @@ let synchr
  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 iter
+  (n : int)
+  (f : ('a, unit) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let 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 }
+  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 iter_t
+  (n : int)
+  (f : (time * 'a, unit) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let 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 }
+  Node { alloc; step; reset }

 (** Sample a value [n] times and assert [f] on the samples. *)
 let check
@ -214,7 +214,7 @@ let check
    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 })
+  iter_t n (Node { alloc; reset; step })

 (** Sample a value [n] times and assert [f] on the dated samples. *)
 let check_t
@ -226,3 +226,50 @@ let check_t
    with Assert_failure _ ->
      (Format.eprintf "Assertion failed at time %.10e\n" now; exit 1) in
  iter_t n (Ztypes.Node { alloc; reset; step })
+
+let period'
+  (p : float)
+  (Node { alloc; step; reset } : ('a, unit) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let alloc () = ref (0.0, alloc ()) in
+  let step s = Option.iter @@ fun (v, _) ->
+    let offset, st = !s in
+    let l, o = Utils.period ~offset v p in
+    List.iter (fun (_, v) -> step st v) l; s := o, st in
+  let reset s = let _, st = !s in reset st; s := 0.0, st in
+  Node { alloc; step; reset }
+
+let period'_t
+  (p : float)
+  (Node { alloc; step; reset } : (time * 'a, unit) Ztypes.node)
+: ('a signal_t, unit) Ztypes.node
+= let alloc () = ref (0.0, alloc ()) in
+  let step s = Option.iter @@ fun (v, h) ->
+    let offset, st = !s in
+    let l, o = Utils.period ~offset v p in
+    List.iter (fun (t, v) -> step st (t +. h, v)) l; s := o, st in
+  let reset s = let _, st = !s in reset st; s := 0.0, st in
+  Node { alloc; step; reset }
+
+let assertion_samples = ref 100
+
+let build_assertion
+  (solver : solver)
+  (assertion : Ztypes.cstate -> (time * 'a, bool) Ztypes.node)
+: ('a signal, bool) Ztypes.node
+= let n = build_sim solver assertion in
+  let n () =
+    let step st = function
+      | Some i ->
+          let v = ref true in
+          let st = Utils.run_on st i (fun (b, now) ->
+            let l = Utils.sample_tracked (b, now) !assertion_samples in
+            List.iter (fun (_, b) -> v := b && !v) l) in
+          !v, st
+      | None -> true, st in
+    let reset p (DNode n) = DNode { n with state=n.reset p n.state } in
+    DNode { state=n (); step; reset } in
+  Lift.wrap n
+
+let build_assertion_rk45 a     = build_assertion RK45 a
+let build_assertion_sundials a = build_assertion Sundials a
--- a/src/lib/std/solve.mli
+++ b/src/lib/std/solve.mli
@ -1,4 +1,6 @@

+open Hsim.Types
+
 type time = float
 type 'a value = 'a Hsim.Types.value
 type 'a signal = 'a value option
@ -10,13 +12,15 @@ val horizon : 'a value -> time
 val make    : time * (time -> 'a) -> 'a value
 val apply   : 'a value * time -> 'a

+(* val sustain : 'a -> 'a value *)
+
 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
+   'a signal, 'b signal_t) Hsim.Types.dnode def

 val build_sim_2024 :
  solver ->
@ -24,7 +28,7 @@ val build_sim_2024 :
  (unit *
   ((Ztypes.cvec, Ztypes.dvec) Hsim.Solver.ivp *
    (Ztypes.cvec, Ztypes.zoutvec) Hsim.Solver.zc),
-   'a signal, 'b signal_t) Hsim.Types.dnode
+   'a signal, 'b signal_t) Hsim.Types.dnode def

 val solve :
  solver ->
@ -54,8 +58,29 @@ val synchr :
  ('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 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
+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
+
+val period' :
+  float -> ('a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node
+val period'_t :
+  float -> (time * 'a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node
+
+val assertion_samples : int ref
+val build_assertion :
+  solver ->
+  (Ztypes.cstate -> (time * 'a, bool) Ztypes.node) ->
+  ('a signal, bool) Ztypes.node
+val build_assertion_sundials :
+  (Ztypes.cstate -> (time * 'a, bool) Ztypes.node) ->
+  ('a signal, bool) Ztypes.node
+val build_assertion_rk45 :
+  (Ztypes.cstate -> (time * 'a, bool) Ztypes.node) ->
+  ('a signal, bool) Ztypes.node
--- a/src/lib/std/solve.zli
+++ b/src/lib/std/solve.zli
@ -16,8 +16,11 @@ val synchr :
  ('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 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
+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
+
+val period'   : float -S-> ('a        -D-> bool) -S-> 'a signal_t -D-> unit
+val period'_t : float -S-> (time * 'a -D-> bool) -S-> 'a signal_t -D-> unit
--- a/src/lib/std/ztypes.ml
+++ b/src/lib/std/ztypes.ml
@ -25,12 +25,14 @@ type zero = bool

 (* a synchronous stream function with type 'a -D-> 'b *)
 (* is represented by an OCaml value of type ('a, 'b) node *)
+type ('s, 'a, 'b) node_rec = {
+  alloc : unit -> 's; (* allocate the state *)
+  step : 's -> 'a -> 'b; (* compute a step *)
+  reset : 's -> unit; (* reset/inialize the state *)
+}
+
 type ('a, 'b) node =
-    Node:
-      { alloc : unit -> 's; (* allocate the state *)
-        step : 's -> 'a -> 'b; (* compute a step *)
-        reset : 's -> unit; (* reset/inialize the state *)
-      } -> ('a, 'b) node
+    Node: ('s, 'a, 'b) node_rec -> ('a, 'b) node

 (* the same with a method copy *)
 type ('a, 'b) cnode =