diff --git a/doc/notes.typ b/doc/notes.typ index 3bc8777..fa0b780 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -129,13 +129,11 @@ cases: #columns(4, [#align(center, { import cetz: * import cetz-plot.plot: * - let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true) canvas({ plot( 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, - y-min: 1, y-ticks: ((4, none),), y-grid: "both", - style: p, { + y-min: 1, y-ticks: ((4, none),), y-grid: "both", { add(((0,2),(.3,2.1),(.7,2.9),(1,3)), line: "spline") add(((1,4),), mark: "o", mark-size: .03) 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, { import cetz: * import cetz-plot.plot: * - let p = palette.new(colors: (blue,blue, blue)).with(stroke: true, fill: true) canvas({ plot( 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 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 ==== 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 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 #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. ]) +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 === Signals diff --git a/src/bin/main.ml b/src/bin/main.ml index 3521d79..044e3ba 100644 --- a/src/bin/main.ml +++ b/src/bin/main.ml @@ -27,7 +27,7 @@ let opts = [ "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)"; "-debug", Arg.Set Debug.debug, "\tPrint debug information"; "-accelerate", Arg.Set accel, "\tConcatenate continuous functions"; - "-sundials", Arg.Set sundials, "\tUse sundials (does not support acceleration)"; + "-sundials", Arg.Set sundials, "\tUse sundials (doesn't support -accelerate)"; "-inplace", Arg.Set inplace, "\tUse imperative solvers"; "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)"; ] @@ -46,8 +46,6 @@ let m = | 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 st = if !inplace then (module State.InPlaceSimState : State.SimState) else (module State.FunctionalSimState : State.SimState) @@ -55,12 +53,16 @@ let sim = if !sundials then let open StatefulSundials in let c = if !inplace then InPlace.csolve else Functional.csolve in + let open StatefulZ in + let z = if !inplace then InPlace.zsolve else Functional.zsolve in let s = Solver.solver c (d_of_dc z) in let open Sim.Sim(val st) in run_until_n (Output.print !sample (run m s)) else let open StatefulRK45 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 open Sim.Sim(val st) in let n = if !accel then accelerate m s else run m (d_of_dc s) in diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 894aae5..87406bc 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -26,14 +26,16 @@ module Sim (S : SimState) = 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 - Some { h=0.0; c=Discontinuous; u=fun _ -> o }, s + Utils.dot o, s let step_continuous s step cset fout zset = 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 (h, f, z), ss = step ss stop in let ms = cset ms (f h) in - let fout t = fout ms (i.u (now +. t)) (f (now +. t)) 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 -> let s, c = if h >= stop @@ -43,48 +45,102 @@ module Sim (S : SimState) = | Some z -> let s = set_running ~mode:Discrete ~now:h s in update (zset ms z) ss s, Discontinuous in - Some { h=h -. now; u=fout; c }, s + let h = h -. now in + { h; u=fout; c }, s, { h; c; u=fms } (** Simulation of a model with any solver. *) let run - (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) - (DNode solver : ('y, 'yder, 'zin, 'zout) solver) - : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim - = let state = get_init model.state solver.state in - let step_discrete s = - step_discrete s model.step model.horizon model.fder model.fzer - model.cget model.csize model.zsize model.jump solver.reset in + (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) + (DNode s : ('y, 'yder, 'zin, 'zout) solver) + : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim + = let state = get_init m.state s.state in + let step_discrete st = + let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget + 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.zset in + Some o, s in - let step_continuous s = - step_continuous s solver.step model.cset model.fout model.zset in - - let step s = function + let step st = function | Some i -> let mode, now, stop = Discrete, 0.0, i.h in - step_discrete (set_running ~mode ~input:i ~now ~stop s) + step_discrete (set_running ~mode ~input:i ~now ~stop st) | None -> - if is_running s then match get_mode s with - | Discrete -> step_discrete s - | Continuous -> step_continuous s - else None, s in + if is_running st then match get_mode st with + | Discrete -> step_discrete st + | Continuous -> step_continuous st + else None, st in - let reset (pm, ps) s = - let ms = model.reset pm (get_mstate s) in - let ss = solver.reset ps (get_sstate s) in - update ms ss (set_idle s) in + let reset (pm, ps) st = + let ms = m.reset pm (get_mstate st) in + let ss = s.reset ps (get_sstate st) in + update ms ss (set_idle st) in + + DNode { state; step; reset } + + let rec run_assert : + 'a 'b. ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode_a -> + (unit -> ('y, 'yder, 'zin, 'zout) solver) -> + ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim = + fun (HNodeA m) get_s -> + let DNode s = get_s () in + let al = List.map (fun a -> run_assert a get_s) m.assertions in + let state = get_init m.body.state s.state, al in + + let step_discrete (st, al) = + let m=m.body in + let o, st = + step_discrete st m.step m.horizon m.fder m.fzer m.cget 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.zset 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) lazy_sim + : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim = let state = get_init m.state s.state in let step_discrete st = - step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize - m.jump s.reset in + let o, st = step_discrete st m.step m.horizon m.fder m.fzer m.cget + m.csize m.zsize m.jump s.reset in + Some o, st in let step_continuous st = - step_continuous st s.step m.cset m.fout m.zset in + let o, st, _ = step_continuous st s.step m.cset m.fout m.zset in + o, st in let rec step st = function | Some i -> @@ -95,14 +151,11 @@ module Sim (S : SimState) = | Discrete -> step_discrete st | Continuous -> let o, st = step_continuous st in - match o with - | None -> None, st - | Some { c=Discontinuous; _ } -> o, st - | Some ({ c=Continuous; _ } as o) -> + match o.c with + | Discontinuous -> Some o, st + | Continuous -> let o', st = step st None in - match o' with - | None -> assert false - | Some o' -> Some (Utils.concat [o;o']), st + Some (Utils.concat [o; Option.get o']), st else None, st in let reset (pm, ps) st = diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index 75f7543..705f20d 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -6,7 +6,7 @@ type continuity = Continuous | Discontinuous type 'a value = { h : time; u : time -> 'a; (* Defined on [[0, h]]. *) - c : continuity } + c : continuity } (* Continuity w.r.t. the next value. *) (** A time signal is a sequence of possibly absent α-values [{ h; u }] where: @@ -14,57 +14,64 @@ type 'a value = - [u: [0, h] -> α] *) 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. *) type ('p, 'a, 'b) dnode = - DNode : - { state : 's; - step : 's -> 'a -> 'b * 's; - reset : 'p -> 's -> 's; - } -> ('p, 'a, 'b) dnode + DNode : ('s, 'p, 'a, 'b) drec -> ('p, 'a, 'b) dnode + +type ('s, 'p, 'a, 'b) drec_c = + { state : 's; + step : 's -> 'a -> 'b * 's; + reset : 'p -> 's -> 's; + copy : 's -> 's } (** A discrete node which supports a state copy. *) type ('p, 'a, 'b) dnode_c = - DNodeC : - { state : 's; - step : 's -> 'a -> 'b * 's; - reset : 'p -> 's -> 's; - copy : 's -> 's; - } -> ('p, 'a, 'b) dnode_c + DNodeC : ('s, 'p, 'a, 'b) drec_c -> ('p, 'a, 'b) dnode_c -(** A continuous node. *) -type ('a, 'b, 'y, 'yder) cnode = - CNode : - { lsty : 'y; - fder : 'a -> 'y -> 'yder; - fout : 'a -> 'y -> 'b; - } -> ('a, 'b, 'y, 'yder) cnode +type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec = + { state: 's; + step : 's -> 'a -> 'b * 's; (** Discrete step function. *) + fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function. *) + fout : 's -> 'a -> 'y -> 'b; (** Continuous output function. *) + fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *) + 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. *) type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode = - HNode : - { state: 's; - step : 's -> 'a -> 'b * 's; (** Discrete step function. *) - fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function. *) - fout : 's -> 'a -> 'y -> 'b; (** Continuous output function. *) - fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *) - 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; - } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode + HNode : ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec -> + ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode + +(** A hybrid node with assertions. *) +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 (** The simulation of a hybrid system is a synchronous function on streams of 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 - subsystem steps as needed to reach the input's horizon. *) -type ('p, 'a, 'b) greedy_sim = ('p, 'a value, 'b value list) dnode - -(** Utils *) +(* 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 } + +(** 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 diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml index ddeb9bb..76313a1 100644 --- a/src/lib/hsim/utils.ml +++ b/src/lib/hsim/utils.ml @@ -1,18 +1,13 @@ open Types +let dot v = { h=0.0; c=Discontinuous; u=fun _ -> v } + (** Offset the [u] function by [now]. *) let offset (u : time -> 'a) (now : time) : time -> 'a = fun t -> u (t +. now) -(** -Concatenate functions. [ -^ ^ -| ---, | ---, -| ___ `--- = | _ `--- -| --' | --' -+--------------> +-------------->] -*) +(** Concatenate functions. *) let rec concat = function | [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [f] -> f @@ -22,6 +17,7 @@ let rec concat = function let { h=hr; u=ur; c } = concat l in { 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 = @@ -44,9 +40,10 @@ let compose (DNode m) (DNode n) = (ms, ns) in DNode { state; step; reset } -let compose_lazy - (DNode m : ('p, 'a, 'b) lazy_sim) - (DNode n : ('q, 'b, 'c) lazy_sim) +(** 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 -> diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml index 2b65544..cc8c255 100644 --- a/src/lib/solvers/statefulZ.ml +++ b/src/lib/solvers/statefulZ.ml @@ -58,7 +58,9 @@ module InPlace = (h, Some vec), s 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 } end