diff --git a/LICENSE.txt b/LICENSE.txt index e7c9f89..951a1e0 100644 --- a/LICENSE.txt +++ b/LICENSE.txt @@ -17,8 +17,8 @@ the two main principles guiding its drafting: property law, and the protection that it offers to both authors and holders of the economic rights over software. -The authors of the CeCILL (for Ce[a] C[nrs] I[nria] L[ogiciel] L[ibre]) -license are: +The authors of the CeCILL (for Ce[a] C[nrs] I[nria] L[ogiciel] L[ibre]) +license are: Commissariat à l'énergie atomique et aux énergies alternatives - CEA, a public scientific, technical and industrial research establishment, @@ -66,7 +66,7 @@ This Agreement may apply to any or all software for which the holder of the economic rights decides to submit the use thereof to its provisions. Frequently asked questions can be found on the official website of the -CeCILL licenses family (http://www.cecill.info/index.en.html) for any +CeCILL licenses family (http://www.cecill.info/index.en.html) for any necessary clarification. diff --git a/doc/notes.typ b/doc/notes.typ index 9a3ba54..59b1e61 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -13,17 +13,48 @@ = Notes on hybrid system simulation -== TODO +== Questions -Multiple different axes of research: -- abstract over the notion of internal state and of our ability to copy some or - part of it -- make free variables contained in solver closures explicit, and require the - ability to copy them for greedy simulation (without the solver itself knowing - about them?) -- think about the notion of input values with overlapping domains -- think about successive input values carrying information about possible - discontinuities +=== Internal state abstraction + +In `LazySim(S).accelerate`, we need access to the type of the simulation's +internal state in order to use the functions provided by `S`. This requires +another definition of the discrete and continuous node types where the inner +state is visible. Is the existential type really needed in the original +definitions? What guarantees does it provide? Could we do with explicit type +parameters for the state? On the other hand, could we somehow implement +`LazySim(S).accelerate` without neeeding this access? + +Should the simulation state abstraction be parameterized by a solver and model +state abstractions? What does it mean to have a functional simulation state when +the solver and/or model states are modified in place? + +We provide a state copy function which is assumed to ensure that the closures +returned by the solver remain valid after another call to the solver (i.e., that +no underlying data structure has been modified). In particular, the state copy +function has no obligation to copy the entire state, only what is needed to +ensure the preservation of this property. Its signature `sstate -> sstate` does +not accurately reflect this. Should the state abstraction be decomposed into two +parts, one relevant for this problem, and which should be copied by the `copy` +function, and another, which does not matter? + +=== Simulation semantics + +What does providing an input whose starting date is not equal to the end date of +the previous input mean (and in particular, in the case where the domains +overlap: $"end"_((n-1)) > "start"_n$)? + +We immediately reset the solvers when obtaining a new input. This is (most +likely) very inneficient, as adaptive solvers usually start with a very small +step size, and increase this step size as needed as they gain more knowledge of +the underlying function. Testing this, and thinking about how we might avoid +this automatic reset, would be beneficial. + +All tests have, up to now, been made only with "stateless" +#smallcaps([Runge-Kutta]) solvers that do not really suffer loss of precision +from resets. Attempting to run some examples with solvers such as +#smallcaps([Sundials CVODE]) would be beneficial to measure the impact that +resets have on the accuracy of the results. == The problem @@ -41,7 +72,7 @@ The main idea is as follows: a solver can be seen as a synchronous function (an inner state, a step function, and a reset function). ```ocaml type ('p, 'a, 'b) dnode = DNode : - { s: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } -> + { state: 's; step: 's -> 'a -> 'b * 's; reset : 'p -> 's -> 's } -> ('p, 'a, 'b) dnode ``` An ODE solver is, in fact, a special case of a discrete node, which is @@ -55,7 +86,7 @@ Similarly, a zero-crossing solver is initialized with a zero-crossing problem, takes as input a horizon and dense function, and returns an actual time reached and optional zero-crossing. ```ocaml -type ('y, 'zout) zc = { fzer : time -> 'y -> 'zout; yc : 'y } +type ('y, 'zout) zc = { fzer : time -> 'y -> 'zout; yc : 'y; size: int } type ('y, 'zin, 'zout) zsolver = (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode ``` @@ -81,70 +112,6 @@ type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode == The simulation -=== Modes and output format - -There are two simulation modes, depending on whether the solver provides a way -to copy its internal state as part of its interface, or more precisely, whether -the dense solution returned by the solver's step function remains valid after we -modify the solver's state (through another integration step or a full reset). - -In *_lazy_* mode, the output has format ```ocaml 'b signal = 'b value option```. -The input stream first provides a value for the simulation to run, and must then -provide as many `None` values as needed for the solver to finish integrating as -much as possible (that is, until it reaches a zero-crossing or the horizon -provided by the current input value). This mode is the one used when the solver -modifies its solution in place (for example, #smallcaps([Sundials])). Each step -of the simulation runs _one_ step of the solver, and then runs all subsystems as -much as they need to in order for them to finish using the output of the solver. -The solver output is then not needed anymore, and we can call the solver again -to progress with the simulation. - -In *_greedy_* mode, the output has format ```ocaml 'b value list```. The input -stream provides a value for the simulation to run, and the simulation then calls -the solver as much as it needs to until it reaches the end of the possible -integration period (zero-crossing or horizon). We concatenate the list of -results (all functions returned from the solver steps are continuous at their -edges and defined on more than a single instant, and we can thus create a large, -piecewise-defined function from all the results up to this point). The solver -then performs as many discrete steps as needed (i.e. we cascade). These are kept -in the output list _as is_. We cannot concatenate a discrete step (i.e. a -function defined on a singleton interval $[t,t]$) with anything else, as it may -be hidden by the (inclusive) border of the other function. For instance, - -#align(center, ```ocaml -f = concat { start = 0.0; length = 1.0; u = fun t -> t } - { start = 1.0; length = 0.0; u = fun t -> 2.0 } -```) - -#columns(2, [#{ - align(center, cetz.canvas({ - import cetz-plot.plot: * - plot(size: (1, 0.5), axis-style: "left", x-tick-step: 1, x-label: "time", - y-tick-step: 1, y-label: none, { - add(((0,0),(1,1))) - add(((1,2),), mark: "o", mark-size: .03) - }) - })) - colbreak() - $ f(t) = cases( - t & "if" 0.0 <= t <= 1.0, - 2.0 & "if" t = 1.0 + partial, - bot & "otherwise", - ) $ -}]) - -Here, what is `f.u 1.0`? We cannot differentiate between the multiple -infinitesimal steps at a given real time `t` using only floating-point numbers -as input. Therefore, we keep the list of inputs separated, apart from the -adjacent continuous functions, which we can concatenate together. Until we reach -the horizon given by the input function, we keep alternating continuous and -discrete steps. The output is a list of functions, which we can interpret as a -stream, and the subsystems (the assertions, for instance) can then be called -repeatedly until they reach the same horizon as the parent system. This has the -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 -need to reach the requested horizon. - === Steps A simulation step is either a discrete step or a continuous step. @@ -239,6 +206,118 @@ interval `[now, now + h]`, and end in one of three possible cases: }) })]) +=== Modes and output format + +There are two simulation modes, depending on whether the solver provides a way +to copy its internal state as part of its interface, or more precisely, whether +the dense solution returned by the solver's step function remains valid after we +modify the solver's state (through another integration step or a full reset). +Both modes output values representing functions defined on real intervals +```ocaml +type continuity = Continuous | Discontinuous +type 'a value = { start: time; length: time; u: time -> 'a; cont: continuity } +type 'a signal = 'a value option +``` +where `u` is defined on the interval `[start, start + length]` and `cont` +represents the relationship of a value with the next one in the stream: if +`cont = Continuous`, the next value in the stream simply extends the current +one; if `cont = Discontinuous`, the next value does not extend the current one, +either due to a jump, a discrete step or the horizon being reached. + +In *_lazy_* mode, the output has format ```ocaml 'b signal = 'b value option```. +The input stream first provides a value for the simulation to run, and must then +provide as many `None` values as needed for the solver to finish integrating as +much as possible (that is, until it reaches a zero-crossing or the horizon +provided by the current input value). This mode is the one used when the solver +modifies its solution in place (for example, #smallcaps([Sundials CVODE])). Each +step of the simulation runs _one_ step of the solver, and then runs all +subsystems as much as they need to in order for them to finish using the output +of the solver. The solver output is then not needed anymore, and we can call the +solver again to progress with the simulation. + +We can also optionally "accelerate" the simulation if given a solver supporting +state copies. That is, we can keep calling the solver as long as we keep making +continuous steps. When we reach a zero-crossing or the horizon, we stop calling +the solver and return the concatenation of all the intermediate functions. The +caller can then keep making simulation steps, which will run through all the +discrete steps, until we reach another continuous step sequence, where the next +simulation step will again run through as many continuous steps as it can, and +so on and so forth. + +#columns(3, { + import cetz: * + import cetz-plot.plot: * + align(right, canvas({ + plot( + size: (2, 1), axis-style: "left", x-tick-step: none, x-label: none, + y-tick-step: none, y-label: none, x-ticks: ((1,none),(2, none)), + x-grid: "both", { + add(((0,2),(1,3))) + add(((1,3),(2,4))) + add(((2,3.5),), mark: "o", mark-size: .03) + add(((2,3),(3,3))) + }) + })) + colbreak(); linebreak() + align(center, $ ="accelerate"=> $) + colbreak(); align(left, canvas({ + plot( + size: (2, 1), axis-style: "left", x-tick-step: none, x-label: none, + y-tick-step: none, y-label: none, x-ticks: ((2, none),), + x-grid: "both", { + add(((0,2),(2,4))); add(((0,2),)) + add(((2,3.5),), mark: "o", mark-size: .03) + add(((2,3),(3,3))) + }) + })) +}) + +In *_greedy_* mode, the output has format ```ocaml 'b value list```. The input +stream provides a value for the simulation to run, and the simulation then calls +the solver as much as it needs to until it reaches the end of the possible +integration period (zero-crossing or horizon). We concatenate the list of +results (all functions returned from the solver steps are continuous at their +edges and defined on more than a single instant, and we can thus create a large, +piecewise-defined function from all the results up to this point). The solver +then performs as many discrete steps as needed (i.e. we cascade). These are kept +in the output list _as is_. We cannot concatenate a discrete step (i.e. a +function defined on a singleton interval $[t,t]$) with anything else, as it may +be hidden by the (inclusive) border of the other function. For instance, + +#align(center, ```ocaml +f = concat { start=0.0; length=1.0; u=fun t -> t } + { start=1.0; length=0.0; u=fun t -> 2.0 } +```) + +#columns(2, [#{ + align(center, cetz.canvas({ + import cetz-plot.plot: * + plot(size: (1, 0.5), axis-style: "left", x-tick-step: 1, x-label: "time", + y-tick-step: 1, y-label: none, { + add(((0,0),(1,1))) + add(((1,2),), mark: "o", mark-size: .03) + }) + })) + colbreak() + $ f(t) = cases( + t & "if" t in [0, 1], + 2 & "if" t = 1 + partial, + bot & "otherwise", + ) $ +}]) + +Here, what is `f.u 1.0`? We cannot differentiate between the multiple +infinitesimal steps at a given real time `t` using only floating-point numbers +as input. Therefore, we keep the list of inputs separated, apart from the +adjacent continuous functions, which we can concatenate together. Until we reach +the horizon given by the input function, we keep alternating continuous and +discrete steps. The output is a list of functions, which we can interpret as a +stream, and the subsystems (the assertions, for instance) can then be called +repeatedly until they reach the same horizon as the parent system. This has the +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 +need to reach the requested horizon. + === Resets ==== Solver reset @@ -254,6 +333,37 @@ although this is less satisfactory. If a discrete step results in a state jump, we should reset the solver before continuing the integration. +===== Continuity + +When defined in terms of limits, the function $f$ is _continuous at some point_ +$c$ of its domain if $ lim_(x -> c) f(x) = f(c) $ + +The Weierstrass and Jordan definition of continuous function is as follows: +given a function $f : D -> bb(R)$ and an element $x_0 in D$, $f$ is said to be +continuous at the point $x_0$ when the following holds: +$ forall epsilon in bb(R)^star_+. exists delta in bb(R)^star_+. forall x in D. + x_0 - delta < x < x_0 + delta ==> f(x_0) - epsilon < f(x) < f(x_0) + epsilon $ + +Every differentiable function $f : (a, b) -> bb(R)$ is continuous. The +derivative $f'(x)$ of a differentiable function $f(x)$ need not be continuous. +If $f'(x)$ is continuous, $f(x)$ is said to be _continuously differentiable_. +The set of such functions is denoted $C^1((a,b))$. More generally, the set of +functions $f : Omega -> bb(R)$ (from an open interval $Omega$ to the reals) such +that $f$ is $n$ times differentiable and such that the $n$-th derivative of $f$ +is continuous is denoted $C^n (Omega)$. + +Every continuous function $f : [a,b] -> bb(R)$ is integrable. + +===== Differentiability classes + +Consider an open set $U$ on the real line and a function $f$ defined on $U$ with +real values. Let $k$ be a non-negative integer. The function $f$ is said to be +of differentiability class $C^k$ if the derivatives $f', f'', ..., f^((k))$ +exist and are continuous on $U$. If $f$ is $k$-differentiable on $U$, then it is +at least in the class $C^(k-1)$ since $f', f'', ..., f^((k-1))$ are continuous +on $U$. The function $f$ is said to be *infinitely differentiable*, *smooth*, or +of *class* $C^infinity$, if it has derivatives of all orders on $U$. + ==== Simulation reset Two possible options for the simulation reset: @@ -300,113 +410,6 @@ Two possible options for the simulation reset: makes this impossible. We thus need reset parameters for both the model and solver. -=== Steps - -The _lazy simulation_'s step function is as follows: - -```ocaml -let step s i = - let ms, ss = S.get_mstate s, S.get_sstate s in - match i, S.is_running s with - | Some i, _ -> - let mode, now, stop = Discrete, 0.0, i.length in - None, S.set_running ~mode ~input:i ~now ~stop s - | None, false -> None, s - | None, true -> - let i, now, stop = S.get_input s, S.get_now s, S.get_stop s in - match S.get_mode s with - | Discrete -> - let o, ms = model.step ms (i.u now) in - let s = - let h = model.horizon ms in - if h <= 0.0 then S.set_mstate ms s - else if now >= stop then S.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 } 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 - S.update ms ss (S.set_running ~mode ~input:i ~stop ~now s) - else S.set_running ~mode:Continuous s in - Some { start = i.start+. now; length = 0.0; u = fun _ -> o }, s - | Continuous -> - let (h, f, z), ss = solver.step ss stop in - let ms = model.cset ms (f h) in - let start = i.start +. now in - let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in - let out = { start; length=h -. now; u=fout } in - let s = match z with - | None -> - let s = if h >= stop - then S.set_running ~mode:Discrete ~now:h s - else S.set_running ~now:h s in - S.update ms ss s - | Some z -> - let s = S.set_running ~mode:Discrete ~now:h s in - S.update (model.zset ms z) ss s in - Some out, s in -``` - -The _greedy simulation_'s step function is as follows: - -```ocaml -let rec step s i = - let ms, ss = S.get_mstate s, S.get_sstate s in - if not (S.is_running s) then - let mode, now, stop = Discrete, 0.0, i.length in - step (S.set_running ~mode ~input:i ~now ~stop s) i - else let now, stop = S.get_now s, S.get_stop s in - match S.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 (S.set_mstate ms s) i - else if now >= stop then [], 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 } 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 = S.set_running ~mode ~input:i ~stop ~now s in - step (S.update ms ss s) i - else step (S.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 = S.set_running ~mode:Discrete ~now:h s in - let rest, s = step (S.update ms ss s) i in - out::rest, s - else - let s = S.set_running ~now:h s in - let rest, s = step (S.update ms ss s) i in - (match rest with - | [] -> [out], s - | f::rest -> Utils.compose [out;f] :: rest, s) - | Some z -> - let s = S.set_running ~mode:Discrete ~now:h s in - let ms = model.zset ms z in - let rest, s = step (S.update ms ss s) i in - out::rest, s in -``` - == Mathematical model #link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the diff --git a/exm/ball.ml b/exm/ball.ml index 97da4d4..9208e83 100644 --- a/exm/ball.ml +++ b/exm/ball.ml @@ -33,13 +33,13 @@ let bouncing_ball () = let step ({ zin; lx; _ } as s) _ : Zls.carray * state = let lx = if zin.{0} = 1l then 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 } in + of_array [| s.lx.{1} |], { zin=zfalse; lx; i=false } in let cget s = s.lx in let cset s lx = { s with lx } in let zset s zin = { s with zin } in let yd = Zls.cmake 4 in let zout = Zls.cmake 1 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 jump _ = true in let zsize = 1 in diff --git a/exm/sincos.ml b/exm/sincos.ml index be1d6c6..7abac36 100644 --- a/exm/sincos.ml +++ b/exm/sincos.ml @@ -13,7 +13,7 @@ let zsize = 1 let fzer _ _ _ = zout let fout _ _ y = of_array [| y.{0}; y.{1}; y.{2} |] let cget s = s.sx -let cset s lx = { s with sx = lx } +let cset s lx = { s with sx=lx } let zset s _ = s let jump _ = true let horizon _ = max_float @@ -25,11 +25,11 @@ let sinus_cosinus theta0 omega = yd.{0} <- omega *. y.{1}; yd.{1} <- -.omega *. y.{0}; yd.{2} <- 1.0; yd in let step { si; sx } _ = let sx = if si then of_array [| sin0; cos0; 0.0 |] else sx in - of_array [| sx.{0}; sx.{1}; sx.{2} |], { sx; si = false } in - let state = { sx = of_array [| sin0; cos0; 0.0 |]; si = true } in + of_array [| sx.{0}; sx.{1}; sx.{2} |], { sx; si=false } in + let state = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in let reset _ _ = state in HNode { - state; fder; fzer; fout; step; horizon; cset; cget; zset; zsize; reset; jump + state; fder; fzer; fout; step; horizon; cset; cget; zset; zsize; reset; jump } let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)" diff --git a/exm/sqrt.ml b/exm/sqrt.ml index 19cdb55..d297be8 100644 --- a/exm/sqrt.ml +++ b/exm/sqrt.ml @@ -8,7 +8,7 @@ type s = Good | Bad let with_nan = false -type state = +type state = { s_auto: s; (* active state of the automaton *) s_pos : carray; s_zin : zarray; @@ -19,7 +19,7 @@ let sqrt () = let fder s y yd = yd.{0} <- -1.0; match s.s_auto with - | Good -> + | Good -> let o = if with_nan then sqrt y.{0} else if y.{0} >= 0.0 then sqrt y.{0} else 0.0 in yd.{1} <- o @@ -27,9 +27,9 @@ let sqrt () = yd.{1} <- 0.0 in let fzero _ y zout = zout.{0} <- -. y.{0} in let fout state y = - let o = + let o = match state.s_auto with - | Good -> let o = + | Good -> let o = if with_nan then sqrt y.{0} else if y.{0} >= 0.0 then sqrt y.{0} else 0.0 in o @@ -40,28 +40,27 @@ let sqrt () = | Good -> let o = if with_nan then sqrt s_pos.{0} else if s_pos.{0} >= 0.0 then sqrt s_pos.{0} else 0.0 in - let state = - if s_zin.{0} = 1l then { state with s_auto = Bad; s_encore = true } + let state = + if s_zin.{0} = 1l then { state with s_auto=Bad; s_encore=true } else state in let pos = of_array [| state.s_pos.{0}; state.s_pos.{1} |] in - of_array [| o; state.s_pos.{0}; state.s_pos.{1} |], - { state with s_zin = zfalse; s_pos = pos } - | Bad -> + of_array [| o; state.s_pos.{0}; state.s_pos.{1} |], + { state with s_zin=zfalse; s_pos=pos } + | Bad -> let o = 42.0 in - let state = { state with s_encore = false; - s_pos = - of_array [| s_pos.{0}; 0.0 |] } in + let state = { state with s_encore=false; + s_pos=of_array [| s_pos.{0}; 0.0 |] } in of_array [| o; state.s_pos.{0}; state.s_pos.{1} |], state in let cget { s_pos; _ } = s_pos in - let cset s l_x = { s with s_pos = l_x } in - let zset s zin = { s with s_zin = zin } in + let cset s l_x = { s with s_pos=l_x } in + let zset s zin = { s with s_zin=zin } in let yd = cmake 2 in let zout = cmake 1 in let zsize = 1 in let s_init = { s_encore = false; s_auto = Good; - s_pos = of_array [| 13.3; 0.0 |]; + s_pos = of_array [| 13.3; 0.0 |]; s_zin = zmake 1 } in let reset _ _ = s_init in let jump _ = true in diff --git a/exm/vdp.ml b/exm/vdp.ml index 568e9c5..8b22067 100644 --- a/exm/vdp.ml +++ b/exm/vdp.ml @@ -23,11 +23,11 @@ let fzer _ _ _ = zout let fout _ _ y = of_array [| y.{0}; y.{1} |] let step { i; lx } _ = let lx = if i then of_array [| x0; y0 |] else lx in - of_array [| lx.{0}; lx.{1} |], { lx; i = false } + of_array [| lx.{0}; lx.{1} |], { lx; i=false } let cget s = s.lx let cset s lx = { s with lx } let zset s _ = s -let state = { lx = of_array [| x0; y0 |]; i = true } +let state = { lx=of_array [| x0; y0 |]; i=true } let reset _ _ = state let jump _ = true let horizon _ = max_float