From ffc583985ada3463e4f5cb0a796f239ede24f756 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Fri, 11 Jul 2025 11:21:07 +0200 Subject: [PATCH] feat: lift runtime into language, start of zelus 2024 compatibility --- doc/notes.typ | 17 ++ doc/rep.typ | 115 +++++++--- doc/sources.bib | 149 ++++++++++-- exm/builtins/main.ml | 12 +- exm/zelus/{ballz/ballz.zls => ball/ball.zls} | 0 exm/zelus/{ballz => ball}/dune | 4 +- exm/zelus/{ballz => ball}/main.ml | 2 +- exm/zelus/ball/ztypes.ml | 12 + exm/zelus/ballz/ztypes.ml | 21 -- exm/zelus/cradle/cradle.zls | 51 +++++ exm/zelus/cradle/dune | 17 ++ exm/zelus/cradle/format.zli | 2 + exm/zelus/cradle/main.ml | 30 +++ exm/zelus/cradle/ztypes.ml | 12 + exm/zelus/sincos/ztypes.ml | 11 +- exm/zelus/solve/dune | 17 ++ exm/zelus/solve/main.ml | 10 + exm/zelus/solve/solve.zli | 23 ++ exm/zelus/solve/time.zls | 59 +++++ exm/zelus/solve/ztypes.ml | 16 ++ exm/zelus_2024/ball/ball.ml | 58 +++++ exm/zelus_2024/ball/dune | 9 + exm/zelus_2024/ball/main.ml | 7 + exm/zelus_2024/ball/ztypes.ml | 12 + src/lib/hsim/sim.ml | 30 --- src/lib/hsim/types.ml | 4 + src/lib/hsim/utils.ml | 58 ++++- src/lib/solvers/statefulRK45.ml | 6 +- src/lib/solvers/statefulSundials.ml | 6 +- src/lib/solvers/statefulZ.ml | 4 +- src/lib/std/lift.ml | 112 ++++++++- src/lib/std/runtime.ml | 50 ++-- src/lib/std/runtime.mli | 24 ++ src/lib/std/solve.ml | 228 +++++++++++++++++++ src/lib/std/solve.mli | 61 +++++ src/lib/std/solve.zli | 23 ++ src/lib/std/ztypes.ml | 25 +- 37 files changed, 1154 insertions(+), 143 deletions(-) rename exm/zelus/{ballz/ballz.zls => ball/ball.zls} (100%) rename exm/zelus/{ballz => ball}/dune (76%) rename exm/zelus/{ballz => ball}/main.ml (68%) create mode 100644 exm/zelus/ball/ztypes.ml delete mode 100644 exm/zelus/ballz/ztypes.ml create mode 100644 exm/zelus/cradle/cradle.zls create mode 100644 exm/zelus/cradle/dune create mode 100644 exm/zelus/cradle/format.zli create mode 100644 exm/zelus/cradle/main.ml create mode 100644 exm/zelus/cradle/ztypes.ml create mode 100644 exm/zelus/solve/dune create mode 100644 exm/zelus/solve/main.ml create mode 100644 exm/zelus/solve/solve.zli create mode 100644 exm/zelus/solve/time.zls create mode 100644 exm/zelus/solve/ztypes.ml create mode 100644 exm/zelus_2024/ball/ball.ml create mode 100644 exm/zelus_2024/ball/dune create mode 100644 exm/zelus_2024/ball/main.ml create mode 100644 exm/zelus_2024/ball/ztypes.ml create mode 100644 src/lib/std/runtime.mli create mode 100644 src/lib/std/solve.ml create mode 100644 src/lib/std/solve.mli create mode 100644 src/lib/std/solve.zli diff --git a/doc/notes.typ b/doc/notes.typ index ecc9b8f..e65a4a7 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -751,3 +751,20 @@ let hybrid ball () = y where and assert (y >= 0) ``` Is this an issue? + += Ideas + +== Continuous assertions with lifted runtime + +Given a synchronous node simulating a hybrid one : +```zelus +let hybrid sincos () = (sin, cos) where + rec der sin = cos init 0.0 + and der cos = -. sin init 1.0 + +let node sincos_sim = Solve.solve_ode45 sincos +``` +We could create a primitive `Solve.assert_continuous` which would take as input +an `'a value` and a function `'a -> Solve.cond`, where `Solve.cond` is a +zero-crossing expression, for instance `fun v -> Solve.up(-. v)`. This could be +passed along to a zero-crossing solver during continuous phases. diff --git a/doc/rep.typ b/doc/rep.typ index fe4061f..27f2899 100644 --- a/doc/rep.typ +++ b/doc/rep.typ @@ -11,10 +11,12 @@ #let simulink = smallcaps[Simulink] #let sundials = smallcaps[Sundials CVODE] #let zelus = smallcaps[Zélus] -#let TODO(..what) = { +#let note(color, prefix, ..what) = { let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") } - block(fill: red, width: 100%, inset: 5pt, align(center, raw("TODO" + msg))) + block(fill: color, width: 100%, inset: 5pt, align(center, raw(prefix + msg))) } +#let TODO(..what) = note(red, "TODO", ..what) +#let MENTION(..what) = note(gray, "MENTION", ..what) #let adot(s) = $accent(#s, dot)$ #let addot(s) = $accent(#s, dot.double)$ @@ -97,14 +99,14 @@ physical systems. Continuous phases are described using ordinary differential equations (ODEs), while discrete phases can be represented as a reactive program in a synchronous language such as #lustre or #esterel. -As a first example, say we wish to model a bouncing ball. We could start by -describing its distance from the ground $y$ with a second-order differential -equation -$ addot(y) = -9.81 $ +As an illustration, say we wished to model an extensively studied system: a +bouncing ball. We could start by describing its distance from the ground $y$ as +a function of time, with a second-order differential equation +$ addot(y) = -9.81, $ where $addot(y)$ denotes the second order derivative of $y$ with -respect to time (the acceleration of the ball), and $9.81$ is the gravitational -constant $g$: the acceleration of the ball is its negation. We now give the -initial position and speed of the ball: +respect to time $(d^2y)/(d t^2)$ (the acceleration of the ball), and $9.81$ is +the gravitational constant $g$: the acceleration of the ball is its negation. We +now give the initial position and speed of the ball: $ y(0) = 50 space space space adot(y)(0) = 0 $ We have just described an initial value problem: given an ODE and an initial value for its dependent variable, its solution is a function $y(t)$ returning @@ -113,17 +115,20 @@ this function using an ODE solver, such as #sundials. As of right now, our ball will fall until the end of time; we have not said anything about how it bounces when it hits the floor. To do so, we need a notion -of discrete _events_. These are modelled by zero-crossings: we monitor a certain -value and stop when it goes from strictly negative to positive or null. For our -purposes, we choose $-y$ as the monitored value, and call the zero-crossing -event $z$. When $z$ occurs (i.e., when the ball touches the ground), we set the -speed $adot(y)$ to $-k dot "last"(adot(y))$, where $"last"(y)$ denotes the left -limit of $y$ (we cannot specify $adot(y)$ in terms of itself), and $k$ is a -factor modelling the loss of inertia due to the collision (say, $0.8$). We can -then resume the approximation of the solution. +of _events_: we need to identify exactly when the ball hits the ground, so that +we may take action to make it bounce. These events are modelled by +zero-crossings: we monitor a certain value and stop when it goes from strictly +negative to positive or null. For our purposes, we choose $-y$ as the monitored +value, and call the zero-crossing event $z$. When $z$ occurs (i.e., when the +ball touches the ground), we set the speed $adot(y)$ to +$-k dot #raw(lang: "zelus", "last")\(adot(y))$, where +$#raw(lang: "zelus", "last")\(y)$ denotes the left limit of $y$ (we cannot +specify $adot(y)$ in terms of itself), and $k$ is a factor modelling the loss of +inertia due to the collision (say, $0.8$). We can then resume the approximation +of the solution. @lst:ball.zls shows how such a model might be expressed in the concrete syntax -of #zelus. +of #zelus @cit:zelus_sync_lng_with_ode. #figure(placement: top, gap: 2em, ```zelus @@ -135,22 +140,68 @@ of #zelus. caption: [The bouncing ball in #zelus] ) -More formally, a hybrid system can be described as an automaton +More formally, and as done in @cit:alg_ana_hyb_sys, a hybrid system $cal(H)$ can +be defined as a graph whose nodes represent continuous behaviour, and whose +edges represent discrete transitions: +$ cal(H) = (L o c, V a r, E d g, A c t, I n v, I n i) $ +where: +- $L o c$ is a finite set of _locations_; +- $V a r$ is a finite set of _variables_; +- $E d g subset.eq L o c times F times L o c$ is a finite set of _transitions_ + == Executing models -Executing such a model is quite simple. There are two modes of execution: -discrete and continuous. In continuous mode, we call the ODE solver in order to -obtain an approximation of the variables defined through ODEs, and monitor for -zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in which -we perform computation steps as needed, until no other zero-crossing occurs, in -which case we go back to the continuous mode, and repeat, as seen in @automaton. +To execute such a model, we first compile it into a synchronous function, as +described in @cit:sync_based_codegen_hyb_sys_lng. The details of this +compilation step are not particularly relevant to our purposes, and can be +ignored. What is more interesting is the output of this compilation step: a +single synchronous function. The simulation loop is then itself described as a +synchronous function operating on -#figure(finite.automaton( - (D: (D: "cascade", C: "no cascade"), - C: (C: "no zero-crossing", D: "zero-crossing")), - initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3) -), caption: [High-level overview of the runtime], placement: top) +#MENTION("Use of a single solver") + +#pagebreak() + +// The compilation of a hybrid model into a synchronous function is described in +// detail in @cit:sync_based_codegen_hyb_sys_lng, but can be summarized quite +// succintly as follows. By pairing this synchronous function with an +// off-the-shelf ODE solver like #sundials, we can then simulate the dynamics of +// the system. This is done by repeatedly performing execution steps according to +// two different modes: discrete and continuous. + +// The continuous mode operates as follows: we first call the ODE solver in order +// to approximate the dynamics of the model's continuous state. + +// Continuous steps first call the ODE solver to approximate the dynamics of the +// model's continuous variables. The solver will return a function defined on a +// time interval, which we then provide as input to the zero-crossing solver, which +// will monitor the evolution of zero-crossing values along this interval. After +// both solvers have been called, we choose what the next step's mode will be: +// - if no zero-crossings have been detected, we output the entire solution +// provided by the ODE solver, and the next step remains continuous; +// - if a zero-crossing occurs, we return the solution provided by the ODE solver +// up to the zero-crossing instant, and the next step becomes a discrete step. + +// Discrete steps perform state changes and side effects. We first call the model's +// step function, which updates the state and outputs a value. We then decide what +// the next step is. If a zero-crossing event occured due to the current step, the +// next step is another discrete step. If no new event occured, we perform a +// continuous step. + +// Executing such a model is quite simple. There are two modes of execution: +// discrete and continuous. In continuous mode, we call the ODE solver in order +// to obtain an approximation of the variables defined through ODEs, and monitor +// for zero-crossings. If a zero-crossing occurs, we enter the discrete mode, in +// which we perform computation steps as needed, until no other zero-crossing +// occurs, in which case we go back to the continuous mode, and repeat, as seen +// in @automaton. + +// #figure(finite.automaton( +// (D: (D: "cascade", C: "no cascade"), +// C: (C: "no zero-crossing", D: "zero-crossing")), +// initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3) +// ), caption: [High-level overview of the runtime], placement: top) = Runtime To solve this issue, we need to redefine what the runtime of our hybrid system @@ -180,10 +231,10 @@ required by the assertion becomes a state variable. == Solvers as synchronous nodes == Simulations as synchronous nodes -#TODO("talk about the new runtime") +#MENTION("the new runtime") = Assertions -#TODO("talk about how assertions are done") +#MENTION("how assertions are done") #pagebreak() = Annex diff --git a/doc/sources.bib b/doc/sources.bib index afd15fa..4c8f0a4 100644 --- a/doc/sources.bib +++ b/doc/sources.bib @@ -1,16 +1,137 @@ -@article{ - ns_sem_benveniste_bourke_caillaud_pouzet, - title={Non-Standard Semantics of Hybrid Systems Modelers}, - author={Benveniste, Albert and Bourke, Timothy and - Caillaud, Benoıt and Pouzet, Marc}, - year={2011}, - language={en} +@article{cit:nonstd_sem_hyb_sys_mod, + title = {Non-standard semantics of hybrid systems modelers}, + journal = {Journal of Computer and System Sciences}, + volume = {78}, + number = {3}, + pages = {877-910}, + year = {2012}, + note = {In Commemoration of Amir Pnueli}, + issn = {0022-0000}, + doi = {https://doi.org/10.1016/j.jcss.2011.08.009}, + url = {https://www.sciencedirect.com/science/article/pii/S0022000011001061}, + author = {Albert Benveniste and Timothy Bourke and Benoît Caillaud and Marc + Pouzet}, + keywords = {Hybrid systems, Hybrid systems modelers, Non-standard analysis, + Non-standard semantics, Constructive semantics, Kahn process + networks, Compilation of hybrid systems}, + abstract = {Hybrid system modelers have become a corner stone of complex + embedded system development. Embedded systems include not only + control components and software, but also physical devices. In + this area, Simulink is a de facto standard design framework, and + Modelica a new player. However, such tools raise several issues + related to the lack of reproducibility of simulations + (sensitivity to simulation parameters and to the choice of a + simulation engine). In this paper we propose using techniques + from non-standard analysis to define a semantic domain for hybrid + systems. Non-standard analysis is an extension of classical + analysis in which infinitesimal (the ε and η in the celebrated + generic sentence ∀ε∃η… of college maths) can be manipulated as + first class citizens. This approach allows us to define both a + denotational semantics, a constructive semantics, and a Kahn + Process Network semantics for hybrid systems, thus establishing + simulation engines on a sound but flexible mathematical + foundation. These semantics offer a clear distinction between the + concerns of the numerical analyst (solving differential + equations) and those of the computer scientist (generating + execution schemes). We also discuss a number of practical and + fundamental issues in hybrid system modelers that give rise to + non-reproducibility of results, non-determinism, and undesirable + side effects. Of particular importance are cascaded mode changes + (also called “zero-crossings” in the context of hybrid systems + modelers).}, } -@inbook{ - opsem_lee_zheng, - title={Operational Semantics of Hybrid Systems}, - ISBN={978-3-540-25108-8}, - author={Lee, Edward A. and Zheng, Haiyang}, - year={2005}, - language={en} +@inbook{cit:op_sem_hyb_sys, + address = {Berlin, Heidelberg}, + series = {Lecture Notes in Computer Science}, + title = {Operational Semantics of Hybrid Systems}, + volume = {3414}, + ISBN = {978-3-540-25108-8}, + url = {http://link.springer.com/10.1007/978-3-540-31954-2_2}, + DOI = {10.1007/978-3-540-31954-2_2}, + abstractNote = {This paper discusses an interpretation of hybrid systems as + executable models. A specification of a hybrid system for this + purpose can be viewed as a program in a domain-specific + programming language. We describe the semantics of HyVisual, + which is such a domain-specific programming language. The + semantic properties of such a language affect our ability to + understand, execute, and analyze a model. We discuss several + semantic issues that come in defining such a programming + language, such as the interpretation of discontinuities in + continuous-time signals, and the interpretation of + discrete-event signals in hybrid systems, and the + consequences of numerical ODE solver techniques. We describe + the solution in HyVisual by giving its operational semantics. + }, + booktitle = {Hybrid Systems: Computation and Control}, + publisher = {Springer Berlin Heidelberg}, + author = {Lee, Edward A. and Zheng, Haiyang}, + editor = {Morari, Manfred and Thiele, Lothar}, + year = {2005}, + pages = {25–53}, + collection = {Lecture Notes in Computer Science}, + language = {en}, +} +@inproceedings{cit:zelus_sync_lng_with_ode, + address = {Philadelphia Pennsylvania USA}, + title = {Zélus: a synchronous language with ODEs}, + ISBN = {978-1-4503-1567-8}, + url = {https://dl.acm.org/doi/10.1145/2461328.2461348}, + DOI = {10.1145/2461328.2461348}, + abstractNote = { Z´elus is a new programming language for modeling systems + that mix discrete logical time and continuous time behaviors. + From a user’s perspective, its main originality is to extend + an existing Lustre-like synchronous language with Ordinary + Differential Equations (ODEs). The extension is conservative: + any synchronous program expressed as dataflow equations and + hierarchical automata can be composed arbitrarily with ODEs + in the same source code. }, + booktitle = { Proceedings of the 16th international conference on Hybrid + systems: computation and control }, + publisher = {ACM}, + author = {Bourke, Timothy and Pouzet, Marc}, + year = {2013}, + month = apr, + pages = {113–118}, + language = {en}, +} +@inbook{cit:sync_based_codegen_hyb_sys_lng, + address = {Berlin, Heidelberg}, + series = {Lecture Notes in Computer Science}, + title = {A Synchronous-Based Code Generator for Explicit Hybrid Systems + Languages}, + volume = {9031}, + rights = {http://www.springer.com/tdm}, + ISBN = {978-3-662-46662-9}, + url = {http://link.springer.com/10.1007/978-3-662-46663-6_4}, + DOI = {10.1007/978-3-662-46663-6_4}, + abstractNote = {Modeling languages for hybrid systems are cornerstones of + embedded systems development in which software interacts with + a physical environment. Sequential code generation from such + languages is important for simulation efficiency and for + producing code for embedded targets. Despite being routinely + used in industrial compilers, code generation is rarely, if + ever, described in full detail, much less formalized. Yet + formalization is an essential step in building trustable + compilers for critical embedded software development.}, + booktitle = {Compiler Construction}, + publisher = {Springer Berlin Heidelberg}, + author = {Bourke, Timothy and Colaço, Jean-Louis and Pagano, Bruno and + Pasteur, Cédric and Pouzet, Marc}, + editor = {Franke, Björn}, + year = {2015}, + pages = {69–88}, + collection = {Lecture Notes in Computer Science}, + language = {en}, +} +@article{cit:alg_ana_hyb_sys, + title = {The algorithmic analysis of hybrid systems}, + author = {Alur, Rajeev and Courcoubetis, Costas and Halbwachs, Nicolas and + Henzinger, Thomas A and Ho, P-H and Nicollin, Xavier and Olivero, + Alfredo and Sifakis, Joseph and Yovine, Sergio}, + journal = {Theoretical computer science}, + volume = {138}, + number = {1}, + pages = {3--34}, + year = {1995}, + publisher = {Elsevier}, } diff --git a/exm/builtins/main.ml b/exm/builtins/main.ml index 43520dd..95e06a8 100644 --- a/exm/builtins/main.ml +++ b/exm/builtins/main.ml @@ -73,21 +73,21 @@ let output = let sim = if !sundials then let open StatefulSundials in - let c = if !inplace then InPlace.csolve else Functional.csolve in + let c = if !inplace then InPlace.csolve () else Functional.csolve () in let open StatefulZ in - let z = if !inplace then InPlace.zsolve else Functional.zsolve in + let z = if !inplace then InPlace.zsolve () else Functional.zsolve () in let s = Solver.solver c (d_of_dc z) in let open Sim.Sim(val st) in - run_until_n (output !sample (run m s)) + Hsim.Utils.run_until_n (output !sample (run m s)) else let open StatefulRK45 in - let c = if !inplace then InPlace.csolve else Functional.csolve in + let c = if !inplace then InPlace.csolve () else Functional.csolve () in let open StatefulZ in - let z = if !inplace then InPlace.zsolve else Functional.zsolve in + let z = if !inplace then InPlace.zsolve () else Functional.zsolve () in let s = Solver.solver_c c z in let open Sim.Sim(val st) in let n = if !accel then accelerate m s else run m (d_of_dc s) in - run_until_n (output !sample n) + Hsim.Utils.run_until_n (output !sample n) let () = sim !stop !steps ignore diff --git a/exm/zelus/ballz/ballz.zls b/exm/zelus/ball/ball.zls similarity index 100% rename from exm/zelus/ballz/ballz.zls rename to exm/zelus/ball/ball.zls diff --git a/exm/zelus/ballz/dune b/exm/zelus/ball/dune similarity index 76% rename from exm/zelus/ballz/dune rename to exm/zelus/ball/dune index 2718a9b..87598d9 100644 --- a/exm/zelus/ballz/dune +++ b/exm/zelus/ball/dune @@ -4,9 +4,9 @@ (:standard -w -a)))) (rule - (targets ballz.ml ballz.zci) + (targets ball.ml ball.zci) (deps - (:zl ballz.zls)) + (:zl ball.zls)) (action (run zeluc %{zl}))) diff --git a/exm/zelus/ballz/main.ml b/exm/zelus/ball/main.ml similarity index 68% rename from exm/zelus/ballz/main.ml rename to exm/zelus/ball/main.ml index 14ad133..85d403e 100644 --- a/exm/zelus/ballz/main.ml +++ b/exm/zelus/ball/main.ml @@ -3,4 +3,4 @@ open Std let input _ = () let output (now, (y, _, _)) = Format.printf "%.10e\t%.10e\n" now y -let () = Runtime.go input Ballz.ball output +let () = Runtime.go input Ball.ball output diff --git a/exm/zelus/ball/ztypes.ml b/exm/zelus/ball/ztypes.ml new file mode 100644 index 0000000..f53e331 --- /dev/null +++ b/exm/zelus/ball/ztypes.ml @@ -0,0 +1,12 @@ +include Std +include Ztypes +include Solvers + +module type IGNORE = sig end +module Defaultsolver : IGNORE = struct end + +module Zlsrun = struct + module Make (S : IGNORE) = struct + let go _ = () + end +end diff --git a/exm/zelus/ballz/ztypes.ml b/exm/zelus/ballz/ztypes.ml deleted file mode 100644 index 7dcf58c..0000000 --- a/exm/zelus/ballz/ztypes.ml +++ /dev/null @@ -1,21 +0,0 @@ -include Std -include Ztypes -include Solvers - -module type IGNORE = sig end -module Defaultsolver : IGNORE = struct end - -module Zlsrun = struct - module Make (S : IGNORE) = struct - let go s = - let s = Lift.lift_hsim s in - let open Hsim in - let state = (module State.InPlaceSimState : State.SimState) in - let solver = - Solver.solver (StatefulSundials.InPlace.csolve) - (Types.d_of_dc StatefulZ.InPlace.zsolve) in - let open Sim.Sim(val state) in - () - (* run_until_n (Utils.ignore 0 (run s solver)) 30. 1 ignore *) - end -end diff --git a/exm/zelus/cradle/cradle.zls b/exm/zelus/cradle/cradle.zls new file mode 100644 index 0000000..8da18a7 --- /dev/null +++ b/exm/zelus/cradle/cradle.zls @@ -0,0 +1,51 @@ + +let mp6 = -. (3.1416 /. 6.) +let g = 9.80665 +let l = 0.2 + +let pi0 = mp6 +let pi1 = 0. +let pi2 = 0. + +let acc x = -. g /. l *. (sin x) + +let hybrid cradle2() = + let rec der p0 = v0 init pi0 reset h01 -> last p1 + and der v0 = acc(p0) init 0.0 reset h01 -> last v1 + and der p1 = v1 init pi1 reset h01 -> last p0 + and der v1 = acc(p1) init 0.0 reset h01 -> last v0 + and h01 = up(last p0 -. last p1) + and init h = -0.1 + and present h01 -> do h = -1.0 *. last h done + else do der h = 0.0 done + in (h, (p0, v0 /. 10.) , (p1, v1 /. 10.)) + +let hybrid cradle3() = + let rec der p0 = v0 init pi0 reset h01 -> last p1 + and der v0 = acc(p0) init 0.0 reset h01 -> last v1 + and der p1 = v1 init pi1 reset h01 -> last p0 | h12 -> last p2 + and der v1 = acc(p1) init 0.0 reset h01 -> last v0 | h12 -> last v2 + and der p2 = v2 init pi2 reset h12 -> last p1 + and der v2 = acc(p2) init 0.0 reset h12 -> last v1 + and h01 = up(last p0 -. last p1) + and h12 = up(last p1 -. last p2) + and init h1 = -0.1 + and present h01 -> do h1 = -1.0 *. last h1 done else do der h1 = 0.0 done + and init h2 = -0.1 + and present h12 -> do h2 = -1.0 *. last h2 done else do der h2 = 0.0 done + in (p0, p1, p2, h1, h2) + +let node print(v, s) = + Format.printf "% .10e%s" v s + +let hybrid main() = + let der t = 1.0 init 0.0 in + let (p0, p1, p2, h1, h2) = cradle3() in + present (period(0.05)) -> ( + print(t, "\t"); + print(p0, "\t"); + print(p1, "\t"); + print(p2, "\t"); + print(h1, "\t"); + print(h2, "\n") + ); () diff --git a/exm/zelus/cradle/dune b/exm/zelus/cradle/dune new file mode 100644 index 0000000..f49eee8 --- /dev/null +++ b/exm/zelus/cradle/dune @@ -0,0 +1,17 @@ +(env + (dev + (flags + (:standard -w -a)))) + +(rule + (targets cradle.ml cradle.zci format.zci) + (deps + (:zl cradle.zls) + (:zli format.zli)) + (action + (run zeluc %{zli} %{zl}))) + +(executable + (public_name cradle.exe) + (name main) + (libraries std)) diff --git a/exm/zelus/cradle/format.zli b/exm/zelus/cradle/format.zli new file mode 100644 index 0000000..c66577b --- /dev/null +++ b/exm/zelus/cradle/format.zli @@ -0,0 +1,2 @@ + +val printf : string -> float -> string -> unit diff --git a/exm/zelus/cradle/main.ml b/exm/zelus/cradle/main.ml new file mode 100644 index 0000000..c82080e --- /dev/null +++ b/exm/zelus/cradle/main.ml @@ -0,0 +1,30 @@ + +open Std + +let input2 _ = () +let output2 (now, (h, (p0, v0), (p1, v1))) = + Format.printf "%.10e\t%.10e\t%.10e\n" now p0 p1 + +let input3 _ = () +let output3 (now, (p0, p1, p2, h1, h2)) = + Format.printf "%.10e\t%.10e\t%.10e\t%.10e\t%.10e\t%.10e\n" + now p0 (p1 +. 1.0) (p2 +. 2.0) (h1 +. 3.0) (h2 +. 4.0) + +let input_main _ = () +let output_main (now, ()) = () + +let three = ref false +let main = ref false + +let toggle y n () = + y := true; + List.iter (fun n -> n := false) n + +let () = + Runtime.register_args [ + "-three", Arg.Unit (toggle three [main]), "\tUse the third model"; + "-main", Arg.Unit (toggle main [three]), "\tUse the main model"; + ]; + if !main then Runtime.go input_main Cradle.main output_main + else if !three then Runtime.go input3 Cradle.cradle3 output3 + else Runtime.go input2 Cradle.cradle2 output2 diff --git a/exm/zelus/cradle/ztypes.ml b/exm/zelus/cradle/ztypes.ml new file mode 100644 index 0000000..f53e331 --- /dev/null +++ b/exm/zelus/cradle/ztypes.ml @@ -0,0 +1,12 @@ +include Std +include Ztypes +include Solvers + +module type IGNORE = sig end +module Defaultsolver : IGNORE = struct end + +module Zlsrun = struct + module Make (S : IGNORE) = struct + let go _ = () + end +end diff --git a/exm/zelus/sincos/ztypes.ml b/exm/zelus/sincos/ztypes.ml index 7dcf58c..f53e331 100644 --- a/exm/zelus/sincos/ztypes.ml +++ b/exm/zelus/sincos/ztypes.ml @@ -7,15 +7,6 @@ module Defaultsolver : IGNORE = struct end module Zlsrun = struct module Make (S : IGNORE) = struct - let go s = - let s = Lift.lift_hsim s in - let open Hsim in - let state = (module State.InPlaceSimState : State.SimState) in - let solver = - Solver.solver (StatefulSundials.InPlace.csolve) - (Types.d_of_dc StatefulZ.InPlace.zsolve) in - let open Sim.Sim(val state) in - () - (* run_until_n (Utils.ignore 0 (run s solver)) 30. 1 ignore *) + let go _ = () end end diff --git a/exm/zelus/solve/dune b/exm/zelus/solve/dune new file mode 100644 index 0000000..7fd0b6e --- /dev/null +++ b/exm/zelus/solve/dune @@ -0,0 +1,17 @@ +(env + (dev + (flags + (:standard -w -a)))) + +(rule + (targets time.ml time.zci) + (deps + (:zl time.zls) + (:zli solve.zli)) + (action + (run zeluc %{zli} %{zl}))) + +(executable + (public_name time.exe) + (name main) + (libraries std)) diff --git a/exm/zelus/solve/main.ml b/exm/zelus/solve/main.ml new file mode 100644 index 0000000..e519cf7 --- /dev/null +++ b/exm/zelus/solve/main.ml @@ -0,0 +1,10 @@ + +open Std + +let input () = () +let output () = flush stdout + +let () = + Runtime.parse_args (); + Runtime.go_discrete input Time.main output + diff --git a/exm/zelus/solve/solve.zli b/exm/zelus/solve/solve.zli new file mode 100644 index 0000000..d53fee6 --- /dev/null +++ b/exm/zelus/solve/solve.zli @@ -0,0 +1,23 @@ + +type time = float +type 'a value +type 'a signal = 'a value option +type 'a signal_t = ('a value * time) option + +val horizon : 'a value -> time +val make : time * (time -> 'a) -> 'a value +val apply : 'a value * time -> 'a + +val solve_ode45 : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t +val solve_sundials : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t + +val synchr : + ('a signal -D-> 'b signal_t) -S-> + ('a signal -D-> 'c signal_t) -S-> + 'a signal -D-> ('b * 'c) signal_t + +val iter : int -S-> ('a -D-> unit) -S-> 'a signal_t -D-> unit +val iter_t : int -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit + +val check : int -S-> ('a -D-> bool) -S-> 'a signal_t -D-> unit +val check_t : int -S-> (time * 'a -D-> bool) -S-> 'a signal_t -D-> unit diff --git a/exm/zelus/solve/time.zls b/exm/zelus/solve/time.zls new file mode 100644 index 0000000..9a0c32d --- /dev/null +++ b/exm/zelus/solve/time.zls @@ -0,0 +1,59 @@ + +let epsilon = 0.0001 + +let input _ = () + +let hybrid sincos() = + let rec der sin = cos init 0.0 + and der cos = -. sin init 1.0 + in (sin, cos) + +let sincos_ode45 = Solve.solve_ode45(sincos) +let sincos_sundials = Solve.solve_sundials(sincos) +let sincos_both = Solve.synchr(sincos_ode45)(sincos_sundials) + +let hybrid ball () = + let rec der y = y' init 50.0 reset z -> 0.0 + and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y') + and z = up(-. y) + in y + +let ball_ode45 = Solve.solve_ode45(ball) +let ball_sundials = Solve.solve_sundials(ball) +let ball_both = Solve.synchr(ball_ode45)(ball_sundials) + +let node print_ball_both (now, (y1, y2)) = + print_float(now); print_string("\t"); + print_float(y1); print_string("\t"); + print_float(y2); print_string("\n"); + () + +let node print_sincos (now, (sin, cos)) = + print_float now; print_string "\t"; + print_float sin; print_string "\t"; + print_float cos; print_string "\n" + +let node print_sincos2 (now, ((sin1, cos1), (sin2, cos2))) = + print_float now; print_string "\t"; + print_float sin1; print_string "\t"; + print_float sin2; print_string "\t"; + print_float cos1; print_string "\t"; + print_float cos2; print_string "\n" + +let node check_sincos (now, (sin, cos)) = + print_sincos (now, (sin, cos)); + sin <= 1.0 +. epsilon && sin >= -1.0 -. epsilon && + cos <= 1.0 +. epsilon && cos >= -1.0 -. epsilon + +let node check_sincos2 (now, ((sin1, cos1), (sin2, cos2))) = + print_sincos2 (now, ((sin1, cos1), (sin2, cos2))); + sin1 <= 1.0 +. epsilon && sin1 >= -1.0 -. epsilon && + cos1 <= 1.0 +. epsilon && cos1 >= -1.0 -. epsilon && + sin2 <= 1.0 +. epsilon && sin2 >= -1.0 -. epsilon && + cos2 <= 1.0 +. epsilon && cos2 >= -1.0 -. epsilon + +let node main() = + let input = Some (Solve.make (30.0, input)) fby None in + let o = run sincos_sundials input in + Solve.check_t 100 check_sincos o + diff --git a/exm/zelus/solve/ztypes.ml b/exm/zelus/solve/ztypes.ml new file mode 100644 index 0000000..ef5ae18 --- /dev/null +++ b/exm/zelus/solve/ztypes.ml @@ -0,0 +1,16 @@ +include Std +include Ztypes +include Solvers + +module type IGNORE = sig end +module Defaultsolver : IGNORE = struct end + +module Zlsrun = struct + module Make (S : IGNORE) = struct + let go _ = () + end +end + +module Stdlib = struct + type nonrec 'a option = 'a option +end diff --git a/exm/zelus_2024/ball/ball.ml b/exm/zelus_2024/ball/ball.ml new file mode 100644 index 0000000..b60cdf3 --- /dev/null +++ b/exm/zelus_2024/ball/ball.ml @@ -0,0 +1,58 @@ +(* The Zelus compiler, version 2024-dev + (2025-06-4-15:49) *) +open Ztypes + +type ('e, 'd, 'c, 'b, 'a) ball = + { mutable time: 'e; mutable major: 'd; mutable up: 'c; + mutable y': 'b; mutable y: 'a } + +let ball = + let machine cstate = + let alloc _ = + cstate.cmax <- cstate.cmax + 1; + cstate.zmax <- cstate.zmax + 1; + { time = -1.; + major = false; + up = { zin = false; zout = 1. }; + y' = -1.; + y = { pos = -1.; der = 0. }; + } in + let step self _ = + let cindex = cstate.cindex in + let cpos = ref cindex in + let zindex = cstate.zindex in + let zpos = ref zindex in + cstate.cindex <- cstate.cindex + 1; + cstate.zindex <- cstate.zindex + 1; + self.major <- cstate.major; + self.time <- cstate.time; + if cstate.major then + for i = cindex to 0 do Zls.set cstate.dvec i 0. done + else begin + self.y.pos <- Zls.get cstate.cvec !cpos; + cpos := !cpos + 1 + end; + let result = + self.up.zout <- -. self.y.pos; + if self.up.zin then self.y' <- -0.8 *. self.y'; + self.y.der <- self.y'; + self.y.pos, self.y', self.up.zin in + cpos := cindex; + if cstate.major then begin + Zls.set cstate.cvec !cpos self.y.pos; + cpos := !cpos + 1; + self.up.zin <- false + end else begin + self.up.zin <- Zls.get_zin cstate.zinvec !zpos; + zpos := !zpos + 1 + end; + zpos := zindex; + Zls.set cstate.zoutvec !zpos self.up.zout; + zpos := !zpos + 1; + Zls.set cstate.dvec !cpos self.y.der; + cpos := !cpos + 1; + result in + let reset self = + self.y.pos <- 50.; self.y' <- 0. in + Node { alloc; step; reset } in + machine diff --git a/exm/zelus_2024/ball/dune b/exm/zelus_2024/ball/dune new file mode 100644 index 0000000..fe72b56 --- /dev/null +++ b/exm/zelus_2024/ball/dune @@ -0,0 +1,9 @@ +(env + (dev + (flags + (:standard -w -a)))) + +(executable + (public_name newball.exe) + (name main) + (libraries std)) diff --git a/exm/zelus_2024/ball/main.ml b/exm/zelus_2024/ball/main.ml new file mode 100644 index 0000000..9413200 --- /dev/null +++ b/exm/zelus_2024/ball/main.ml @@ -0,0 +1,7 @@ + +open Std + +let input _ = () +let output (now, (y, _, _)) = Format.printf "%.10e\t%.10e\n" now y +let () = Runtime.go_2024 input Ball.ball output + diff --git a/exm/zelus_2024/ball/ztypes.ml b/exm/zelus_2024/ball/ztypes.ml new file mode 100644 index 0000000..f53e331 --- /dev/null +++ b/exm/zelus_2024/ball/ztypes.ml @@ -0,0 +1,12 @@ +include Std +include Ztypes +include Solvers + +module type IGNORE = sig end +module Defaultsolver : IGNORE = struct end + +module Zlsrun = struct + module Make (S : IGNORE) = struct + let go _ = () + end +end diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 2662a7e..5ab2dc9 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -171,34 +171,4 @@ module Sim (S : SimState) = update ms ss (set_idle st) in DNode { state; step; reset } - - (** Run the model on the given input until the end of the input or until the - model stops answering. *) - let run_on (DNode n) input use = - let out = n.step n.state (Some input) in - let state = match out with None, s -> s | Some o, s -> use o; s in - let rec loop state = - let o, state = n.step state None in - match o with None -> () | Some o -> use o; loop state in - loop state - - (** Run the model on multiple inputs. *) - let run_on_n (DNode n) inputs use = - ignore @@ List.fold_left (fun state i -> - let o, state = n.step state (Some i) in - begin match o with None -> () | Some o -> use o end; - let rec loop state = - let o, state = n.step state None in - match o with None -> state | Some o -> use o; loop state in - loop state) n.state inputs - - (** Run the model autonomously until [h], or until the model stops - answering. *) - let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () } - - (** Run the model autonomously until [h], split in [k] steps. *) - let run_until_n n h k = - let h = h /. float_of_int k in - run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () })) - end diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index aac5cd4..0f0e589 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -14,6 +14,10 @@ type 'a value = - [u: [0, h] -> α] *) type 'a signal = 'a value option +(** A time signal with absolute timestamps added. + These represent the starting date for the value. *) +type 'a signal_t = ('a value * time) option + type ('s, 'p, 'a, 'b) drec = { state : 's; step : 's -> 'a -> 'b * 's; diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml index e7b8190..38c555a 100644 --- a/src/lib/hsim/utils.ml +++ b/src/lib/hsim/utils.ml @@ -7,6 +7,34 @@ let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v } let offset (u : time -> 'a) (now : time) : time -> 'a = fun t -> u (t +. now) +(** Cut a value into two at a specified timestamp. *) +let cutoff { h; u; c } t = + if t < 0.0 then + raise (Invalid_argument "Cutoff point cannot be negative"); + if t > h then + raise (Invalid_argument "Cutoff point cannot be greater than horizon"); + { h=t; c=Continuous; u }, { h=h -. t; c; u=fun n -> u (t +. n) } + +(** Join two values. *) +let join { h=hl; u=ul; c=cl } { h=hr; u=ur; c=cr } = + let h = min hl hr in + let u = fun t -> ul t, ur t in + let c = match cl, cr with + | Continuous, Continuous -> Continuous + | _ -> Discontinuous in + { h; u; c } + +(** Map a function. *) +let map_value f ({ u; _ } as v) = + { v with u=fun t -> f (u t) } + +(** Swap a pair. *) +let swap v = map_value (fun (a, b) -> b, a) v + +let map_signal f v = Option.map (map_value f) v + +let swap_signal v = Option.map swap v + (** Concatenate functions. *) let rec concat = function | [] -> raise (Invalid_argument "Cannot concatenate an empty value list") @@ -67,7 +95,7 @@ let compose_sim DNode { state; step; reset } (** Track the evolution of a signal in time. *) -let track : (unit, 'a signal, ('a value * time) option) dnode = +let track : (unit, 'a signal, 'a signal_t) dnode = let state = 0.0 in let step now = function | None -> None, now @@ -101,3 +129,31 @@ let do_and_reset (DNode m) (DNode n) f = m.reset ms mp, n.reset ns np in DNode { state; step; reset } +(** Run a model on the given input until the end of the input or until the model + stops answering. *) +let run_on (DNode n) input use = + let out = n.step n.state (Some input) in + let state = match out with None, s -> s | Some o, s -> use o; s in + let rec loop state = + let o, state = n.step state None in + match o with None -> () | Some o -> use o; loop state in + loop state + +(** Run the model on multiple inputs. *) +let run_on_n (DNode n) inputs use = + Stdlib.ignore @@ List.fold_left (fun state i -> + let o, state = n.step state (Some i) in + begin match o with None -> () | Some o -> use o end; + let rec loop state = + let o, state = n.step state None in + match o with None -> state | Some o -> use o; loop state in + loop state) n.state inputs + +(** Run the model autonomously until [h], or until the model stops answering. *) +let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () } + +(** Run the model autonomously until [h], split in [k] steps. *) +let run_until_n n h k = + let h = h /. float_of_int k in + run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () })) + diff --git a/src/lib/solvers/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml index d85658e..48a2069 100644 --- a/src/lib/solvers/statefulRK45.ml +++ b/src/lib/solvers/statefulRK45.ml @@ -7,7 +7,8 @@ module Functional = struct type ('state, 'vec) state = { state: 'state; vec: 'vec } - let csolve : (carray, carray) csolver_c = + let csolve () : (carray, carray) csolver_c = + Common.Debug.print "Instantiating RK45"; let open Odexx.Ode45 in let state = @@ -37,7 +38,8 @@ module InPlace = struct type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec } - let csolve : (carray, carray) csolver_c = + let csolve () : (carray, carray) csolver_c = + Common.Debug.print "Instantiating RK45"; let open Odexx.Ode45 in let state = diff --git a/src/lib/solvers/statefulSundials.ml b/src/lib/solvers/statefulSundials.ml index 30156ab..e733da6 100644 --- a/src/lib/solvers/statefulSundials.ml +++ b/src/lib/solvers/statefulSundials.ml @@ -7,7 +7,8 @@ module Functional = struct type ('state, 'vec) state = { state : 'state; vec : 'vec } - let csolve : (carray, carray) csolver = + let csolve () : (carray, carray) csolver = + Format.printf "Instantiating Sundials"; let open Cvode in let state = @@ -37,7 +38,8 @@ module InPlace = struct type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec } - let csolve : (carray, carray) csolver = + let csolve () : (carray, carray) csolver = + Common.Debug.print "Instantiating Sundials"; let open Cvode in let state = diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml index b642a08..aed275f 100644 --- a/src/lib/solvers/statefulZ.ml +++ b/src/lib/solvers/statefulZ.ml @@ -7,7 +7,7 @@ module Functional = struct type ('state, 'vec) state = { state: 'state; vec: 'vec } - let zsolve : (carray, zarray, carray) zsolver_c = + let zsolve () : (carray, zarray, carray) zsolver_c = let open Illinois in let state = @@ -38,7 +38,7 @@ module InPlace = struct type ('state, 'vec) state = { mutable state : 'state; mutable vec : 'vec } - let zsolve : (carray, zarray, carray) zsolver_c = + let zsolve () : (carray, zarray, carray) zsolver_c = let open Illinois in let state = diff --git a/src/lib/std/lift.ml b/src/lib/std/lift.ml index 222b7a6..ee72f9e 100644 --- a/src/lib/std/lift.ml +++ b/src/lib/std/lift.ml @@ -6,8 +6,10 @@ open Ztypes type ('s, 'a) state = { mutable state : 's; mutable input : 'a option; mutable time : time } -let lift f = - let cstate = +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; @@ -61,7 +63,7 @@ let lift f = let o = f_step state (st.time, input) in o, st in - let reset _ ({ state; _ } as st) = f_reset state; st in + let reset () ({ state; _ } as st) = f_reset state; st in (* horizon *) let horizon { time; _ } = @@ -140,3 +142,107 @@ let lift_hsim n = derivative state cstates ignore_der no_roots_in no_roots_out no_time; cstates in HNode { state; fder; fzer; fout; step; reset; horizon; jump; cget; cset; zset; csize; zsize } + +let lift_2024 + (f : Ztypes.cstate_new -> (time * 'a, 'b) node) +: (unit, 'a, 'b, cvec, dvec, zinvec, zoutvec) Hsim.Types.hnode += let cstate = + { cvec = cmake 0; dvec = cmake 0; cindex = 0; zindex = 0; + cend = 0; zend = 0; cmax = 0; zmax = 0; + zinvec = zmake 0; zoutvec = cmake 0; + major = false; horizon = max_float; time=0.0 } in + let Node { alloc=f_alloc; step=f_step; reset=f_reset } = f cstate in + let state = { state = f_alloc (); input = None; time = 0.0 } in + let csize, zsize = cstate.cmax, cstate.zmax in + let no_roots_in = zmake zsize in + let no_roots_out = cmake zsize in + let ignore_der = cmake csize in + let cstates = cmake csize in + cstate.cvec <- cstates; + f_reset state.state; + + let no_time = -1.0 in + + (* the function that compute the derivatives *) + let fder { state; time; _ } offset input y = + cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der; + cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out; + cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time; + ignore (f_step state (time +. offset, input)); + cstate.dvec in + + (* the function that compute the zero-crossings *) + let fzer { state; time; _ } offset input y = + cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der; + cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out; + cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time; + ignore (f_step state (time +. offset, input)); + cstate.zoutvec in + + (* the function which compute the output during integration *) + let fout { state; time; _ } offset input y = + cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der; + cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out; + cstate.cindex <- 0; cstate.zindex <- 0; cstate.time <- time; + f_step state (time +. offset, input) in + + (* the function which compute a discrete step *) + let step ({ state; time; _ } as st) offset input = + st.input <- Some input; + st.time <- time +. offset; + cstate.time <- time; + cstate.major <- true; + cstate.horizon <- infinity; + cstate.zinvec <- no_roots_in; + cstate.zoutvec <- no_roots_out; + cstate.dvec <- ignore_der; + cstate.cindex <- 0; + cstate.zindex <- 0; + let o = f_step state (st.time, input) in + o, st in + + let reset () ({ state; _ } as st) = f_reset state; st in + + (* horizon *) + let horizon { time; _ } = + cstate.horizon -. time in + + let jump _ = true in + + (* the function which sets the zinvector into the *) + (* internal zero-crossing variables *) + let zset ({ state; input; _ } as st) zinvec = + cstate.major <- false; + cstate.zoutvec <- no_roots_out; + cstate.dvec <- ignore_der; + cstate.zinvec <- zinvec; + cstate.cindex <- 0; + cstate.zindex <- 0; + ignore (f_step state (no_time, Option.get input)); + st in + + let cset ({ state; input; _ } as st) _ = + cstate.major <- false; + cstate.horizon <- infinity; + cstate.zinvec <- no_roots_in; + cstate.zoutvec <- no_roots_out; + cstate.dvec <- ignore_der; + cstate.cindex <- 0; + cstate.zindex <- 0; + ignore (f_step state (no_time, Option.get input)); + st in + + let cget { state; input; _ } = + cstate.major <- false; + cstate.horizon <- infinity; + cstate.zinvec <- no_roots_in; + cstate.zoutvec <- no_roots_out; + cstate.dvec <- ignore_der; + cstate.cindex <- 0; + cstate.zindex <- 0; + ignore (f_step state (no_time, Option.get input)); + cstate.cvec in + + HNode + { state; fder; fzer; step; fout; reset; + horizon; cset; cget; zset; zsize; csize; jump } diff --git a/src/lib/std/runtime.ml b/src/lib/std/runtime.ml index 8ebe08e..b998de3 100644 --- a/src/lib/std/runtime.ml +++ b/src/lib/std/runtime.ml @@ -1,31 +1,55 @@ open Hsim.Types -let sample = ref 0 -let stop = ref 10.0 +let sample = ref 0 +let stop = ref 10.0 +let sundials = ref false -let options = [ +let opts = ref [ "-sample", Arg.Set_int sample, "\tSampling frequency (default=0)"; "-stop", Arg.Set_float stop, "\tStop time (default=10.0)"; "-debug", Arg.Set Common.Debug.debug, "\tShow debug information"; + "-sundials", Arg.Set sundials, "\tUse sundials cvode"; ] -let arg s = - Format.eprintf "Unexpected argument: %s\n" s; exit 1 +let anon = ref (fun s -> Format.eprintf "Unexpected argument: %s\n" s; exit 1) let usage = "" +let register_args l = opts := !opts @ l +let register_anon f = anon := f +let parse_args () = Arg.parse (Arg.align !opts) !anon usage + let go (input : time -> 'a) (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) (output : (time * 'b) -> unit) -= Arg.parse options arg usage; +: unit += parse_args (); let input = { h=(!stop); c=Discontinuous; u=input } in let output o = List.iter output @@ Hsim.Utils.sample_tracked o !sample in - let model = Lift.lift model in - let open Hsim in - let solver = Solver.solver_c Solvers.StatefulRK45.InPlace.csolve - Solvers.StatefulZ.InPlace.zsolve in - let open Sim.Sim(State.InPlaceSimState) in - let sim = Hsim.Utils.(compose (run model (d_of_dc solver)) track) in - run_on sim input output + let solver = Solve.(if !sundials then Sundials else RK45) in + Hsim.Utils.run_on (Solve.build_sim solver model) input output + +let go_discrete + (input : unit -> 'a) + (Ztypes.Node { alloc; step; reset } : ('a, 'b) Ztypes.node) + (output : 'b -> unit) +: unit += parse_args (); + let mem = alloc () in + reset mem; + while true do + input () |> step mem |> output + done; () + +let go_2024 + (input : time -> 'a) + (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) + (output : (time * 'b) -> unit) +: unit += parse_args (); + let input = { h=(!stop); c=Discontinuous; u=input } in + let output o = List.iter output @@ Hsim.Utils.sample_tracked o !sample in + let solver = Solve.(if !sundials then Sundials else RK45) in + Hsim.Utils.run_on (Solve.build_sim_2024 solver model) input output diff --git a/src/lib/std/runtime.mli b/src/lib/std/runtime.mli new file mode 100644 index 0000000..e43ec52 --- /dev/null +++ b/src/lib/std/runtime.mli @@ -0,0 +1,24 @@ + +open Hsim.Types + +val register_args : (string * Arg.spec * string) list -> unit +val register_anon : (string -> unit) -> unit +val parse_args : unit -> unit + +val go : + (time -> 'a) -> + (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> + ((time * 'b) -> unit) -> + unit + +val go_2024 : + (time -> 'a) -> + (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> + ((time * 'b) -> unit) -> + unit + +val go_discrete : + (unit -> 'a) -> + ('a, 'b) Ztypes.node -> + ('b -> unit) -> + unit diff --git a/src/lib/std/solve.ml b/src/lib/std/solve.ml new file mode 100644 index 0000000..1072a6d --- /dev/null +++ b/src/lib/std/solve.ml @@ -0,0 +1,228 @@ + +open Hsim +open Types + +type nonrec 'a value = 'a value +type nonrec 'a signal = 'a signal +type nonrec 'a signal_t = 'a signal_t + +type time = float + +type solver = RK45 | Sundials + +(** Get a value's horizon [h] (reminder: a value is defined on [[0,h]]). *) +let horizon { h; _ } = h + +(** Create a value from a horizon and function. *) +let make (h, u) = { h; u; c=Discontinuous } + +(** Apply a value at a time t. *) +let apply ({ u; h; _ }, t) = + if t > h then raise (Invalid_argument (Format.sprintf + "Requested time t=%.10e is greater than the horizon h=%.10e" t h)); + u t + +let build_sim + (solver : solver) + (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) +: (unit * + ((Ztypes.cvec, Ztypes.dvec) Solver.ivp * + (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode += let model = Lift.lift model in + let solver = Hsim.Solver.solver + (match solver with + | RK45 -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve () + | Sundials -> Solvers.StatefulSundials.InPlace.csolve ()) + (d_of_dc @@ Solvers.StatefulZ.InPlace.zsolve ()) in + let open Hsim.Sim.Sim(Hsim.State.InPlaceSimState) in + let DNode s = Hsim.Utils.(compose (run model solver) track) in + DNode { s with reset=fun p -> s.reset (p, ())} + +let build_sim_2024 + (solver : solver) + (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) +: (unit * + ((Ztypes.cvec, Ztypes.dvec) Solver.ivp * + (Ztypes.cvec, Ztypes.zoutvec) Solver.zc), 'a signal, 'b signal_t) dnode += let model = Lift.lift_2024 model in + let solver = Hsim.Solver.solver + (match solver with + | RK45 -> d_of_dc @@ Solvers.StatefulRK45.InPlace.csolve () + | Sundials -> Solvers.StatefulSundials.InPlace.csolve ()) + (d_of_dc @@ Solvers.StatefulZ.InPlace.zsolve ()) in + let open Hsim.Sim.Sim(Hsim.State.InPlaceSimState) in + let DNode s = Hsim.Utils.(compose (run model solver) track) in + DNode { s with reset=fun p -> s.reset (p, ())} + +(** Lift a hybrid node into a discrete node on streams of functions. *) +let solve + (solver : solver) + (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) +: ('a signal, 'b signal_t) Ztypes.node += let DNode sim = build_sim solver model in + let alloc () = ref sim.state in + let step s a = let b, s' = sim.step !s a in s := s'; b in + let reset _ = () in + Ztypes.Node { alloc; step; reset } + +let solve_2024 + (solver : solver) + (model : Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) +: ('a signal, 'b signal_t) Ztypes.node += let DNode sim = build_sim_2024 solver model in + let alloc () = ref sim.state in + let step s a = let b, s' = sim.step !s a in s := s'; b in + let reset _ = () in + Ztypes.Node { alloc; step; reset } + +let solve_ode45 m = solve RK45 m +let solve_ode45_2024 m = solve_2024 RK45 m +let solve_sundials m = solve Sundials m +let solve_sundials_2024 m = solve_2024 Sundials m + +(** Utility function for [synchr]. + + During synchronization, step the simulation that is lagging behind ([m]) and + join it with the stored value for the other ([n]). + Takes as arguments: + - The step method for [m]; + - The input; + - The last stop times for [m] and [n]; + - The state of [m]; + - The stored value for [n]. + + Returns: + - The common output value up to the common reached date; + - The new reached date of [m]; + - The stored value for [m]; + - The stored value for [n]. *) +let synchr_neq + (m_step : 'ms -> 'a signal -> 'b signal_t) + (input : 'a signal) + (m_stop : time) (n_stop : time) (m_state : 'ms) (n_value : 'c value) +: ('b * 'c) signal_t * time * 'b signal * 'c signal += match m_step m_state input with + | None -> None, m_stop, None, Some n_value + | Some (m_value, m_start) -> + let m_stop = m_start +. m_value.h in + let m_value, n_value, m_rest, n_rest = + (* Three possible scenarios: *) + if m_stop < n_stop then begin + (* [m] is still behind [n]: cut off [n_value] at [m_stop'] *) + let n_value, n_rest = Utils.cutoff n_value m_value.h in + m_value, n_value, None, Some n_rest + end else if n_stop < m_stop then begin + (* [m] overtakes [n]: cut off [m_value] at [n_stop] *) + let m_value, m_rest = Utils.cutoff m_value (n_stop -. m_start) in + m_value, n_value, Some m_rest, None + end else + (* [m] reaches [n] exactly: *) + m_value, n_value, None, None in + let mn_value = Utils.join m_value n_value in + Some (mn_value, m_start), m_stop, m_rest, n_rest + +(** Utility function for [synchr]. + + During synchronization, step both simulations at the same time. + Takes as arguments: + - The step functions for both simulations; + - The input; + - The states of both simulations; + - The last stop times of both simulations. + + Returns: + - The common output value up to the common reached date; + - The new stop times for both simulations; + - The new stored values for both simulations. *) +let synchr_eq + (m_step : 'ms -> 'a signal -> 'b signal_t) + (n_step : 'ns -> 'a signal -> 'c signal_t) + (input : 'a signal) (m_state : 'ms) (n_state : 'ns) + (m_stop : time) (n_stop : time) +: ('b * 'c) signal_t * time * time * 'b signal * 'c signal += match m_step m_state input, n_step n_state input with + | Some (m_value, m_start), Some (n_value, n_start) -> + let m_stop, n_stop = m_start +. m_value.h, n_start +. n_value.h in + let m_value, n_value, m_rest, n_rest = + if m_stop < n_stop then + let n_value, n_rest = Utils.cutoff n_value m_value.h in + m_value, n_value, None, Some n_rest + else if m_stop > n_stop then + let m_value, m_rest = Utils.cutoff m_value n_value.h in + m_value, n_value, Some m_rest, None + else m_value, n_value, None, None in + let mn_value = Utils.join m_value n_value in + Some (mn_value, min m_stop n_stop), m_stop, n_stop, m_rest, n_rest + | None, None -> None, m_stop, n_stop, None, None + | _ -> assert false + +(** Synchronize two simulations as much as possible. *) +let synchr + (m : ('a signal, 'b signal_t) Ztypes.node) + (n : ('a signal, 'c signal_t) Ztypes.node) +: ('a signal, ('b * 'c) signal_t) Ztypes.node += let Ztypes.Node { alloc=m_alloc; step=m_step; reset=m_reset } = m in + let Ztypes.Node { alloc=n_alloc; step=n_step; reset=n_reset } = n in + let alloc () = + ref ((0.0, None, m_alloc ()), (0.0, None, n_alloc ())) in + let step state input = + let (m_stop, m_value, m_state), (n_stop, n_value, n_state) = !state in + let m_stop, m_rest, m_state, n_stop, n_rest, n_state, output = + if m_stop < n_stop then + let n_value = Option.get n_value in + let output, m_stop, m_rest, n_rest = + synchr_neq m_step input m_stop n_stop m_state n_value in + m_stop, m_rest, m_state, n_stop, n_rest, n_state, output + else if m_stop > n_stop then + let m_value = Option.get m_value in + let output, n_stop, n_rest, m_rest = + synchr_neq n_step input n_stop m_stop n_state m_value in + let output = Option.map (fun (u, t) -> Utils.swap u, t) output in + m_stop, m_rest, m_state, n_stop, n_rest, n_state, output + else + let output, m_stop, n_stop, m_rest, n_rest = + synchr_eq m_step n_step input m_state n_state m_stop n_stop in + m_stop, m_rest, m_state, n_stop, n_rest, n_state, output in + state := (m_stop, m_rest, m_state), (n_stop, n_rest, n_state); + output in + let reset ({ contents=((_, _, ms), (_, _, ns)) } as s) = + n_reset ns; m_reset ms; s := (0.0, None, ms), (0.0, None, ns) in + Ztypes.Node { alloc; step; reset } + +(** Sample a value [n] times and iterate [f] on the samples. *) +let iter n f = + let Ztypes.Node { alloc; step; reset } = f in + let step s = + Option.iter @@ fun (v, _) -> + List.iter (fun (_, v) -> step s v) @@ Utils.sample v n in + Ztypes.Node { alloc; step; reset } + +(** Sample a value [n] times and iterate [f] on the dated samples. *) +let iter_t n f = + let Ztypes.Node { alloc; step; reset } = f in + let step s = + Option.iter @@ fun (v, h) -> + List.iter (fun (t, v) -> step s (t +. h, v)) @@ Utils.sample v n in + Ztypes.Node { alloc; step; reset } + +(** Sample a value [n] times and assert [f] on the samples. *) +let check + (n : int) + (Ztypes.Node { alloc; step; reset } : ('a, bool) Ztypes.node) +: ('a signal_t, unit) Ztypes.node += let step s (now, v) = + try assert (step s v) + with Assert_failure _ -> + (Format.eprintf "Assertion failed at time %.10e\n" now; exit 1) in + iter_t n (Ztypes.Node { alloc; reset; step }) + +(** Sample a value [n] times and assert [f] on the dated samples. *) +let check_t + (n : int) + (Ztypes.Node { alloc; step; reset } : (time * 'a, bool) Ztypes.node) +: ('a signal_t, unit) Ztypes.node += let step s (now, v) = + try assert (step s (now, v)) + with Assert_failure _ -> + (Format.eprintf "Assertion failed at time %.10e\n" now; exit 1) in + iter_t n (Ztypes.Node { alloc; reset; step }) diff --git a/src/lib/std/solve.mli b/src/lib/std/solve.mli new file mode 100644 index 0000000..f936e45 --- /dev/null +++ b/src/lib/std/solve.mli @@ -0,0 +1,61 @@ + +type time = float +type 'a value = 'a Hsim.Types.value +type 'a signal = 'a value option +type 'a signal_t = ('a value * time) option + +type solver = RK45 | Sundials + +val horizon : 'a value -> time +val make : time * (time -> 'a) -> 'a value +val apply : 'a value * time -> 'a + +val build_sim : + solver -> + (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> + (unit * + ((Ztypes.cvec, Ztypes.dvec) Hsim.Solver.ivp * + (Ztypes.cvec, Ztypes.zoutvec) Hsim.Solver.zc), + 'a signal, 'b signal_t) Hsim.Types.dnode + +val build_sim_2024 : + solver -> + (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> + (unit * + ((Ztypes.cvec, Ztypes.dvec) Hsim.Solver.ivp * + (Ztypes.cvec, Ztypes.zoutvec) Hsim.Solver.zc), + 'a signal, 'b signal_t) Hsim.Types.dnode + +val solve : + solver -> + (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node + +val solve_2024 : + solver -> + (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node + +val solve_ode45 : + (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node +val solve_ode45_2024 : + (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node +val solve_sundials : + (Ztypes.cstate -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node +val solve_sundials_2024 : + (Ztypes.cstate_new -> (time * 'a, 'b) Ztypes.node) -> + ('a signal, 'b signal_t) Ztypes.node + +val synchr : + ('a signal, 'b signal_t) Ztypes.node -> + ('a signal, 'c signal_t) Ztypes.node -> + ('a signal, ('b * 'c) signal_t) Ztypes.node + +val iter : int -> ('a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node +val iter_t : int -> (time * 'a, unit) Ztypes.node -> ('a signal_t, unit) Ztypes.node + +val check : int -> ('a, bool) Ztypes.node -> ('a signal_t, unit) Ztypes.node +val check_t : int -> (time * 'a, bool) Ztypes.node -> ('a signal_t, unit) Ztypes.node diff --git a/src/lib/std/solve.zli b/src/lib/std/solve.zli new file mode 100644 index 0000000..7670f25 --- /dev/null +++ b/src/lib/std/solve.zli @@ -0,0 +1,23 @@ + +type time = float +type 'a value +type 'a signal = 'a value option +type 'a signal_t = ('a value * time) option + +val horizon : 'a value -> time +val make : time -> (time -> 'a) -> 'a value +val apply : 'a value -> time -> 'a + +val solve_ode45 : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t +val solve_sundials : ('a -C-> 'b) -S-> 'a signal -D-> 'b signal_t + +val synchr : + ('a signal -D-> 'b signal_t) -S-> + ('a signal -D-> 'c signal_t) -S-> + 'a signal -D-> ('b * 'c) signal_t + +val iter : int -S-> ('a -D-> unit) -S-> 'a signal_t -D-> unit +val iter_t : int -S-> (time * 'a -D-> unit) -S-> 'a signal_t -D-> unit + +val check : int -S-> ('a -D-> bool) -S-> 'a signal_t -D-> unit +val check_t : int -S-> (time * 'a -D-> bool) -S-> 'a signal_t -D-> unit diff --git a/src/lib/std/ztypes.ml b/src/lib/std/ztypes.ml index 82ee092..6ce393f 100644 --- a/src/lib/std/ztypes.ml +++ b/src/lib/std/ztypes.ml @@ -43,10 +43,10 @@ type ('a, 'b) cnode = open Bigarray -type time = float -type cvec = (float, float64_elt, c_layout) Array1.t -type dvec = (float, float64_elt, c_layout) Array1.t -type zinvec = (int32, int32_elt, c_layout) Array1.t +type time = float +type cvec = (float, float64_elt, c_layout) Array1.t +type dvec = (float, float64_elt, c_layout) Array1.t +type zinvec = (int32, int32_elt, c_layout) Array1.t type zoutvec = (float, float64_elt, c_layout) Array1.t (* The interface with the ODE solver *) @@ -67,6 +67,23 @@ type cstate = mutable major : bool; (* integration iff [major = false] *) } +(* The interface with the ODE solver (new Zélus version). *) +type cstate_new = + { mutable dvec : dvec; (* Derivative vector. *) + mutable cvec : cvec; (* Position vector. *) + mutable zinvec : zinvec; (* Zero-crossing result vector. *) + mutable zoutvec : zoutvec; (* Zero-crossing value vector. *) + mutable cindex : int; (* Position in position vector. *) + mutable zindex : int; (* Position in zero-crossing vector. *) + mutable cend : int; (* End of position vector. *) + mutable zend : int; (* End of zero-crossing vector. *) + mutable cmax : int; (* Maximum size of position vector. *) + mutable zmax : int; (* Maximum size of zero-crossing vector. *) + mutable horizon : float; (* Next horizon. *) + mutable major : bool; (* Step mode: major <-> discrete. *) + mutable time : float; (* Simulation time. *) + } + (* A hybrid node is a node that is parameterised by a continuous state *) (* all instances points to this global parameter and read/write on it *) type ('a, 'b) hnode = cstate -> (time * 'a, 'b) node