hms/hsim
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge pull request 'no-start is the new main' (#1) from no-start into main

Browse source 
Reviewed-on: https://codeberg.org/17maiga/hsim/pulls/1
This commit is contained in:
17maiga 2025-06-23 15:53:22 +02:00
commit d16e61cc60
29 changed files with 1423 additions and 334 deletions
Download patch file Download diff file Expand all files Collapse all files

BIN
doc/img/ball.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 28 KiB

197
doc/notes.typ
View file

@ -129,13 +129,11 @@ cases:
#columns(4, [#align(center, { #columns(4, [#align(center, {
import cetz: * import cetz: *
import cetz-plot.plot: * import cetz-plot.plot: *
let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
canvas({ canvas({
plot( plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none, size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
x-ticks: ((1,none),), x-grid: "both", y-tick-step: none, y-label: none, x-ticks: ((1,none),), x-grid: "both", y-tick-step: none, y-label: none,
y-min: 1, y-ticks: ((4, none),), y-grid: "both", y-min: 1, y-ticks: ((4, none),), y-grid: "both", {
style: p, {
add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline") add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline")
add(((1,4),), mark: "o", mark-size: .03) add(((1,4),), mark: "o", mark-size: .03)
add(((1,5),(1.3,4.7),(1.7,3.2),(2,3)), line: "spline") add(((1,5),(1.3,4.7),(1.7,3.2),(2,3)), line: "spline")
@ -180,7 +178,6 @@ interval `[now, now + h]`, and end in one of three possible cases:
#columns(3, [#align(center, { #columns(3, [#align(center, {
import cetz: * import cetz: *
import cetz-plot.plot: * import cetz-plot.plot: *
let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true)
canvas({ canvas({
plot( plot(
size: (2,2), axis-style: "left", x-tick-step: none, x-label: none, size: (2,2), axis-style: "left", x-tick-step: none, x-label: none,
@ -319,6 +316,71 @@ advantage of not needing/having optional values in the input and output streams,
and the calling system does not need to guess how many steps the simulation will and the calling system does not need to guess how many steps the simulation will
need to reach the requested horizon. need to reach the requested horizon.
=== Composing simulations
The composition of two simulations must maintain the following invariant: the
output is only ever absent if the solvers have finished integrating up to the
requested horizon. This is because we need to keep providing absent values
(`None`) to the first simulation until it reaches the requested horizon,
otherwise we reset beforehand, and the only information we have on whether the
simulation has reached its horizon is the nature of its output. Naïve
composition breaks this:
#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
import cetz.draw: *
rect((0, 0), (1, 1))
rect((2, 0), (3, 1))
content((0.5, 0.5), $M$)
content((2.5, 0.5), $N$)
line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
line((1, 0.5), (2, 0.5), mark: (end: "straight"))
line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
}),
$ M_i & : & a_1 & | & bot & |
\ M_o & : & b_1 & | & b_2 & |
\ N_i & : & b_1 & | & b_2 & |
\ N_o & : & c_1 & | & #text(fill: red, $c_2$) & |
$)
Instead, we need to make sure we only provide $N$ with a new value when
needed. This is done through the following composition:
#grid(columns: (1fr, 1fr), align: center + horizon, cetz.canvas({
import cetz.draw: *
rect((0, 0), (1, 1))
rect((2, 0), (3, 1))
content((0.5, 0.5), $M$)
content((2.5, 0.5), $N$)
line((-0.5, 0.5), (0, 0.5), mark: (end: "straight"))
line((1, 0.5), (2, 0.5), mark: (end: "straight"))
line((3, 0.5), (3.5, 0.5), mark: (end: "straight"))
}),
$ I &: & a_1 & | & bot & | & bot & & & | & bot & | & bot & & & |
\ M_i &: & a_1 & | & & | & & | & bot & | & & | & & | & bot & |
\ M_o &: & b_1 & | & & | & & | & b_2 & | & & | & & | & bot & |
\ N_i &: & b_1 & | & bot & | & bot & | & b_2 & | & bot & | & bot & | & bot & |
\ N_o &: & c_1 & | & c_2 & | & bot & | & c_3 & | & c_4 & | & bot & | & bot & |
\ O &: & c_1 & | & c_2 & | & & & c_3 & | & c_4 & | & & & bot & |
$)
```ocaml
let step (ms, ns) = function
| Some i ->
let v, ms = m.step ms (Some i) in
let o, ns = n.step ns v in
o, (ms, ns)
| None ->
let o, ns = n.step ns None in
match o with
| Some o -> Some o, (ms, ns)
| None ->
let v, ms = m.step ms None in
match v with
| None -> None, (ms, ns)
| Some v ->
let o, ns = n.step ns (Some v) in
o, (ms, ns)
```
=== Resets === Resets
==== Solver reset ==== Solver reset
@ -411,6 +473,67 @@ Two possible options for the simulation reset:
makes this impossible. We thus need reset parameters for both the model and makes this impossible. We thus need reset parameters for both the model and
solver. solver.
=== Assertions
Assertions are themselves hybrid systems, with their own internal state, `fder`,
`fzer`, `fout` and `reset` functions, etc.
A hybrid model with assertions is then a recursive datatype
```ocaml
type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
HNodeA : {
body : ('s, 'p,' a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
assertions : ('p, 's, float, 'y, 'yder, 'zin, 'zout) hnode_a list;
}
```
Assertions in sub-models need to be "lifted" to the parent model in order to be
visible for the simulation. Since the inner state of sub-models is contained in
the state of the parent model, these assertions can be composed with projectors
on the inner state of the parent model in order to provide them with the correct
state.
What does it mean for an assertion to run in continuous time ?
There are two possible approaches:
- Sampling: take "snapshots" of the state at various timestamps and check the
assertions at these snapshots.
- Problems: expensive memory usage, and requires state copies for the model.
Can also impact stop times for the parent model, which defeats the point of
transparent assertions.
- Real-time: consider the solution provided by the solver as a continuous
function and use it as input to monitor for zero-crossings
- Problems: limits the expressions we can use in continuous assertions to
continuous functions; boolean expressions are inherently discontinuous. We
can model inequalities using zero-crossing constructs like `up`
(for instance, `a <= b` $-->$ `up(a -. b)`).
Building a continuous function of the entire state is possible. Since the entire
state can be considered as the sum of a discrete and a continuous part, and that
the discrete part is considered piecewise-constant, we can build a function
`time -> 's` using the model's `cset` function (which updates the state's
continuous part) as follows:
#grid(columns: (1fr, 1fr), align: center + horizon, ```ocaml
let f_st (ms : 's) (f : time -> 'y) =
fun t -> cset ms (f t)
```, cetz.canvas({
import cetz-plot.plot: *
plot(
axis-style: "left", x-label: none, x-tick-step: none,
y-label: none, y-tick-step: none, y-min: 0.1, {
add(((0, 1), (1, 1)), style: (stroke: blue), label: $D$)
add(((1, 2), (2.5, 2)), style: (stroke: blue))
add(domain: (0, 1), t => - calc.sin(t + 1) + 1.5, style: (stroke: red))
add(domain: (1, 2.5), t => calc.sin(t) + 2, style: (stroke: red), label: $C$)
})
}))
If `cset` is pure (that is, it returns a copy of the state, rather than modify
it in-place), this function poses no problem. If `cset` is not pure and modifies
the state in-place, this is still useable as long as `f_state` is never used in
parallel: since `cset` only modifies the continuous part of the state, we can
always recover the original state by calling `f_state` with the appropriate
horizon. However, the model's continuous state needs to be re-updated to the
reached horizon after all assertions/submodels have finished running.
== Mathematical model == Mathematical model
#link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the #link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the
@ -497,6 +620,70 @@ simulation loop as follows:
models has much in common with code generation for synchronous languages. models has much in common with code generation for synchronous languages.
]) ])
In [Branicky, 1995b], dynamical systems are defined as follows:
#block(fill: rgb(230, 230, 230), stroke: black, inset: 5pt, [
A continuous (resp. discrete) dynamical system defined on the topological
space $X$ over the semigroup $bb(R)^+$ (resp. $bb(Z)^+$) is a function
$f : X times bb(R)^+ -> X$ (resp. $f : X times bb(Z)^+ -> X$) with the
following three properties:
1. initial condition: $f(p, 0) = p$,
2. continuity on both arguments,
3. semigroup property: $f(f(p, t_1), t_2) = f(p, t_1 + t_2)$
for any point $p in X$ and any $t_1, t_2 in bb(R)^+$ (resp $bb(Z)^+$).
Technically, such functions are referred to as semidynamical systems, with the
term dynamical system reserved for those which the semigroups $bb(R)^+$,
$bb(Z)^+$ above may be replaced by the groups $bb(R)$, $bb(Z)$. However, the
more "popular" notion of dynamical system in math and engineering - and the
one used here - requires only the semigroup property. Thus, the term
_reversible_ dynamical system is used when it is necessary to distinguish the
group from semigroup case.
The shorthand $[X, S, f]$ denotes the dynamical system $f$ defined on $X$ over
the semigroup $S$; $X$ is referred to as its _state space_ and points in $X$
are called _states_. If a dynamical system is defined on a subset of $X$, we
say it is a dynamical system _in_ $X$.
For every fixed value of the parameter $s$, the function $f(dot, s)$ defines
a mapping of the space $X$ into itself. Given $[X, bb(Z)^+, f]$, $f(dot, 1)$
is its _transition function_. Thus if $[X, bb(Z)^+, f]$ is reversible, its
transition function is invertible, with inverse given by $f(dot, -1)$.
The set $f(p, S) = { f(p, i) | i in S }$ is called the _orbit_ or _trajectory_
of the point $p$. A _fixed point_ of $[X, S, f]$ is a point $p$ such that
$f(p, s) = p$ for all $s in S$. A set $A subset X$ is _invariant w.r.t._ $f$,
or simply _invariant_, if $f(A, s) subset A$ for all $s in S$.
The notions of equivalence and homomorphism are crucial. Two dynamical systems
$[X, S, f]$, $[Y, S, g]$ are said to be _isomorphic_ (also
_topologically equivalent_ or simply _equivalent_) if there exists a
homeomorphism $psi : X -> Y$ such that
$ psi(f(p, s)) = g(psi(p), s) $
for all $p in X$ and $s in S$. If the mapping $psi$ is only continuous, then
$[X, S, f]$ is said to be _homomorphic_ to $[Y, S, g]$. Homomorphisms preserve
trajectories, fixed points, and invariant sets.
In this paper, the continuous dynamical systems dealt with are defined by the
solutions of ordinary differential equations (ODEs):
$ x'(t) = f(x(t)) $
where $x(t) in X subset bb(R)^n$. The function $f : X -> bb(R)^n$ is called a
_vector field_ on $bb(R)^n$. The resulting dynamical system is then given by
$phi(x_0, t) = x(t)$ where $x(dot)$ is the solution to the above equation
starting at $x_0$ at $t = 0$. (We assume existence and uniqueness of
solutions; see [Hirsch, Smalle] for conditions). A system of ODEs is called
_autonomous_ or _time-invariant_ if its vector field does not depend
explicitly on time. Throughout, the shorthand _continuous_ (resp. _Lipschitz_)
_ODEs_ denotes ODEs with continuous (resp. Lipschitz) vector fields.
An ODE _with input and outputs_ is given by
$ x'(t) & = f(x(t), u(t)), \ y(t) & = h(x(t)) $
where $x(t) in X subset bb(R)^n$, $u(t) in U subset bb(R)^m$,
$y in Y subset bb(R)^p$, $f : bb(R)^n times bb(R)^m -> bb(R)^n$, and
$h : bb(R)^n -> bb(R)^p$. The functions $u(dot)$ and $y(dot)$ are the _inputs_
and _outputs_, respectively.
])
== Miscellaneous == Miscellaneous
=== Signals === Signals
@ -507,7 +694,7 @@ $f : S u p e r d e n s e(bb(V))$, $f(t, n) = f(t + n partial)$ in the
non-standard semantics. non-standard semantics.
Our representation instead uses Our representation instead uses
$S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we $S i g n a l(bb(V)) = S t r e a m((h : bb(R)) * ([0, h] -> bb(V)))$. Can we
convert between the two as we want? convert between the two as we want?
#breakl(```ocaml #breakl(```ocaml

16
doc/pres.typ
View file

@ -145,7 +145,7 @@ Solution: simulate assertions with their own solver.
Nodes take an additional reset parameter `'p` for their reset function: Nodes take an additional reset parameter `'p` for their reset function:
#align(top)[ #align(top)[
#grid(columns: (3fr, 2fr), align: (left, left), #grid(columns: (3fr, 2fr), align: (left, left),
```ocaml ```ocaml
type ('p, 'a, 'b) dnode = type ('p, 'a, 'b) dnode =
DNode : { DNode : {
@ -200,7 +200,7 @@ Hybrid nodes are a bit more complex:
#slide(repeat: 3, self => [ #slide(repeat: 3, self => [
#let (uncover, only) = utils.methods(self) #let (uncover, only) = utils.methods(self)
ODE solvers are discrete nodes producing streams of functions defined on ODE solvers are discrete nodes producing streams of functions defined on
successive intervals: successive intervals:
$(h: bb(R)_+) -->^D (h': bb(R)_+) times (u: [0,h'] -> bb(V))$. $(h: bb(R)_+) -->^D (h': bb(R)_+) times (u: [0,h'] -> bb(V))$.
@ -249,7 +249,7 @@ Hybrid nodes are a bit more complex:
#slide(repeat: 3, self => [ #slide(repeat: 3, self => [
#let (uncover, only) = utils.methods(self) #let (uncover, only) = utils.methods(self)
Zero-crossing solvers are discrete nodes too, looking for zero-crossings on Zero-crossing solvers are discrete nodes too, looking for zero-crossings on
functions: functions:
$(h: bb(R)_+) times (u: [0, h] -> bb(V)) -->^D $(h: bb(R)_+) times (u: [0, h] -> bb(V)) -->^D
(h': bb(R)_+) times (z: bb(Z)?)$. (h': bb(R)_+) times (z: bb(Z)?)$.
@ -456,7 +456,7 @@ When receiving a new value, store it
box((10, 1.25), (3, 1.5), double: true, content: `state`) box((10, 1.25), (3, 1.5), double: true, content: `state`)
// (h, u) -> state // (h, u) -> state
content((-2.25, 3), $(h, u)$) content((-2.25, 3), $(h, u)$)
arrow((-1, 3), (r, 1.5), (d, 2.5), (r, 5), (u, 1.5), (r, 4.5)) arrow((-1, 3), (r, 1.5), (d, 2.5), (r, 5), (u, 1.5), (r, 4.5))
// sim -> bot // sim -> bot
@ -583,7 +583,7 @@ When receiving a new value, store it
// solver.step // solver.step
box((6.5, 1), (3, 1.5), content: `step`) box((6.5, 1), (3, 1.5), content: `step`)
// solver.state // solver.state
box((6.5, 2.75), (3, 1.5), double: true, content: `state`) box((6.5, 2.75), (3, 1.5), double: true, content: `state`)
@ -757,7 +757,7 @@ If no zero-crossing is found, there is no discontinuity.
content((-3, 1.5), $"Signal"(alpha)$) content((-3, 1.5), $"Signal"(alpha)$)
content((5.5, 2), $"Signal"(beta)$) content((5.5, 2), $"Signal"(beta)$)
content((14, 1.5), $"Signal"(gamma)$) content((14, 1.5), $"Signal"(gamma)$)
content((-4, -1), $alpha:$) content((-4, -1), $alpha:$)
content((-4, -2.5), $beta:$) content((-4, -2.5), $beta:$)
content((-4, -4), $gamma:$) content((-4, -4), $gamma:$)
@ -795,7 +795,7 @@ If no zero-crossing is found, there is no discontinuity.
The output signal stops at least as often as the input signal does. The output signal stops at least as often as the input signal does.
#linebreak() #linebreak()
This calls for a "lazy" composition: we request rather than provide data. This calls for a "lazy" composition: we request rather than provide data.
#align(center, cetz-canvas({ #align(center, cetz-canvas({
import cetz.draw: * import cetz.draw: *
box((0, 0), (3, 3), content: $M_1$) box((0, 0), (3, 3), content: $M_1$)
@ -809,7 +809,7 @@ If no zero-crossing is found, there is no discontinuity.
content((-3, 1.5), $"Signal"(alpha)$) content((-3, 1.5), $"Signal"(alpha)$)
content((5.5, 2), $"Signal"(beta)$) content((5.5, 2), $"Signal"(beta)$)
content((14, 1.5), $"Signal"(gamma)$) content((14, 1.5), $"Signal"(gamma)$)
content((-4, -1), $alpha:$) content((-4, -1), $alpha:$)
content((-4, -2.5), $beta:$) content((-4, -2.5), $beta:$)
content((-4, -4), $gamma:$) content((-4, -4), $gamma:$)

192
doc/rep.typ Normal file
View file

@ -0,0 +1,192 @@
#import "@preview/cetz:0.4.0"
#import "@preview/finite:0.4.1"
#set page(paper: "a4")
#set par(justify: true)
#set raw(syntaxes: "zelus.sublime-syntax")
#show raw: set text(font: "CaskaydiaCove NF")
#let lustre = smallcaps[Lustre]
#let esterel = smallcaps[Esterel]
#let simulink = smallcaps[Simulink]
#let sundials = smallcaps[Sundials CVODE]
#let zelus = smallcaps[Zélus]
#let TODO(..what) = {
let msg = if what.pos().len() == 0 { "" } else { ": " + what.pos().join("") }
block(fill: red, width: 100%, inset: 5pt, align(center, raw("TODO" + msg)))
}
#let adot(s) = $accent(#s, dot)$
#let addot(s) = $accent(#s, dot.double)$
#align(center)[
#text(16pt)[
*Internship Report --- MPRI M2 \
Hybrid System Models with Transparent Assertions*
]
Henri Saudubray, supervised by Marc Pouzet, Inria PARKAS
]\
#heading(outlined: false)[General Context]
What is the report about? Where does the problem come from? What is the state of
the art in that area?
Hybrid systems modelers such as #simulink are essential tools in the development
of embedded systems evolving in physical environments. They allow for precise
descriptions of hybrid systems through both continuous-time behaviour defined
through Ordinary Differential Equations (ODEs) and discrete-time reactive
behaviour similar to what is found in the synchronous languages such as #lustre.
The #zelus language aims to reconcile synchronous languages with hybrid systems,
by taking a synchronous language kernel and extending it with continuous-time
constructs similar to what is found in tools like #simulink. Continuous-time
behaviour is computed through the use of an off-the-shelf ODE solver such as
#sundials.
#heading(outlined: false)[Research Problem]
What specific question(s) have you studied? What motivates the question? What
are the applications/consequences? Is it a new problem? Why did you choose this
problem?
The simulation of hybrid system models, as done in #simulink and #zelus, uses a
single ODE solver instance to simulate the entire model. This has various
advantages: it is less computationally intensive, and simplifies the work of the
compiler. Unfortunately, it also raises a difficult problem: sub-systems which
seemingly should not interfere with each other end up affecting each other's
results. This is due to the chosen integration method: an adaptive solver like
#sundials will vary its step length throughout the integration process, and the
addition of ODEs in the overall system can influence these step lengths,
affecting the results obtained for pre-existing ODEs. This is particularly
problematic in the case of runtime assertions, which are typically expected to
be transparent: they should not affect the final result of the computation.
#heading(outlined: false)[Proposed Contributions]
What is your answer to the research problem? Please remain at a high level;
technical details should be given in the main body of the report. Pay special
attention to the description of the _scientific_ approach.
To solve this, we propose a new runtime for the #zelus language, which simulates
assertions with their own solvers in order to maintain the separation between
assertions and the model they operate on.
#TODO()
#heading(outlined: false)[Arguments Supporting Their Validity]
What is the evidence that your solution is a good solution? (Experiments?
Proofs?) Comment on the robustness of your solution: does it rely on
assumptions, and if so, to what extent?
#TODO("Justify")
#heading(outlined: false)[Summary and Future Work]
What did you contribute to the area? What comes next? What is a good _next_ step
or question?
#TODO("Next steps")
#pagebreak()
#outline()
#pagebreak()
= Background
== Hybrid systems
Hybrid systems are dynamical systems whose behaviour evolves through both
continuous and discrete phases. They are used to model software interacting with
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 $
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:
$ 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
the value of the variable at a given time $t$. We can find an approximation of
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.
@lst:ball.zls shows how such a model might be expressed in the concrete syntax
of #zelus.
#figure(placement: top, gap: 2em,
```zelus
let hybrid ball () = y where
rec der y = y' init 50.0
and der y' = -9.81 init 0.0 reset z -> -0.8 *. (last y')
and z = up (-. y)
```,
caption: [The bouncing ball in #zelus]
) <lst:ball.zls>
More formally, a hybrid system can be described as an automaton
== 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.
#figure(finite.automaton(
(D: (D: "cascade", C: "no cascade"),
C: (C: "no zero-crossing", D: "zero-crossing")),
initial: "D", final: (), layout: finite.layout.linear.with(spacing: 3)
), caption: [High-level overview of the runtime], placement: top) <automaton>
= Runtime
To solve this issue, we need to redefine what the runtime of our hybrid system
models is. The runtime is the piece of software that handles the pairing of a
hybrid model with a solver and the different execution phases that need to take
place.
To allow for independent simulation of the assertions, we need to simulate them
with their own ODE solvers. Since assertions can themselves contain ODEs, they
can be viewed as hybrid systems themselves. Assertions can also contain their
own sub-assertions, and so recursively: our runtime needs to handle this
accordingly.
== Hybrid models with assertions
A hybrid model with assertions can be defined using a recursive data type:
```ocaml
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;
} -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
```
Notice that the list of assertions takes as input the inner state of the parent
model: assertions are checked on the model state, and any variable which is
required by the assertion becomes a state variable.
== Solvers as synchronous nodes
== Simulations as synchronous nodes
#TODO("talk about the new runtime")
= Assertions
#TODO("talk about how assertions are done")
#pagebreak()
= Annex
#set heading(level: 2)
#bibliography("sources.bib", full: true)
#set heading(level: 1)

12
exm/ball.ml
View file

@ -29,7 +29,7 @@ let fder y yd =
else begin yd.{0} <- 0.0; yd.{1} <- 0.0; yd.{2} <- 0.0; yd.{3} <- 0.0 end; else begin yd.{0} <- 0.0; yd.{1} <- 0.0; yd.{2} <- 0.0; yd.{3} <- 0.0 end;
yd yd
let fzer y zo = zo.{0} <- -. y.{1}; zo let fzer y zo = zo.{0} <- -. y.{1}; zo
let fout _ _ y = of_array [| y.{1} |] let fout _ _ y = of_array [| y.{1}; y.{0} |]
let jump _ = true let jump _ = true
let horizon _ = max_float let horizon _ = max_float
let cget s = s.lx let cget s = s.lx
@ -38,9 +38,9 @@ let zset s zin = { s with zin }
let step ({ zin; lx; _ } as s) zfalse = let step ({ zin; lx; _ } as s) zfalse =
let lx = if zin.{0} = 1l then let lx = if zin.{0} = 1l then
of_array [| -. 0.8 *. lx.{0}; 0.0; lx.{2}; lx.{3} |] else lx in of_array [| -. 0.8 *. lx.{0}; 0.0; lx.{2}; lx.{3} |] else lx in
of_array [| s.lx.{1} |], { zin=zfalse; lx; i=false } of_array [| s.lx.{1}; s.lx.{0} |], { zin=zfalse; lx; i=false }
let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode = let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec =
let yd = cmake csize in let yd = cmake csize in
let zout = cmake zsize in let zout = cmake zsize in
let zfalse = zmake 1 in let zfalse = zmake 1 in
@ -49,10 +49,10 @@ let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode =
let step s _ = step s zfalse in let step s _ = step s zfalse in
let state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in let state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
let reset _ _ = state in let reset _ _ = state in
HNode { state; fder; fzer; fout; step; reset; horizon; { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset;
jump; cset; cget; zset; csize; zsize } csize; zsize }
let errmsg = "Too many arguments for the model (needed: 0)" let errmsg = "Too many arguments for the model (needed: 0)"
let init = function let init = function
| [] -> bouncing_ball () | [] -> HNode (bouncing_ball ())
| _ -> raise (Invalid_argument errmsg) | _ -> raise (Invalid_argument errmsg)

11
exm/dune
View file

@ -1,3 +1,10 @@
(env
(dev
(flags
(:standard -w -a))))
(library (library
(name examples) (name examples)
(libraries hsim solvers)) (libraries hsim solvers))
(include_subdirs unqualified)

72
exm/sin1x.ml Normal file
View file

@ -0,0 +1,72 @@
(* Example due to JB. Jeannin - zero-crossing detection *)
(* the signal x(t) is expected to have a zero-crossing at [t = 0.3] *)
(* to draw: ./main.exe -s sin_1_x -stop 2 -sample 10000 | feedgnuplot --stream --domain --lines *)
open Hsim.Types
open Solvers.Zls
(* let hybrid sin_1_x() =
let der t = 1.0 init 0.0 in
let t = t -. 0.3 in
let v = sin(1/t) in
let o = t *. (v *. v) in [that is: t *. (sin(1/t)^2]
present up(o) -> 1.0 else 0.0, o *)
let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a
type ('a, 'b) state = { s_x: 'a; zin: 'b }
let csize = 1
let zsize = 1
let fder _ _ yd = yd.{0} <- 1.0
let fzer _ y zout =
let t = y.{0} in
let t = t -. 0.3 in
let v = sin(1.0 /. t) in
let o = t *. (v *. v) in
zout.{0} <- o
let fout s y =
let z = if s.zin.{0} = 1l then 1.0 else 0.0 in
let t = y.{0} in
let t = t -. 0.3 in
let v = sin(1.0 /. t) in
let o = t *. (v *. v) in
of_array [| z; o |]
let step ({ s_x; zin } as s) zfalse =
let z = if zin.{0} = 1l then 1.0 else 0.0 in
let t = s_x.{0} in
let t = t -. 0.3 in
let v = sin(1.0 /. t) in
let o = t *. (v *. v) in
of_array [| z; o |], { s with zin = zfalse }
let reset _ s = s
let cget { s_x; _ } = s_x
let cset s l_x = { s with s_x = l_x }
let zset s zin = { s with zin }
let jump _ = true
let horizon _ = max_float
let sin_1_x () =
let zfalse = zmake 1 in
let yd = cmake 1 in
let zout = cmake 1 in
let fder _ _ y = fder 0.0 y yd; yd in
let fzer _ _ y = fzer 0.0 y zout; zout in
let fout s _ y = fout s y in
let step s _ = step s zfalse in
let state = { s_x = of_array [| 0.0 |]; zin = zfalse } in
{ state; fder; fzer; fout; step; reset; horizon;
cset; cget; zset; csize; zsize; jump }
let errmsg = "Too many arguments for the model (needed: 0)"
let init = function
| [] -> HNode (sin_1_x ())
| _ -> raise (Invalid_argument errmsg)

79
exm/sin1x_der.ml Normal file
View file

@ -0,0 +1,79 @@
(* Example due to JB. Jeannin - zero-crossing detection. *)
(* The same as [sin1x.ml] but with an integrator.
d(t . (sin(1/t))^2))/dt = sin(1/t)(sin(1/t) - (2 . cos(1/t)) / t)
d(sin(f(x)))/dx = cos(f(x)) f'(x) *)
open Hsim.Types
open Solvers.Zls
(* let hybrid sin_1_x() =
let der t0 = 1.0 init 0.0 in
let d = 0.6 in
let t = t0 -. d in
let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in
present up(o) -> 1.0 else 0.0, o *)
let of_array a = Bigarray.Array1.of_array Bigarray.Float64 Bigarray.c_layout a
type ('a, 'b) state = { s_x: 'a; zin: 'b }
let d = 0.66
let y0 = -. d *. let v = (sin(1.0 /. (-. d))) in v *. v
let csize = 2
let zsize = 1
let fder d y yd =
yd.{0} <- 1.0;
let t = y.{0} -. d in
yd.{1} <- sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
let fzer z y zout =
let o = y.{1} in
zout.{0} <- if z then o else 1.0
let fout s y =
let z = if s.zin.{0} = 1l then 1.0 else 0.0 in
let o = y.{1} in
of_array [| z; o |]
let step ({ s_x; zin } as s) zfalse =
let z = if zin.{0} = 1l then 1.0 else -1.0 in
let o = s_x.{1} in
of_array [| z; o |], { s with zin = zfalse }
let reset _ s = s
let cget { s_x; _ } = s_x
let cset s l_x = { s with s_x = l_x }
let zset s zin = { s with zin }
let horizon _ = max_float
let jump _ = true
let sin_1_x z d =
let zfalse = zmake 1 in
let yd = cmake 2 in
let zout = cmake 1 in
let fder _ _ y = fder d y yd; yd in
let fzer _ _ y = fzer z y zout; zout in
let fout s _ y = fout s y in
let step s _ = step s zfalse in
let state = { s_x = of_array [| 0.0; y0 |]; zin = zfalse } in
{ state; fder; fzer; fout; step; horizon;
cset; cget; zset; csize; zsize; reset; jump }
let errmsg_many = "Too many arguments to model (needed: [bool, float])"
let errmsg_few = "Too few arguments to model (needed: [bool, float])"
let errmsg_invalid = "Invalid arguments to model (needed: [bool, float])"
let init = function
| [zer; d] ->
let zer, d = try bool_of_string zer, float_of_string d
with _ -> raise (Invalid_argument errmsg_invalid) in
HNode (sin_1_x zer d)
| [] | [_] -> raise (Invalid_argument errmsg_few)
| _ -> raise (Invalid_argument errmsg_many)

8
exm/sincos.ml
View file

@ -34,11 +34,11 @@ let sinus_cosinus theta0 omega =
HNode { state; fder; fzer; fout; step; reset; horizon; HNode { state; fder; fzer; fout; step; reset; horizon;
jump; cset; cget; zset; csize; zsize } jump; cset; cget; zset; csize; zsize }
let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)" let errmsg_invalid = "Invalid arguments to model (needed: [float, float])"
let errmsg_few = "Too few arguments to model (needed: 2 floats)" let errmsg_few = "Too few arguments to model (needed: [float, float])"
let errmsg_many = "Too many arguments to model (needed: 2 floats)" let errmsg_many = "Too many arguments to model (needed: [float, float])"
let init = function let init = function
| [t0; om] -> | [om; t0] ->
let t0, om = try float_of_string t0, float_of_string om let t0, om = try float_of_string t0, float_of_string om
with Failure _ -> raise (Invalid_argument errmsg_invalid) in with Failure _ -> raise (Invalid_argument errmsg_invalid) in
sinus_cosinus t0 om sinus_cosinus t0 om

22
exm/zelus/ballz/ballz.zls Normal file
View file

@ -0,0 +1,22 @@
let g = 9.81
let y0 = 50.0
let y'0 = 0.0
let hybrid ball (y0, y'0) = (y, y', z) where
rec der y = y' init y0
and der y' = -. g init y'0 reset z -> -0.8 *. (last y')
and z = up(-. y)
let hybrid main () =
let der t = 1.0 init 0.0 in
let rec der p = 1.0 init -0.01 reset s -> -0.01
and s = up(p) in
let (y, y', z) = ball (y0, y'0) in
present z | s -> (
print_float t;
print_string "\t";
print_float y;
print_string "\t";
print_float y';
print_newline ()
); ()

6
exm/zelus/ballz/dune Normal file
View file

@ -0,0 +1,6 @@
(rule
(targets ballz.ml ballz.zci)
(deps
(:zl ballz.zls))
(action
(run zeluc %{zl})))

74
exm/zelus/sin_1_x.ml Normal file
View file

@ -0,0 +1,74 @@
(* The Zelus compiler, version 2.2-dev
(2025-06-16-15:24) *)
open Common
open Ztypes
open Solvers
(* open Zls *)
type ('f , 'e , 'd , 'c , 'b , 'a) _sin_1_x =
{ mutable major_22 : 'f ;
mutable i_29 : 'e ;
mutable x_28 : 'd ;
mutable result_27 : 'c ; mutable o_26 : 'b ; mutable t0_23 : 'a }
let sin_1_x (cstate_30:Ztypes.cstate) =
let sin_1_x_alloc _ =
cstate_30.cmax <- (+) cstate_30.cmax 2 ;
cstate_30.zmax <- (+) cstate_30.zmax 1;
{ major_22 = false ;
i_29 = (false:bool) ;
x_28 = { zin = false; zout = 1. } ;
result_27 = (42.:float) ;
o_26 = { pos = 42.; der = 0. } ; t0_23 = { pos = 42.; der = 0. } } in
let sin_1_x_step self ((_time_21:float) , ()) =
((let (cindex_31:int) = cstate_30.cindex in
let cpos_33 = ref (cindex_31:int) in
let (zindex_32:int) = cstate_30.zindex in
let zpos_34 = ref (zindex_32:int) in
cstate_30.cindex <- (+) cstate_30.cindex 2 ;
cstate_30.zindex <- (+) cstate_30.zindex 1 ;
self.major_22 <- cstate_30.major ;
(if cstate_30.major then
for i_1 = cindex_31 to 1 do Zls.set cstate_30.dvec i_1 0. done
else ((self.o_26.pos <- Zls.get cstate_30.cvec !cpos_33 ;
cpos_33 := (+) !cpos_33 1) ;
(self.t0_23.pos <- Zls.get cstate_30.cvec !cpos_33 ;
cpos_33 := (+) !cpos_33 1))) ;
(let (result_35:(float * float)) =
(if self.i_29 then
self.o_26.pos <- ( *. ) ((~-.) 0.6)
(( ** ) (sin ((/.) 1. ((~-.) 0.6))) 2.))
;
self.i_29 <- false ;
self.t0_23.der <- 1. ;
(let (t_25:float) = (-.) self.t0_23.pos 0.6 in
self.o_26.der <- ( *. ) (sin ((/.) 1. t_25))
((-.) (sin ((/.) 1. t_25))
((/.) (( *. ) 2.
(cos ((/.) 1. t_25)))
t_25)) ;
self.x_28.zout <- self.o_26.pos ;
(begin match self.x_28.zin with
| true -> self.result_27 <- 1.
| _ -> self.result_27 <- 0. end) ;
(self.result_27 , self.o_26.pos)) in
cpos_33 := cindex_31 ;
(if cstate_30.major then
(((Zls.set cstate_30.cvec !cpos_33 self.o_26.pos ;
cpos_33 := (+) !cpos_33 1) ;
(Zls.set cstate_30.cvec !cpos_33 self.t0_23.pos ;
cpos_33 := (+) !cpos_33 1)) ; ((self.x_28.zin <- false)))
else (((self.x_28.zin <- Zls.get_zin cstate_30.zinvec !zpos_34 ;
zpos_34 := (+) !zpos_34 1)) ;
zpos_34 := zindex_32 ;
((Zls.set cstate_30.zoutvec !zpos_34 self.x_28.zout ;
zpos_34 := (+) !zpos_34 1)) ;
((Zls.set cstate_30.dvec !cpos_33 self.o_26.der ;
cpos_33 := (+) !cpos_33 1) ;
(Zls.set cstate_30.dvec !cpos_33 self.t0_23.der ;
cpos_33 := (+) !cpos_33 1)))) ; result_35)):float * float) in
let sin_1_x_reset self =
((self.i_29 <- true ; self.t0_23.pos <- 0.):unit) in
Node { alloc = sin_1_x_alloc; step = sin_1_x_step ; reset = sin_1_x_reset }

BIN
exm/zelus/sin_1_x.zci Normal file
View file

Binary file not shown.

7
exm/zelus/sin_1_x.zls Normal file
View file

@ -0,0 +1,7 @@
let hybrid sin_1_x () =
let der t0 = 1.0 init 0.0 in
let d = 0.6 in
let t = t0 -. d in
let der o = sin(1.0 /. t) *. (sin(1.0 /. t) -. (2.0 *. cos(1.0 /. t)) /. t)
init (-. d *. (sin(1.0 /. (-. d))) ** 2.0) in
(present up(o) -> 1.0 else 0.0), o

6
exm/zelus/sincos/dune Normal file
View file

@ -0,0 +1,6 @@
(rule
(targets sincosz.ml sincosz.zci)
(deps
(:zl sincosz.zls))
(action
(run zeluc %{zl})))

4
exm/zelus/sincos/sincosz.zls Normal file
View file

@ -0,0 +1,4 @@
let hybrid f () = (sin, cos) where
rec der sin = cos init 0.0
and der cos = -. sin init 1.0

4
exm/ztypes.ml Normal file
View file

@ -0,0 +1,4 @@
include Common
include Ztypes
include Solvers

118
src/bin/lift.ml Normal file
View file

@ -0,0 +1,118 @@
open Hsim.Types
open Solvers.Zls
open Common.Ztypes
type ('s, 'a) state =
{ mutable state : 's; mutable input : 'a option; mutable time : time }
let lift f =
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
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; _ } input y =
cstate.major <- false;
cstate.zinvec <- no_roots_in;
cstate.zoutvec <- no_roots_out;
cstate.cvec <- y;
cstate.dvec <- ignore_der;
cstate.cindex <- 0;
cstate.zindex <- 0;
ignore (f_step state (no_time, input));
cstate.dvec in
(* the function that compute the zero-crossings *)
let fzer { state; _ } input y =
cstate.major <- false;
cstate.zinvec <- no_roots_in;
cstate.dvec <- ignore_der;
cstate.zoutvec <- no_roots_out;
cstate.cvec <- y;
cstate.cindex <- 0;
cstate.zindex <- 0;
ignore (f_step state (no_time, input));
cstate.zoutvec in
(* the function which compute the output during integration *)
let fout { state; _ } input y =
cstate.major <- false;
cstate.zoutvec <- no_roots_out;
cstate.dvec <- ignore_der;
cstate.zinvec <- no_roots_in;
cstate.cvec <- y;
cstate.cindex <- 0;
cstate.zindex <- 0;
f_step state (no_time, input) in
(* the function which compute a discrete step *)
let step ({ state; time; _ } as st) input =
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 (time, input) in
o, { st with state; input=Some input } in
let reset _ ({ state; _ } as st) = f_reset state; st in
(* horizon *)
let horizon _ = cstate.horizon 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 }

91
src/bin/main.ml
View file

@ -4,17 +4,26 @@ open Solvers
open Examples open Examples
open Common open Common
open Types open Types
open Lift
let sample = ref 10 let sample = ref 1
let stop = ref 30.0 let stop = ref 10.0
let greedy = ref false let accel = ref false
let inplace = ref false let inplace = ref false
let sundials = ref false let sundials = ref false
let speed = ref false
let steps = ref 1 let steps = ref 1
let model = ref None 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 zelus = ref false
let gt0i v i = v := if i <= 0 then 1 else i 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 let gt0f v f = v := if f <= 0.0 then 1.0 else f
let opt r s = r := Some s
let modelargs = ref [] let modelargs = ref []
let set_model s = let set_model s =
@ -26,46 +35,78 @@ let opts = [
"-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)"; "-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)";
"-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)"; "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
"-debug", Arg.Set Debug.debug, "\tPrint debug information"; "-debug", Arg.Set Debug.debug, "\tPrint debug information";
"-greedy", Arg.Set greedy, "\tUse greedy simulation"; "-accelerate", Arg.Set accel, "\tConcatenate continuous functions";
"-sundials", Arg.Set sundials, "\tUse sundials (not compatible with greedy)"; "-sundials", Arg.Set sundials, "\tUse sundials (doesn't support -accelerate)";
"-inplace", Arg.Set inplace, "\tUse imperative solvers"; "-inplace", Arg.Set inplace, "\tUse imperative solvers";
"-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)"; "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)";
"-speed", Arg.Set speed, "\tLog the step length";
"-minstep", Arg.String (opt minstep), "\tSet minimum solver step length";
"-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";
"-zelus", Arg.Set zelus, "\tUse the output of the Zélus compiler";
] ]
let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:" let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2 let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
let m = let args = List.rev !modelargs
try match !model with let () = ignore lift
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ball" -> Ball.init !modelargs
| Some "vdp" -> Vdp.init !modelargs
| Some "sincos" -> Sincos.init !modelargs
| Some "sqrt" -> Sqrt.init !modelargs
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
let z = StatefulZ.(if !inplace then InPlace.zsolve else Functional.zsolve) let wrap_zelus (HNode m) =
let ret = Bigarray.(Array1.create Float64 c_layout 0) in
let fout s a y = ignore (m.fout s a y); ret in
let step s () = let _, s = m.step s () in ret, s in
HNode { m with fout; step }
let m =
try
if !zelus then
match !model with
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ballz" -> wrap_zelus (lift Ballz.main)
| Some "sincosz" -> wrap_zelus (lift Sincosz.f)
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
else
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 s -> Format.eprintf "Unknown model: %s\n" s; exit 2
with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
let st = if !inplace then (module State.InPlaceSimState : State.SimState) let st = if !inplace then (module State.InPlaceSimState : State.SimState)
else (module State.FunctionalSimState : State.SimState) else (module State.FunctionalSimState : State.SimState)
let output =
if !no_print then Output.ignore
else if !speed then Output.print_h
else Output.print (* Output.ignore *)
let sim = let sim =
if !sundials then if !sundials then
if !greedy then let open StatefulSundials in
(Format.eprintf "Sundials does not support greedy simulation\n"; exit 2) let c = if !inplace then InPlace.csolve else Functional.csolve in
else let open StatefulZ in
let open StatefulSundials in let z = if !inplace then InPlace.zsolve else Functional.zsolve in
let c = if !inplace then InPlace.csolve else Functional.csolve in let s = Solver.solver c (d_of_dc z) in
let s = Solver.solver c (d_of_dc z) in let open Sim.Sim(val st) in
let open Sim.LazySim(val st) in run_until_n m s run_until_n (output !sample (run m s))
else else
let open StatefulRK45 in 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 s = Solver.solver_c c z in let s = Solver.solver_c c z in
if !greedy then let open Sim.GreedySim(val st) in run_until_n m s let open Sim.Sim(val st) in
else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s) let n = if !accel then accelerate m s else run m (d_of_dc s) in
run_until_n (output !sample n)
let () = sim !stop !steps (Output.print !sample) let () = sim !stop !steps ignore

58
src/bin/output.ml
View file

@ -1,8 +1,8 @@
open Hsim.Types open Hsim.Types
open Common open Hsim.Utils
let print_entry t y = let print_entry y t =
let n = Bigarray.Array1.dim y in let n = Bigarray.Array1.dim y in
let rec loop i = let rec loop i =
if i = n then () if i = n then ()
@ -12,17 +12,53 @@ let print_entry t y =
Format.printf "\n"; Format.printf "\n";
flush stdout flush stdout
let print samples { start; length; u } = let print_entry_h y t h =
let step = length /. (float_of_int samples) in let n = Bigarray.Array1.dim y in
let rec loop i =
if i = n then ()
else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
Format.printf "% .10e\t% .10e" t h;
loop 0;
Format.printf "\n";
flush stdout
let print_sample samples ({ h; u; _ }, now) =
let step = h /. (float_of_int samples) in
let rec loop i = let rec loop i =
if i > samples then () if i > samples then ()
else if i = samples then print_entry (start +. length) (u length) else if i = samples then print_entry (u h) (now +. h)
else let t = float_of_int i *. step in else let t = float_of_int i *. step in
(print_entry (start +. t) (u t); loop (i+1)) in (print_entry (u t) (now +. t); loop (i+1)) in
if length <= 0.0 then begin Debug.print "D: "; print_entry start (u 0.0) end if h <= 0.0 then print_entry (u 0.0) now else loop 0
else begin Debug.print "C: "; loop 0 end
let print_limits { start; length; _ } = let print_sample_h samples ({ h; u; _ }, now) =
if length <= 0.0 then Format.printf "D: % .10e\n" start let step = h /. (float_of_int samples) in
else Format.printf "C: % .10e\t% .10e\n" start (start +. length) let rec loop i =
if i > samples then ()
else if i = samples then print_entry_h (u h) (now +. h) h
else let t = float_of_int i *. step in
(print_entry_h (u t) (now +. t) h; loop (i+1)) in
if h <= 0.0 then print_entry_h (u 0.0) now h else loop 0
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
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
DNode { m with reset=fun p -> m.reset (p, ((), ())) }
let ignore _ n =
let state = () in
let step () = function
| None -> None, ()
| Some _ -> Some (), () in
let reset () () = () in
let i = DNode { state; step; reset } in
let DNode n = compose n i in
DNode { n with reset=fun p -> n.reset (p, ()) }

7
src/lib/common/debug.ml
View file

@ -2,3 +2,10 @@
let debug = ref false let debug = ref false
let print s = if !debug then Format.printf "%s\n" s else () let print s = if !debug then Format.printf "%s\n" s else ()
let print_entry y =
let n = Bigarray.Array1.dim y in
let rec loop i =
if i = n then ()
else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
if !debug then (loop 0; Format.printf "\n"; flush stdout)

151
src/lib/common/ztypes.ml Normal file
View file

@ -0,0 +1,151 @@
(**************************************************************************)
(* *)
(* Zelus *)
(* A synchronous language for hybrid systems *)
(* http://zelus.di.ens.fr *)
(* *)
(* Marc Pouzet and Timothy Bourke *)
(* *)
(* Copyright 2012 - 2019. All rights reserved. *)
(* *)
(* This file is distributed under the terms of the CeCILL-C licence *)
(* *)
(* Zelus is developed in the INRIA PARKAS team. *)
(* *)
(**************************************************************************)
(* Type declarations and values that must be linked with *)
(* the generated code *)
type 'a continuous = { mutable pos: 'a; mutable der: 'a }
type ('a, 'b) zerocrossing = { mutable zin: 'a; mutable zout: 'b }
type 'a signal = 'a * bool
type zero = bool
(* a synchronous stream function with type 'a -D-> 'b *)
(* is represented by an OCaml value of type ('a, 'b) node *)
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
(* the same with a method copy *)
type ('a, 'b) cnode =
Cnode:
{ alloc : unit -> 's; (* allocate the state *)
copy : 's -> 's -> unit; (* copy the source into the destination *)
step : 's -> 'a -> 'b; (* compute a step *)
reset : 's -> unit; (* reset/inialize the state *)
} -> ('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 zoutvec = (float, float64_elt, c_layout) Array1.t
(* The interface with the ODE solver *)
type cstate =
{ mutable dvec : dvec; (* the vector of derivatives *)
mutable cvec : cvec; (* the vector of positions *)
mutable zinvec : zinvec; (* the vector of boolean; true when the
solver has detected a zero-crossing *)
mutable zoutvec : zoutvec; (* the corresponding vector that define
zero-crossings *)
mutable cindex : int; (* the position in the vector of positions *)
mutable zindex : int; (* the position in the vector of zero-crossings *)
mutable cend : int; (* the end of the vector of positions *)
mutable zend : int; (* the end of the zero-crossing vector *)
mutable cmax : int; (* the maximum size of the vector of positions *)
mutable zmax : int; (* the maximum number of zero-crossings *)
mutable horizon : float; (* the next horizon *)
mutable major : bool; (* integration iff [major = false] *)
}
(* 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
type 'b hsimu =
Hsim:
{ alloc : unit -> 's;
(* allocate the initial state *)
maxsize : 's -> int * int;
(* returns the max length of the *)
(* cvector and zvector *)
csize : 's -> int;
(* returns the current length of the continuous state vector *)
zsize : 's -> int;
(* returns the current length of the zero-crossing vector *)
step : 's -> cvec -> dvec -> zinvec -> time -> 'b;
(* computes a step *)
derivative : 's -> cvec -> dvec -> zinvec -> zoutvec -> time -> unit;
(* computes the derivative *)
crossings : 's -> cvec -> zinvec -> zoutvec -> time -> unit;
(* computes the zero-crossings *)
reset : 's -> unit;
(* resets the state *)
horizon : 's -> time;
(* gives the next time horizon *)
} -> 'b hsimu
(* a function with type 'a -C-> 'b, when given to a solver, is *)
(* represented by an OCaml value of type ('a, 'b) hsnode *)
type ('a, 'b) hsnode =
Hnode:
{ state : 's;
(* the discrete state *)
zsize : int;
(* the maximum size of the zero-crossing vector *)
csize : int;
(* the maximum size of the continuous state vector (positions) *)
derivative : 's -> 'a -> time -> cvec -> dvec -> unit;
(* computes the derivative *)
crossing : 's -> 'a -> time -> cvec -> zoutvec -> unit;
(* computes the derivative *)
output : 's -> 'a -> cvec -> 'b;
(* computes the zero-crossings *)
setroots : 's -> 'a -> cvec -> zinvec -> unit;
(* returns the zero-crossings *)
majorstep : 's -> time -> cvec -> 'a -> 'b;
(* computes a step *)
reset : 's -> unit;
(* resets the state *)
horizon : 's -> time;
(* gives the next time horizon *)
} -> ('a, 'b) hsnode
(* An idea suggested by Adrien Guatto, 26/04/2021 *)
(* provide a means to the type for input/outputs of nodes *)
(* express them with GADT to ensure type safety *)
(* type ('a, 'b) node =
| Fun : { step : 'a -> 'b;
typ_arg: 'a typ;
typ_return: 'b typ
} -> ('a, 'b) node
| Node :
{ state : 's; step : 's -> 'a -> 'b * 's;
typ_arg: 'a typ;
typ_state : 's typ;
typ_return: 'b typ } -> ('a, 'b) node
and 'a typ =
| Tunit : unit typ
| Tarrow : 'a typ * 'b typ -> ('a * 'b) typ
| Tint : int -> int typ
| Ttuple : 'a typlist -> 'a typ
| Tnode : 'a typ * 'b typ -> ('a,'b) node typ
and 'a typlist =
| Tnil : unit typlist
| Tpair : 'a typ * 'b typlist -> ('a * 'b) typlist
Q1: do it for records? sum types ? How?
Q2: provide a "type_of" function for every introduced type?
*)

343
src/lib/hsim/sim.ml
View file

@ -3,202 +3,195 @@ open Types
open Solver open Solver
open State open State
module LazySim (S : SimState) = module Sim (S : SimState) =
struct struct
include S include S
(** "Lazy" simulation of a model with any solver. *) let step_discrete s step hor fder fzer cget zset csize zsize jump reset =
let ms, ss, zin = get_mstate s, get_sstate s, get_zin s in
let ms = match zin with Some z -> zset ms z | None -> ms in
let i, now, stop = get_input s, get_now s, get_stop s in
let o, ms = step ms (i.u now) in
let s =
let s = set_zin None s in
let h = hor ms in
if h <= 0.0 then set_mstate ms s
else if now >= stop then set_idle s
else if jump ms then begin
let init = cget ms and stop = stop -. now in
let fder t = fder ms (Utils.offset i.u now t) in
let fzer t = fzer ms (Utils.offset i.u now t) in
let ivp = { fder; stop; init; size=csize } in
let zc = { init; fzer; size=zsize } in
let ss = reset (ivp, zc) ss in
let i = { i with h=i.h -. now; u=Utils.offset i.u now } in
let mode, stop, now = Continuous, i.h, 0.0 in
update ms ss (set_running ~mode ~input:i ~stop ~now s)
end else set_running ~mode:Continuous s in
Utils.dot o, s
let step_continuous s step cset fout hor =
let ms, ss = get_mstate s, get_sstate s in
let i, now, stop = get_input s, get_now s, get_stop s in
let stop = min stop (hor ms) in
let (h, f, z), ss = step ss (min stop (hor ms)) in
let ms = cset ms (f h) in
let fy t = f (now +. t) in
let fms t = cset ms (fy t) in
let fout t = fout ms (i.u (now +. t)) (fy t) in
let s, c = match z with
| None ->
if h >= stop
then set_running ~mode:Discrete ~now:h s, Discontinuous
else set_running ~now: h s, Continuous
| Some _ -> set_running ~mode:Discrete ~now:h s, Discontinuous in
let h = h -. now in
{ h; u=fout; c }, update ms ss (set_zin z s), { h; c; u=fms }
(** Simulation of a model with any solver. *)
let run let run
(HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
(DNode solver : ('y, 'yder, 'zin, 'zout) solver) (DNode s : ('y, 'yder, 'zin, 'zout) solver)
: ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
= let state = get_init model.state solver.state in = let state = get_init m.state s.state in
let step_discrete st =
let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset
m.csize m.zsize m.jump s.reset in
Some o, s in
let step_continuous st =
let o, s, _ = step_continuous st s.step m.cset m.fout m.horizon in
Some o, s in
let step s i = let step st = function
let ms, ss = get_mstate s, get_sstate s in | Some i ->
match i, is_running s with let mode, now, stop = Discrete, 0.0, i.h in
| Some i, _ -> step_discrete (set_running ~mode ~input:i ~now ~stop st)
let mode, now, stop = Discrete, 0.0, i.length in | None ->
None, set_running ~mode ~input:i ~now ~stop s if is_running st then match get_mode st with
| None, false -> None, s | Discrete -> step_discrete st
| None, true -> | Continuous -> step_continuous st
let i, now, stop = get_input s, get_now s, get_stop s in else None, st in
match get_mode s with
| Discrete -> let reset (pm, ps) st =
let o, ms = model.step ms (i.u now) in let ms = m.reset pm (get_mstate st) in
let s = let ss = s.reset ps (get_sstate st) in
let h = model.horizon ms in update ms ss (set_idle st) in
if h <= 0.0 then set_mstate ms s
else if now >= stop then set_idle s DNode { state; step; reset }
else if model.jump ms then begin
let init = model.cget ms and stop = stop -. now in let rec run_assert :
let fder t = model.fder ms (Utils.offset i now t) in 'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a ->
let fzer t = model.fzer ms (Utils.offset i now t) in (unit -> ('y, 'yder, 'zin, 'zout) solver) ->
let ivp = { fder; stop; init; size=model.csize } in ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim =
let zc = { init; fzer; size=model.zsize } in fun (HNodeA m) get_s ->
let ss = solver.reset (ivp, zc) ss in let DNode s = get_s () in
let i = { start=i.start +. now; length=i.length -. now; let al = List.map (fun a -> run_assert a get_s) m.assertions in
u=Utils.offset i now } in let state = get_init m.body.state s.state, al in
let mode, stop, now = Continuous, i.length, 0.0 in
update ms ss (set_running ~mode ~input:i ~stop ~now s) let step_discrete (st, al) =
end else set_running ~mode:Continuous s in let m=m.body in
Some { start=i.start +. now; length=0.0; u=fun _ -> o }, s 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 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 step_continuous (st, al) =
let ({ h; _ } as o), st, u =
step_continuous st s.step m.body.cset m.body.fout m.body.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 i ->
let mode, now, stop = Discrete, 0.0, i.h in
step_discrete (set_running ~mode ~input:i ~now ~stop st, al)
| None ->
if is_running st then match get_mode st with
| Discrete -> step_discrete (st, al)
| Continuous -> step_continuous (st, al)
else None, (st, al) in
let reset (pm, ps) (st, al) =
let ms = m.body.reset pm (get_mstate st) in
let ss = s.reset ps (get_sstate st) 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
let step_discrete st =
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 in
Some o, st in
let step_continuous st =
let o, st, _ = step_continuous st s.step m.cset m.fout m.horizon in
o, st in
let rec step st = function
| Some i ->
let mode, now, stop = Discrete, 0.0, i.h in
step_discrete (set_running ~mode ~input:i ~now ~stop st)
| None ->
if is_running st then match get_mode st with
| Discrete -> step_discrete st
| Continuous -> | Continuous ->
let (h, f, z), ss = solver.step ss stop in let o, st = step_continuous st in
let ms = model.cset ms (f h) in match o.c with
let start = i.start +. now in | Discontinuous -> Some o, st
let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in | Continuous ->
let out = { start; length=h -. now; u=fout } in let o', st = step st None in
let s = match z with Some (Utils.concat [o; Option.get o']), st
| None -> else None, st in
let s = if h >= stop
then set_running ~mode:Discrete ~now:h s
else set_running ~now:h s in
update ms ss s
| Some z ->
let s = set_running ~mode:Discrete ~now:h s in
update (model.zset ms z) ss s in
Some out, s in
let reset (pm, ps) s = let reset (pm, ps) st =
let ms = model.reset pm (get_mstate s) in let ms = m.reset pm (get_mstate st) in
let ss = solver.reset ps (get_sstate s) in let ss = s.reset ps (get_sstate st) in
update ms ss (set_idle s) in update ms ss (set_idle st) in
DNode { state; step; reset } DNode { state; step; reset }
(** Run the model on the given input until the end of the input or until the (** Run the model on the given input until the end of the input or until the
model stops answering. *) model stops answering. *)
let run_on model solver input use = let run_on (DNode n) input use =
let DNode sim = run model solver in let out = n.step n.state (Some input) in
let state = sim.step sim.state (Some input) in let state = match out with None, s -> s | _ -> assert false in
let state = match state with None, s -> s | _ -> assert false in
let rec loop state = let rec loop state =
let o, state = sim.step state None in let o, state = n.step state None in
match o with None -> () | Some o -> use o; loop state in match o with None -> () | Some o -> use o; loop state in
loop state loop state
(** Run the model on multiple inputs. *) (** Run the model on multiple inputs. *)
let run_on_n model solver inputs use = let run_on_n (DNode n) inputs use =
let DNode sim = run model solver in
ignore @@ List.fold_left (fun state i -> ignore @@ List.fold_left (fun state i ->
let state = match sim.step state (Some i) with let o, state = n.step state (Some i) in
| None, s -> s | _ -> assert false in begin match o with None -> () | Some o -> use o end;
let rec loop state = let rec loop state =
let o, state = sim.step state None in let o, state = n.step state None in
match o with None -> state | Some o -> use o; loop state in match o with None -> state | Some o -> use o; loop state in
loop state) sim.state inputs loop state) n.state inputs
(** Run the model autonomously until [length], or until the model stops (** Run the model autonomously until [h], or until the model stops
answering. *) answering. *)
let run_until model solver length = let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
run_on model solver { start = 0.0; length; 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 _ -> () }))
(** Run the model autonomously until [length], split in multiple [steps]. *)
let run_until_n model solver length steps =
let step = length /. (float_of_int steps) in
let inputs = List.init steps (fun s ->
let start = float_of_int s *. step in
let stop = min (float_of_int (s+1) *. step) length in
{ start; length = stop -. start; u = fun _ -> () }) in
run_on_n model solver inputs
end
module GreedySim (S : SimState) =
struct
include S
(** "Greedy" simulation of a model with an appropriate solver. *)
let run
(HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
(DNodeC solver : ('y, 'yder, 'zin, 'zout) solver_c)
: ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) greedy_sim
= let state = get_init model.state solver.state in
let rec step s i =
let ms, ss = get_mstate s, get_sstate s in
if not (is_running s) then
let mode, now, stop = Discrete, 0.0, i.length in
step (set_running ~mode ~input:i ~now ~stop s) i
else let now, stop = get_now s, get_stop s in
match get_mode s with
| Discrete ->
let o, ms = model.step ms (i.u now) in
let h = model.horizon ms in
let rest, s =
if h <= 0.0 then step (set_mstate ms s) i
else if now >= stop then [], set_idle s
else if model.jump ms then
let init = model.cget ms in
let fder t = model.fder ms (Utils.offset i now t) in
let fzer t = model.fzer ms (Utils.offset i now t) in
let ivp = { fder; stop = stop -. now; init; size = model.csize } in
let zc = { init; fzer; size = model.zsize } in
let ss = solver.reset (ivp, zc) ss in
let i = { start=i.start +. now; length=i.length -. now;
u=Utils.offset i now } in
let mode, stop, now = Continuous, i.length, 0.0 in
let s = set_running ~mode ~input:i ~stop ~now s in
step (update ms ss s) i
else step (set_running ~mode:Continuous s) i in
{ start = i.start +. now; length = 0.0; u = fun _ -> o }::rest, s
| Continuous ->
let (h, f, z), ss = solver.step ss stop in
let ss = solver.copy ss in
let ms = model.cset ms (f h) in
let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in
let out = { start = i.start +. now; length = h -. now; u = fout } in
match z with
| None ->
if h >= stop then
let s = set_running ~mode:Discrete ~now:h s in
let rest, s = step (update ms ss s) i in
out::rest, s
else
let s = set_running ~now:h s in
let rest, s = step (update ms ss s) i in
(match rest with
| [] -> [out], s
| f::rest -> Utils.compose [out;f] :: rest, s)
| Some z ->
let s = set_running ~mode:Discrete ~now:h s in
let ms = model.zset ms z in
let rest, s = step (update ms ss s) i in
out::rest, s in
let reset (mp, sp) s =
let ms = model.reset mp (get_mstate s) in
let ss = solver.reset sp (get_sstate s) in
update ms ss (set_idle s) 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 model solver input use =
let DNode sim = run model solver in
let o, _ = sim.step sim.state input in
List.iter use o
(** Run the model on multiple inputs. *)
let run_on_n model solver inputs use =
let DNode sim = run model solver in
let o, _ = List.fold_left (fun (acc, state) i ->
let o, state = sim.step state i in
o::acc, state) ([], sim.state) inputs in
List.iter use (List.concat (List.rev o))
(** Run the model autonomously until [length], or until the model stops
answering. *)
let run_until model solver length =
run_on model solver { start = 0.0; length; u = fun _ -> () }
(** Run the model autonomously until [length], split in multiple [steps]. *)
let run_until_n model solver length steps =
let step = length /. (float_of_int steps) in
let inputs = List.init steps (fun s ->
let start = float_of_int s *. step in
let stop = min (float_of_int (s+1) *. step) length in
{ start; length = stop -. start; u = fun _ -> () }) in
run_on_n model solver inputs
end end

54
src/lib/hsim/state.ml
View file

@ -15,59 +15,65 @@ module type SimState =
- Idle: waiting for input; - Idle: waiting for input;
- Running: currently integrating; in this case, we have access to the - Running: currently integrating; in this case, we have access to the
step mode, current input, timestamp and stop time. *) step mode, current input, timestamp and stop time. *)
type ('a, 'ms, 'ss) state type ('a, 'ms, 'ss, 'zin) state
(** Get the model state. *) (** Get the model state. *)
val get_mstate : ('a, 'ms, 'ss) state -> 'ms val get_mstate : ('a, 'ms, 'ss, 'zin) state -> 'ms
(** Get the solver state. *) (** Get the solver state. *)
val get_sstate : ('a, 'ms, 'ss) state -> 'ss val get_sstate : ('a, 'ms, 'ss, 'zin) state -> 'ss
(** Get the last zero-crossing value. *)
val get_zin : ('a, 'ms, 'ss, 'zin) state -> 'zin option
(** Get the current step mode. (** Get the current step mode.
⚠ Should only be called when running (see [is_running]). *) ⚠ Should only be called when running (see [is_running]). *)
val get_mode : ('a, 'ms, 'ss) state -> mode val get_mode : ('a, 'ms, 'ss, 'zin) state -> mode
(** Get the current input. (** Get the current input.
⚠ Should only be called when running (see [is_running]). *) ⚠ Should only be called when running (see [is_running]). *)
val get_input : ('a, 'ms, 'ss) state -> 'a value val get_input : ('a, 'ms, 'ss, 'zin) state -> 'a value
(** Get the current timestamp. (** Get the current timestamp.
⚠ Should only be called when running (see [is_running]). *) ⚠ Should only be called when running (see [is_running]). *)
val get_now : ('a, 'ms, 'ss) state -> time val get_now : ('a, 'ms, 'ss, 'zin) state -> time
(** Get the current stop time. (** Get the current stop time.
⚠ Should only be called when running (see [is_running]). *) ⚠ Should only be called when running (see [is_running]). *)
val get_stop : ('a, 'ms, 'ss) state -> time val get_stop : ('a, 'ms, 'ss, 'zin) state -> time
(** Build an initial state. *) (** Build an initial state. *)
val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss) state val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state
(** Is the simulation running or idle ? *) (** Is the simulation running or idle ? *)
val is_running : ('a, 'ms, 'ss) state -> bool val is_running : ('a, 'ms, 'ss, 'zin) state -> bool
(** Update the model state. *) (** Update the model state. *)
val set_mstate : 'ms -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state val set_mstate : 'ms -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the solver state. *) (** Update the solver state. *)
val set_sstate : 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state val set_sstate : 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the zero-crossing value. *)
val set_zin : 'zin option -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update both the solver and model states. *) (** Update both the solver and model states. *)
val update : 'ms -> 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state val update : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the status to running. *) (** Update the status to running. *)
val set_running : val set_running :
?mode:mode -> ?input:'a value -> ?now:time -> ?stop:time -> ?mode:mode -> ?input:'a value -> ?now:time -> ?stop:time ->
('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the status to idle. *) (** Update the status to idle. *)
val set_idle : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state val set_idle : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
end end
module type SimStateCopy = module type SimStateCopy =
sig sig
include SimState include SimState
val copy : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state val copy : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
end end
module FunctionalSimState : SimState = module FunctionalSimState : SimState =
@ -88,15 +94,17 @@ module FunctionalSimState : SimState =
(** Internal state of the simulation node: model state, solver state and (** Internal state of the simulation node: model state, solver state and
current simulation status. *) current simulation status. *)
type ('a, 'ms, 'ss) state = type ('a, 'ms, 'ss, 'zin) state =
{ status : 'a status; (** Current simulation status. *) { status : 'a status; (** Current simulation status. *)
mstate : 'ms; (** Model state. *) mstate : 'ms; (** Model state. *)
sstate : 'ss } (** Solver state. *) sstate : 'ss; (** Solver state. *)
zin : 'zin option; } (** Last zero-crossing vector *)
exception Not_running exception Not_running
let get_mstate state = state.mstate let get_mstate state = state.mstate
let get_sstate state = state.sstate let get_sstate state = state.sstate
let get_zin state = state.zin
let is_running state = let is_running state =
match state.status with Running _ -> true | Idle -> false match state.status with Running _ -> true | Idle -> false
@ -120,6 +128,7 @@ module FunctionalSimState : SimState =
let set_mstate mstate state = { state with mstate } let set_mstate mstate state = { state with mstate }
let set_sstate sstate state = { state with sstate } let set_sstate sstate state = { state with sstate }
let set_zin zin state = { state with zin }
let update mstate sstate state = { state with mstate; sstate } let update mstate sstate state = { state with mstate; sstate }
@ -132,7 +141,7 @@ module FunctionalSimState : SimState =
let get_stop s = let get_stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running match s.status with Running r -> r.stop | Idle -> raise Not_running
let get_init mstate sstate = { status = Idle; mstate; sstate } let get_init mstate sstate = { status = Idle; mstate; sstate; zin = None }
end end
module InPlaceSimState : SimState = module InPlaceSimState : SimState =
@ -146,15 +155,17 @@ module InPlaceSimState : SimState =
mutable stop : time; mutable stop : time;
} -> 'a status } -> 'a status
type ('a, 'ms, 'ss) state = type ('a, 'ms, 'ss, 'zin) state =
{ mutable status : 'a status; { mutable status : 'a status;
mutable mstate : 'ms; mutable mstate : 'ms;
mutable sstate : 'ss } mutable sstate : 'ss;
mutable zin : 'zin option }
exception Not_running exception Not_running
let get_mstate state = state.mstate let get_mstate state = state.mstate
let get_sstate state = state.sstate let get_sstate state = state.sstate
let get_zin state = state.zin
let is_running state = let is_running state =
match state.status with Running _ -> true | Idle -> false match state.status with Running _ -> true | Idle -> false
@ -179,6 +190,7 @@ module InPlaceSimState : SimState =
let set_mstate mstate state = state.mstate <- mstate; state let set_mstate mstate state = state.mstate <- mstate; state
let set_sstate sstate state = state.sstate <- sstate; state let set_sstate sstate state = state.sstate <- sstate; state
let set_zin zin state = state.zin <- zin; state
let update mstate sstate state = let update mstate sstate state =
state.mstate <- mstate; state.sstate <- sstate; state state.mstate <- mstate; state.sstate <- sstate; state
@ -192,6 +204,6 @@ module InPlaceSimState : SimState =
let get_stop s = let get_stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running match s.status with Running r -> r.stop | Idle -> raise Not_running
let get_init mstate sstate = { status = Idle; mstate; sstate } let get_init mstate sstate = { status = Idle; mstate; sstate; zin=None }
end end

98
src/lib/hsim/types.ml
View file

@ -1,69 +1,77 @@
type time = float type time = float
type continuity = Continuous | Discontinuous
(** Input and output values are functions defined on intervals. *) (** Input and output values are functions defined on intervals. *)
type 'a value = type 'a value =
{ start : time; { h : time;
length : time; (* Relative: [end = start + length]. *) u : time -> 'a; (* Defined on [[0, h]]. *)
u : time -> 'a } (* Defined on [[start, end]]. *) c : continuity } (* Continuity w.r.t. the next value. *)
(** A time signal is a sequence of possibly absent α-values (** A time signal is a sequence of possibly absent α-values
[{ start; length; u }] where: [{ h; u }] where:
- [start] and [length] are positive (possibly null) floating-point numbers; - [h : R⁺]
- [u: [0, length] -> α] *) - [u: [0, h] -> α] *)
type 'a signal = 'a value option type 'a signal = 'a value option
type ('s, 'p, 'a, 'b) drec =
{ state : 's;
step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's }
(** A discrete node. *) (** A discrete node. *)
type ('p, 'a, 'b) dnode = type ('p, 'a, 'b) dnode =
DNode : DNode : ('s, 'p, 'a, 'b) drec -> ('p, 'a, 'b) dnode
{ state : 's;
step : 's -> 'a -> 'b * 's; type ('s, 'p, 'a, 'b) drec_c =
reset : 'p -> 's -> 's; { state : 's;
} -> ('p, 'a, 'b) dnode step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's;
copy : 's -> 's }
(** A discrete node which supports a state copy. *) (** A discrete node which supports a state copy. *)
type ('p, 'a, 'b) dnode_c = type ('p, 'a, 'b) dnode_c =
DNodeC : DNodeC : ('s, 'p, 'a, 'b) drec_c -> ('p, 'a, 'b) dnode_c
{ state : 's;
step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's;
copy : 's -> 's;
} -> ('p, 'a, 'b) dnode_c
(** A continuous node. *) type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
type ('a, 'b, 'y, 'yder) cnode = { state: 's;
CNode : step : 's -> 'a -> 'b * 's; (** Discrete step function. *)
{ lsty : 'y; fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function. *)
fder : 'a -> 'y -> 'yder; fout : 's -> 'a -> 'y -> 'b; (** Continuous output function. *)
fout : 'a -> 'y -> 'b; fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *)
} -> ('a, 'b, 'y, 'yder) cnode reset : 'p -> 's -> 's; (** Reset function. *)
horizon : 's -> time; (** Next integration horizon. *)
jump : 's -> bool; (** Discontinuity flag. *)
cget : 's -> 'y; (** Get continuous state. *)
cset : 's -> 'y -> 's; (** Set continuous state. *)
zset : 's -> 'zin -> 's; (** Set zero-crossing state. *)
csize : int;
zsize : int }
(** A hybrid node. *) (** A hybrid node. *)
type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode = type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
HNode : HNode : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec ->
{ state: 's; ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
step : 's -> 'a -> 'b * 's; (** Discrete step function. *)
fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function. *) (** A hybrid node with assertions. *)
fout : 's -> 'a -> 'y -> 'b; (** Continuous output function. *) type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a =
fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *) HNodeA : {
reset : 'p -> 's -> 's; (** Reset function. *) body : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec;
horizon : 's -> time; (** Next integration horizon. *) assertions : ('p, 's, unit, 'y, 'yder, 'zin, 'zout) hnode_a list
jump : 's -> bool; (** Discontinuity flag. *) } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a
cget : 's -> 'y; (** Get continuous state. *)
cset : 's -> 'y -> 's; (** Set continuous state. *)
zset : 's -> 'zin -> 's; (** Set zero-crossing state. *)
csize : int;
zsize : int;
} -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
(** The simulation of a hybrid system is a synchronous function on streams of (** The simulation of a hybrid system is a synchronous function on streams of
functions. *) functions. *)
type ('p, 'a, 'b) lazy_sim = ('p, 'a signal, 'b signal) dnode type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode
(** Greedy simulation takes in an input and computes as many solver and (* Utils *)
subsystem steps as needed to reach the input's horizon. *)
type ('p, 'a, 'b) greedy_sim = ('p, 'a value, 'b value list) dnode
(** Utils *)
(** Consider a node with state copying as a node without state copying. *)
let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset } let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset }
(** Consider a model without assertions as a model with an empty list of
assertions. *)
let a_of_h (HNode body) = HNodeA { body; assertions=[] }
(** Consider a model without its assertions. *)
let h_of_a (HNodeA { body; _ }) = HNode body

90
src/lib/hsim/utils.ml
View file

@ -1,24 +1,80 @@
open Types open Types
(** Offset the [input.u] function by [now]. *) let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v }
let offset (input : 'a value) (now : time) : time -> 'a =
fun t -> input.u ((now -. input.start) +. t)
(** (** Offset the [u] function by [now]. *)
Concatenate functions. [ let offset (u : time -> 'a) (now : time) : time -> 'a =
^ ^ fun t -> u (t +. now)
| ---, | ---,
| ___ `--- = | _ `--- (** Concatenate functions. *)
| --' | --' let rec concat = function
+--------------> +-------------->]
*)
let rec compose = function
| [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
| [f] -> f | [f] -> f
| { start; u; _ } :: l -> | { c=Discontinuous; _ } :: _ ->
let { start=sr; length=lr; u=ur } = compose l in raise (Invalid_argument "Cannot concatenate discontinuous functions")
let sw = sr -. start in | { u; h; c=Continuous } :: l ->
let length = sw +. lr in let { h=hr; u=ur; c } = concat l in
{ start; length; u=fun t -> if t < sw then u t else ur (t -. sw) } { c; h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) }
(** Sample a function at [n] equidistant points. *)
let sample { h; u; _ } n =
let hs = h /. float_of_int n in
let rec step i =
if i > n then []
else if i = n then [(h, u h)]
else let t = float_of_int i *. hs in
(t, u t) :: step (i+1) in
if h <= 0.0 then [(0.0, u 0.0)] else step 0
(** Compose two nodes together. *)
let compose (DNode m) (DNode n) =
let state = m.state, n.state in
let step (ms, ns) i =
let v, ms = m.step ms i in
let o, ns = n.step ns v in
o, (ms, ns) in
let reset (ms, ns) (mp, np) =
let ms = m.reset ms mp in
let ns = n.reset ns np in
(ms, ns) in
DNode { state; step; reset }
(** 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
let step (ms, ns) = function
| Some i ->
let v, ms = m.step ms (Some i) in
let o, ns = n.step ns v in
o, (ms, ns)
| None ->
let o, ns = n.step ns None in
match o with Some o -> Some o, (ms, ns)
| None -> let v, ms = m.step ms None in
match v with None -> None, (ms, ns)
| Some v -> let o, ns = n.step ns (Some v) in
o, (ms, ns) in
let reset (ms, ns) (mp, np) =
let ms = m.reset ms mp in
let ns = n.reset ns np in
(ms, ns) in
DNode { state; step; reset }
(** Track the evolution of a signal in time. *)
let track : (unit, 'a signal, ('a value * time) option) dnode =
let state = 0.0 in
let step now = function
| None -> None, now
| Some i -> Some (i, now), now +. i.h in
let reset () _ = 0.0 in
DNode { state; step; reset }
(** Apply a function to a signal. *)
let map f =
let state = () in
let step () = function None -> None, () | Some v -> Some (f v), () in
let reset () () = () in
DNode { state; step; reset }

20
src/lib/solvers/illinois.ml
View file

@ -3,7 +3,9 @@
(* part of the Zelus standard library. *) (* part of the Zelus standard library. *)
(* It is implemented with in-place modification of arrays. *) (* It is implemented with in-place modification of arrays. *)
let debug = ref false let debug () =
(* false *)
!Common.Debug.debug
let printf x = Format.printf x let printf x = Format.printf x
@ -121,7 +123,7 @@ type t = {
let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c = let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c =
s.t1 <- t; s.t1 <- t;
g t1 c f1; (* fill f1, because it is immediately copied into f0 by next_mesh *) g t1 c f1; (* fill f1, because it is immediately copied into f0 by next_mesh *)
if !debug then (printf "z|---------- init(%.24e, ... ----------@." t; if debug () then (printf "z|---------- init(%.24e, ... ----------@." t;
log_limit s.f1); log_limit s.f1);
s.bothf_valid <- false s.bothf_valid <- false
@ -162,7 +164,7 @@ let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
g t c s.f1; g t c s.f1;
s.bothf_valid <- true; s.bothf_valid <- true;
if !debug then if debug () then
(printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1; (printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1;
log_limits s.f0 s.f1) log_limits s.f0 s.f1)
@ -212,7 +214,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
dky t_right 0; (* c = dky_0(t_right); update state *) dky t_right 0; (* c = dky_0(t_right); update state *)
ignore (update_roots calc_zc f_left (get_f_right f_right') roots); ignore (update_roots calc_zc f_left (get_f_right f_right') roots);
if !debug then if debug () then
(printf (printf
"z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@." "z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@."
t_left t_right ttol; t_left t_right ttol;
@ -280,20 +282,20 @@ let find ({ g = g; bothf_valid = bothf_valid;
match check_interval calc_zc f_left f_mid with match check_interval calc_zc f_left f_mid with
| SearchLeft -> | SearchLeft ->
if !debug then printf "z| (%.24e -- %.24e] %.24e@." if debug () then printf "z| (%.24e -- %.24e] %.24e@."
t_left t_mid t_right; t_left t_mid t_right;
let alpha = if i >= 1 then alpha *. 0.5 else alpha in let alpha = if i >= 1 then alpha *. 0.5 else alpha in
let n_mid = f_mid_from_f_right f_right' in let n_mid = f_mid_from_f_right f_right' in
seek (t_left, f_left, n_mid, t_mid, Some f_mid, alpha, i + 1) seek (t_left, f_left, n_mid, t_mid, Some f_mid, alpha, i + 1)
| SearchRight -> | SearchRight ->
if !debug then printf "z| %.24e (%.24e -- %.24e]@." if debug () then printf "z| %.24e (%.24e -- %.24e]@."
t_left t_mid t_right; t_left t_mid t_right;
let alpha = if i >= 1 then alpha *. 2.0 else alpha in let alpha = if i >= 1 then alpha *. 2.0 else alpha in
seek (t_mid, f_mid, f_left, t_right, f_right', alpha, i + 1) seek (t_mid, f_mid, f_left, t_right, f_right', alpha, i + 1)
| FoundMid -> | FoundMid ->
if !debug then printf "z| %.24e [%.24e] %.24e@." if debug () then printf "z| %.24e [%.24e] %.24e@."
t_left t_mid t_right; t_left t_mid t_right;
ignore (update_roots calc_zc f_left f_mid roots); ignore (update_roots calc_zc f_left f_mid roots);
let f_tmp = f_mid_from_f_right f_right' in let f_tmp = f_mid_from_f_right f_right' in
@ -303,7 +305,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
if not bothf_valid then (clear_roots roots; assert false) if not bothf_valid then (clear_roots roots; assert false)
else begin else begin
if !debug then if debug () then
printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1; printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1;
match check_interval calc_zc f0 f1 with match check_interval calc_zc f0 f1 with
@ -314,7 +316,7 @@ let find ({ g = g; bothf_valid = bothf_valid;
end end
| FoundMid -> begin | FoundMid -> begin
if !debug then printf "z| zero-crossing at limit (%.24e)@." t1; if debug () then printf "z| zero-crossing at limit (%.24e)@." t1;
ignore (update_roots calc_zc f0 f1 roots); ignore (update_roots calc_zc f0 f1 roots);
s.bothf_valid <- false; s.bothf_valid <- false;
t1 t1

10
src/lib/solvers/odexx.ml
View file

@ -51,7 +51,9 @@ module GenericODE (Butcher : BUTCHER_TABLEAU) : STATE_ODE_SOLVER =
struct (* {{{1 *) struct (* {{{1 *)
open Bigarray open Bigarray
let debug = ref false (* !Debug.debug *) let debug () =
false
(* !Common.Debug.debug *)
let pow = 1.0 /. float(Butcher.order) let pow = 1.0 /. float(Butcher.order)
@ -274,7 +276,7 @@ struct (* {{{1 *)
"odexx: step size < min step size (\n now=%.24e\n h=%.24e\n< min_step=%.24e)" "odexx: step size < min step size (\n now=%.24e\n h=%.24e\n< min_step=%.24e)"
t h s.min_step); t h s.min_step);
if !debug then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t; if debug () then Printf.printf "s|\ns|----------step(%.24e)----------\n" max_t;
let rec onestep (alreadyfailed: bool) h = let rec onestep (alreadyfailed: bool) h =
@ -288,11 +290,11 @@ struct (* {{{1 *)
let tnew = if finished then max_t else t +. h *. (mA maxK) in let tnew = if finished then max_t else t +. h *. (mA maxK) in
mapinto ynew (make_newval y k maxK); mapinto ynew (make_newval y k maxK);
f tnew ynew k.(maxK); f tnew ynew k.(maxK);
if !debug then log_step t y k.(0) tnew ynew k.(maxK); if debug () then log_step t y k.(0) tnew ynew k.(maxK);
let err = h *. calculate_error (abs_tol /. rel_tol) k y ynew in let err = h *. calculate_error (abs_tol /. rel_tol) k y ynew in
if err > rel_tol then begin if err > rel_tol then begin
if !debug then Printf.printf "s| error exceeds tolerance\n"; if debug () then Printf.printf "s| error exceeds tolerance\n";
if h <= hmin then failwith if h <= hmin then failwith
(Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e" (Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e"

7
src/lib/solvers/statefulZ.ml
View file

@ -15,7 +15,8 @@ module Functional =
vec = zmake 0 } in vec = zmake 0 } in
let reset { fzer; init; size } { vec; _ } = let reset { fzer; init; size } { vec; _ } =
let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in let fzer t cvec zout =
let zout' = fzer t cvec in blit zout' zout in
{ state = initialize size fzer init; { state = initialize size fzer init;
vec = if length vec = size then vec else zmake size } in vec = if length vec = size then vec else zmake size } in
@ -58,7 +59,9 @@ module InPlace =
(h, Some vec), s (h, Some vec), s
else (h, None), s in else (h, None), s in
let copy _ = raise Common.Errors.TODO in let copy s =
let vec = zmake (length s.vec) in
blit s.vec vec; s.vec <- vec; s in
DNodeC { state; step; reset; copy } DNodeC { state; step; reset; copy }
end end