From 1a4f950324824ab4549590f1ace0b39eacc58ab0 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Tue, 3 Jun 2025 16:57:37 +0200 Subject: [PATCH 1/5] feat: consider values without absolute timestamps --- src/bin/output.ml | 16 +++++----- src/lib/hsim/sim.ml | 72 ++++++++++++++++++------------------------- src/lib/hsim/types.ml | 31 +++++++++---------- src/lib/hsim/utils.ml | 16 +++++----- 4 files changed, 60 insertions(+), 75 deletions(-) diff --git a/src/bin/output.ml b/src/bin/output.ml index 715a739..a5aceab 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -12,17 +12,17 @@ let print_entry t y = Format.printf "\n"; flush stdout -let print samples { start; length; u } = - let step = length /. (float_of_int samples) in +let print samples { h; u } = + let step = h /. (float_of_int samples) in let rec loop i = if i > samples then () - else if i = samples then print_entry (start +. length) (u length) + else if i = samples then print_entry h (u h) else let t = float_of_int i *. step in - (print_entry (start +. t) (u t); loop (i+1)) in - if length <= 0.0 then begin Debug.print "D: "; print_entry start (u 0.0) end + (print_entry t (u t); loop (i+1)) in + if h <= 0.0 then begin Debug.print "D: "; print_entry 0.0 (u 0.0) end else begin Debug.print "C: "; loop 0 end -let print_limits { start; length; _ } = - if length <= 0.0 then Format.printf "D: % .10e\n" start - else Format.printf "C: % .10e\t% .10e\n" start (start +. length) +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 diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 346b1db..0066293 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -18,7 +18,7 @@ module LazySim (S : SimState) = let ms, ss = get_mstate s, get_sstate s in match i, is_running s with | Some i, _ -> - let mode, now, stop = Discrete, 0.0, i.length in + let mode, now, stop = Discrete, 0.0, i.h in None, set_running ~mode ~input:i ~now ~stop s | None, false -> None, s | None, true -> @@ -32,23 +32,21 @@ module LazySim (S : SimState) = else if now >= stop then set_idle s else if model.jump ms then begin let init = model.cget ms and stop = stop -. now 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 fder t = model.fder ms (Utils.offset i.u now t) in + let fzer t = model.fzer ms (Utils.offset i.u now t) in let ivp = { fder; stop; 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 i = { 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 - Some { start=i.start +. now; length=0.0; u=fun _ -> o }, s + Some { h=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 out = { h=h -. now; u=fout } in let s = match z with | None -> let s = if h >= stop @@ -71,8 +69,8 @@ module LazySim (S : SimState) = model stops answering. *) let run_on model solver input use = let DNode sim = run model solver in - let state = sim.step sim.state (Some input) in - let state = match state with None, s -> s | _ -> assert false in + let out = sim.step sim.state (Some input) in + let state = match out with None, s -> s | _ -> assert false in let rec loop state = let o, state = sim.step state None in match o with None -> () | Some o -> use o; loop state in @@ -89,19 +87,15 @@ module LazySim (S : SimState) = match o with None -> state | Some o -> use o; loop state in loop state) sim.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. *) - let run_until model solver length = - run_on model solver { start = 0.0; length; u = fun _ -> () } + let run_until model solver h = + run_on model solver { h; 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 + let h = length /. float_of_int steps in + run_on_n model solver (List.init steps (fun _ -> { h; u=fun _ -> () })) end module GreedySim (S : SimState) = @@ -118,7 +112,7 @@ module GreedySim (S : SimState) = 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 + let mode, now, stop = Discrete, 0.0, i.h 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 @@ -130,24 +124,22 @@ module GreedySim (S : SimState) = 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 fder t = model.fder ms (Utils.offset i.u now t) in + let fzer t = model.fzer ms (Utils.offset i.u 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 + let i = { h=i.h -. now; u=Utils.offset i.u now } in + let mode, stop, now = Continuous, i.h, 0.0 in + step (update ms ss (set_running ~mode ~input:i ~stop ~now s)) i else step (set_running ~mode:Continuous s) i in - { start = i.start +. now; length = 0.0; u = fun _ -> o }::rest, s + { h=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 + let out = { h=h -. now; u=fout } in match z with | None -> if h >= stop then @@ -159,7 +151,7 @@ module GreedySim (S : SimState) = let rest, s = step (update ms ss s) i in (match rest with | [] -> [out], s - | f::rest -> Utils.compose [out;f] :: rest, s) + | f::rest -> Utils.concat [out;f] :: rest, s) | Some z -> let s = set_running ~mode:Discrete ~now:h s in let ms = model.zset ms z in @@ -188,17 +180,13 @@ module GreedySim (S : SimState) = 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 + (** Run the model autonomously until [h], or until the model stops answering. *) - let run_until model solver length = - run_on model solver { start = 0.0; length; u = fun _ -> () } + let run_until model solver h = + run_on model solver { h; 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 + (** Run the model autonomously until [h], split in [n] steps. *) + let run_until_n model solver h n = + let h = h /. float_of_int n in + run_on_n model solver (List.init n (fun _ -> { h; u=fun _ -> () })) end diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index e452591..bf9007a 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -3,14 +3,13 @@ type time = float (** Input and output values are functions defined on intervals. *) type 'a value = - { start : time; - length : time; (* Relative: [end = start + length]. *) - u : time -> 'a } (* Defined on [[start, end]]. *) + { h : time; + u : time -> 'a } (* Defined on [[0, h]]. *) (** A time signal is a sequence of possibly absent α-values - [{ start; length; u }] where: - - [start] and [length] are positive (possibly null) floating-point numbers; - - [u: [0, length] -> α] *) + [{ h; u }] where: + - [h : R⁺] + - [u: [0, h] -> α] *) type 'a signal = 'a value option (** A discrete node. *) @@ -42,16 +41,16 @@ type ('a, 'b, 'y, 'yder) cnode = 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. *) + 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 diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml index e4da67f..0286344 100644 --- a/src/lib/hsim/utils.ml +++ b/src/lib/hsim/utils.ml @@ -1,9 +1,9 @@ open Types -(** Offset the [input.u] function by [now]. *) -let offset (input : 'a value) (now : time) : time -> 'a = - fun t -> input.u ((now -. input.start) +. t) +(** Offset the [u] function by [now]. *) +let offset (u : time -> 'a) (now : time) : time -> 'a = + fun t -> u (t +. now) (** Concatenate functions. [ @@ -13,12 +13,10 @@ Concatenate functions. [ | --' | --' +--------------> +-------------->] *) -let rec compose = function +let rec concat = function | [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [f] -> f - | { start; u; _ } :: l -> - let { start=sr; length=lr; u=ur } = compose l in - let sw = sr -. start in - let length = sw +. lr in - { start; length; u=fun t -> if t < sw then u t else ur (t -. sw) } + | { u; h } :: l -> + let { h=hr; u=ur } = concat l in + { h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) } From 65918ab59bb56628145257a1afe67bc23b0b893e Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Wed, 4 Jun 2025 16:26:37 +0200 Subject: [PATCH 2/5] feat: remove greedy sim, acceleration based on continuity, and composition --- src/bin/main.ml | 25 ++-- src/bin/output.ml | 14 ++- src/lib/hsim/sim.ml | 262 +++++++++++++++++------------------------- src/lib/hsim/types.ml | 4 +- src/lib/hsim/utils.ml | 67 ++++++++++- 5 files changed, 195 insertions(+), 177 deletions(-) diff --git a/src/bin/main.ml b/src/bin/main.ml index 983c89d..3521d79 100644 --- a/src/bin/main.ml +++ b/src/bin/main.ml @@ -7,7 +7,7 @@ open Types let sample = ref 10 let stop = ref 30.0 -let greedy = ref false +let accel = ref false let inplace = ref false let sundials = ref false let steps = ref 1 @@ -26,8 +26,8 @@ let opts = [ "-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)"; "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)"; "-debug", Arg.Set Debug.debug, "\tPrint debug information"; - "-greedy", Arg.Set greedy, "\tUse greedy simulation"; - "-sundials", Arg.Set sundials, "\tUse sundials (not compatible with greedy)"; + "-accelerate", Arg.Set accel, "\tConcatenate continuous functions"; + "-sundials", Arg.Set sundials, "\tUse sundials (does not support acceleration)"; "-inplace", Arg.Set inplace, "\tUse imperative solvers"; "-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)"; ] @@ -53,19 +53,18 @@ let st = if !inplace then (module State.InPlaceSimState : State.SimState) let sim = if !sundials then - if !greedy then - (Format.eprintf "Sundials does not support greedy simulation\n"; exit 2) - else - let open StatefulSundials in - let c = if !inplace then InPlace.csolve else Functional.csolve in - let s = Solver.solver c (d_of_dc z) in - let open Sim.LazySim(val st) in run_until_n m s + let open StatefulSundials in + let c = if !inplace then InPlace.csolve else Functional.csolve 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 s = Solver.solver_c c z in - if !greedy then let open Sim.GreedySim(val st) in run_until_n m s - else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s) + let open Sim.Sim(val st) in + let n = if !accel then accelerate m s else run m (d_of_dc s) in + run_until_n (Output.print !sample n) -let () = sim !stop !steps (Output.print !sample) +let () = sim !stop !steps ignore diff --git a/src/bin/output.ml b/src/bin/output.ml index a5aceab..c397aa9 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -1,8 +1,9 @@ open Hsim.Types +open Hsim.Utils open Common -let print_entry t y = +let print_entry y t = let n = Bigarray.Array1.dim y in let rec loop i = if i = n then () @@ -12,17 +13,20 @@ let print_entry t y = Format.printf "\n"; flush stdout -let print samples { h; u } = +let print_sample samples ({ h; u; _ }, now) = let step = h /. (float_of_int samples) in let rec loop i = if i > samples then () - else if i = samples then print_entry h (u h) + else if i = samples then print_entry (u h) (now +. h) else let t = float_of_int i *. step in - (print_entry t (u t); loop (i+1)) in - if h <= 0.0 then begin Debug.print "D: "; print_entry 0.0 (u 0.0) end + (print_entry (u t) (now +. t); loop (i+1)) in + if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end else begin Debug.print "C: "; loop 0 end 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, ((), ())) } diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 0066293..894aae5 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -3,60 +3,70 @@ open Types open Solver open State -module LazySim (S : SimState) = +module Sim (S : SimState) = struct include S - (** "Lazy" simulation of a model with any solver. *) + let step_discrete s step hor fder fzer cget csize zsize jump reset = + 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 o, ms = step ms (i.u now) in + let s = + 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 + Some { h=0.0; c=Discontinuous; u=fun _ -> 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 s, c = match z with + | None -> + let s, c = if h >= stop + then set_running ~mode:Discrete ~now:h s, Discontinuous + else set_running ~now: h s, Continuous in + update ms ss s, c + | 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 + + (** 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 - let step s i = - let ms, ss = get_mstate s, get_sstate s in - match i, is_running s with - | Some i, _ -> + let step_continuous s = + step_continuous s solver.step model.cset model.fout model.zset in + + let step s = function + | Some i -> let mode, now, stop = Discrete, 0.0, i.h in - None, set_running ~mode ~input:i ~now ~stop s - | None, false -> None, s - | None, true -> - let i, now, stop = get_input s, get_now s, get_stop s in - match 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 set_mstate ms s - else if now >= stop then set_idle s - else if model.jump ms then begin - let init = model.cget ms and stop = stop -. now in - let fder t = model.fder ms (Utils.offset i.u now t) in - let fzer t = model.fzer ms (Utils.offset i.u now t) in - let ivp = { fder; stop; init; size=model.csize } in - let zc = { init; fzer; size=model.zsize } in - let ss = solver.reset (ivp, zc) ss in - let i = { 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 - Some { h=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 fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in - let out = { h=h -. now; u=fout } in - let s = match z with - | None -> - 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 + step_discrete (set_running ~mode ~input:i ~now ~stop s) + | None -> + if is_running s then match get_mode s with + | Discrete -> step_discrete s + | Continuous -> step_continuous s + else None, s in let reset (pm, ps) s = let ms = model.reset pm (get_mstate s) in @@ -65,128 +75,70 @@ module LazySim (S : SimState) = 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 out = sim.step sim.state (Some input) in - let state = match out with None, s -> s | _ -> assert false in - let rec loop state = - let o, state = sim.step state None in - match o with None -> () | Some o -> use o; loop state in - loop state + 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 + = 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 step_continuous st = + step_continuous st s.step m.cset m.fout m.zset in - (** Run the model on multiple inputs. *) - let run_on_n model solver inputs use = - let DNode sim = run model solver in - ignore @@ List.fold_left (fun state i -> - let state = match sim.step state (Some i) with - | None, s -> s | _ -> assert false in - let rec loop state = - let o, state = sim.step state None in - match o with None -> state | Some o -> use o; loop state in - loop state) sim.state inputs + 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 -> + let o, st = step_continuous st in + match o with + | None -> None, st + | Some { c=Discontinuous; _ } -> o, st + | Some ({ c=Continuous; _ } as o) -> + let o', st = step st None in + match o' with + | None -> assert false + | Some o' -> Some (Utils.concat [o;o']), st + else None, st in - (** Run the model autonomously until [h], or until the model stops - answering. *) - let run_until model solver h = - run_on model solver { h; u = fun _ -> () } - - (** Run the model autonomously until [length], split in multiple [steps]. *) - let run_until_n model solver length steps = - let h = length /. float_of_int steps in - run_on_n model solver (List.init steps (fun _ -> { h; u=fun _ -> () })) - 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.h 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.u now t) in - let fzer t = model.fzer ms (Utils.offset i.u 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 = { h=i.h -. now; u=Utils.offset i.u now } in - let mode, stop, now = Continuous, i.h, 0.0 in - step (update ms ss (set_running ~mode ~input:i ~stop ~now s)) i - else step (set_running ~mode:Continuous s) i in - { h=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 = { h=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.concat [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 + 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 } (** 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 + let run_on (DNode n) input use = + let out = n.step n.state (Some input) in + let state = match out with None, s -> s | _ -> assert false in + let rec loop state = + let o, state = n.step state None in + match o with None -> () | Some o -> use o; loop state in + loop state (** Run the model on multiple inputs. *) - let run_on_n 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)) + let run_on_n (DNode n) inputs use = + ignore @@ List.fold_left (fun state i -> + let o, state = n.step state (Some i) in + begin match o with None -> () | Some o -> use o end; + let rec loop state = + let o, state = n.step state None in + match o with None -> state | Some o -> use o; loop state in + loop state) n.state inputs (** Run the model autonomously until [h], or until the model stops answering. *) - let run_until model solver h = - run_on model solver { h; u = fun _ -> () } + let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () } + + (** Run the model autonomously until [h], split in [k] steps. *) + let run_until_n n h k = + let h = h /. float_of_int k in + run_on_n n (List.init k (fun _ -> { h; c=Continuous; u=fun _ -> () })) - (** Run the model autonomously until [h], split in [n] steps. *) - let run_until_n model solver h n = - let h = h /. float_of_int n in - run_on_n model solver (List.init n (fun _ -> { h; u=fun _ -> () })) end diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index bf9007a..75f7543 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -1,10 +1,12 @@ type time = float +type continuity = Continuous | Discontinuous (** Input and output values are functions defined on intervals. *) type 'a value = { h : time; - u : time -> 'a } (* Defined on [[0, h]]. *) + u : time -> 'a; (* Defined on [[0, h]]. *) + c : continuity } (** A time signal is a sequence of possibly absent α-values [{ h; u }] where: diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml index 0286344..ddeb9bb 100644 --- a/src/lib/hsim/utils.ml +++ b/src/lib/hsim/utils.ml @@ -16,7 +16,68 @@ Concatenate functions. [ let rec concat = function | [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [f] -> f - | { u; h } :: l -> - let { h=hr; u=ur } = concat l in - { h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) } + | { c=Discontinuous; _ } :: _ -> + raise (Invalid_argument "Cannot concatenate discontinuous functions") + | { u; h; c=Continuous } :: l -> + 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) } +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 } + +let compose_lazy + (DNode m : ('p, 'a, 'b) lazy_sim) + (DNode n : ('q, 'b, 'c) lazy_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 } From 883e5fff01e66b25d41e9ab5a292c839b21f1f13 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Wed, 11 Jun 2025 12:00:36 +0200 Subject: [PATCH 3/5] feat: start of assertions --- doc/notes.typ | 195 ++++++++++++++++++++++++++++++++++- src/bin/main.ml | 8 +- src/lib/hsim/sim.ml | 121 ++++++++++++++++------ src/lib/hsim/types.ml | 87 +++++++++------- src/lib/hsim/utils.ml | 19 ++-- src/lib/solvers/statefulZ.ml | 4 +- 6 files changed, 341 insertions(+), 93 deletions(-) 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 From 589f89c76803eb4b8934858c9f9f0edbe5778a34 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Mon, 23 Jun 2025 10:06:01 +0200 Subject: [PATCH 4/5] feat: start of lift, debugging, cleanup --- doc/img/ball.png | Bin 0 -> 28572 bytes doc/notes.typ | 10 +- doc/pres.typ | 16 +-- doc/rep.typ | 192 ++++++++++++++++++++++++++++++++ exm/ball.ml | 8 +- exm/dune | 6 +- exm/sin1x.ml | 72 ++++++++++++ exm/sin1x_der.ml | 79 +++++++++++++ exm/sincos.ml | 8 +- exm/zelus/ballz/ballz.ml | 131 ++++++++++++++++++++++ exm/zelus/ballz/ballz.ml.bak | 137 +++++++++++++++++++++++ exm/zelus/ballz/ballz.ml.bak2 | 135 ++++++++++++++++++++++ exm/zelus/ballz/ballz.zci | Bin 0 -> 6309 bytes exm/zelus/ballz/ballz.zls | 19 ++++ exm/zelus/sin_1_x.ml | 74 ++++++++++++ exm/zelus/sin_1_x.zci | Bin 0 -> 9262 bytes exm/zelus/sin_1_x.zls | 7 ++ exm/zelus/sincos/sincosz.ml | 45 ++++++++ exm/zelus/sincos/sincosz.zci | Bin 0 -> 1170 bytes exm/zelus/sincos/sincosz.zls | 4 + src/bin/lift.ml | 119 ++++++++++++++++++++ src/bin/main.ml | 42 +++++-- src/bin/output.ml | 38 ++++++- src/lib/common/debug.ml | 7 ++ src/lib/common/ztypes.ml | 151 +++++++++++++++++++++++++ src/lib/hsim/sim.ml | 8 +- src/lib/hsim/types.ml | 2 +- src/lib/solvers/illinois.ml | 21 ++-- src/lib/solvers/odexx.ml | 10 +- src/lib/solvers/statefulRK45.ml | 2 + src/lib/solvers/statefulZ.ml | 5 +- 31 files changed, 1297 insertions(+), 51 deletions(-) create mode 100644 doc/img/ball.png create mode 100644 doc/rep.typ create mode 100644 exm/sin1x.ml create mode 100644 exm/sin1x_der.ml create mode 100644 exm/zelus/ballz/ballz.ml create mode 100644 exm/zelus/ballz/ballz.ml.bak create mode 100644 exm/zelus/ballz/ballz.ml.bak2 create mode 100644 exm/zelus/ballz/ballz.zci create mode 100644 exm/zelus/ballz/ballz.zls create mode 100644 exm/zelus/sin_1_x.ml create mode 100644 exm/zelus/sin_1_x.zci create mode 100644 exm/zelus/sin_1_x.zls create mode 100644 exm/zelus/sincos/sincosz.ml create mode 100644 exm/zelus/sincos/sincosz.zci create mode 100644 exm/zelus/sincos/sincosz.zls create mode 100644 src/bin/lift.ml create mode 100644 src/lib/common/ztypes.ml diff --git a/doc/img/ball.png b/doc/img/ball.png new file mode 100644 index 0000000000000000000000000000000000000000..2d578f451610daa8fe0263d07157906b337dba8a GIT binary patch literal 28572 zcmeAS@N?(olHy`uVBq!ia0y~yU}|7sV0^&A#=yW}dhyN^1_lPs0*}aI1_r((Aj~*b zn@^g7L4m>3#WAE}&YQcHHK9-E+J88I_xhbXZ>kk^I8SJGaR&(seqY+XZE4w@YgxGq zZ$w^+Pha!)!pd7>(bkt&`stSCEZthVDj{~lu14;zhBw{}ii)q#&bxD`dh+gn{}L7~ zYSB5r^SjOSIp=?W{q)tU?a3VdAEq_McF!yOG^ctkb#QTUacN5o+s&Y)q@;98hpz!7 zz;>8%LPtl(i3AG$dYW$Y$rP)n^ZD(5Ffa&oAN}xRaevakKR+WjrEr?%-dgfxDaeA1mNj1U za(;h%Tk+u_yJ7jeJCo(5>d!60WWajo4kL`)>dLIukRqq>qn|-q@JTuI4{)O$sN- z{U$uYT2s9=3l|kw@bn*VmbETRxV9#8`JX5akwa57gFjpe_D}rx=O;tM?)Uqo>-Nik zzgN9pXu7(RlF}lRt*&Fo+(JS_ z=l=OK(>VRnQSo@0Pc{6uUj*j;y*GP9gpS7&p$;pj7scC?4qV~yKFYKIKL7UH%)9Th zF?6^k+dTWIdHi7k2g40@8@bPg{{H-rK3aVG{CVk$_T#s#-@SOza5by-{`={R%V5_!t$skPP*HnU|D$s& zJ(N1vy|G#nkqU~D@0;AGvWId|VVczSB~U6HkU`j~e;>&D0Fu^V(}T1*;@PiTjIGnmd^9eCmGJ)!}=>ww0l#*-R{Q$l~eg&+lBWt`oi5kUKL7up=O0(c|K&1HJ9B_tzNX>%ylTGt zKhKpv_*hYJJ+?gY>8YtZ-tYT;;mVbkK6(cJrBQCCad&?sHMFMeHFJi>HD_v%8 zKgkpG@iSHi-5ec(sOsqwbx(!svgDuSS_@tB-z%i*1Vz8{)us;n|by%9hJ#{PR>5oXA6$A zES=96JT$(k2`tr(*&(npgiG2iN8;-BrIQZK^6r1I=JscUr_bJe5pj2MahbxHJh#1V z@obJCER!X7e0o`U<8S=i7q(KNyJbMuzVxtWa9eEXvUsCm?b_WJnz-NTAA7ysR6#O) zhry*DhXuTA*D`KBt-zbiKT+eGK&Q~m_H!<`(iz^C%?@}LTewpE__fvzg_+!9c43PW zQkRIjy0&IbcaS;`@vW*AKdjtkG-ybV|k(G70+g+PA-O0tpMTIw$-9a_4a@Lsx zfBx$n-NwF4pxWmH+>aoy6qC& z(6ITvmc#XXj$3lMSI2UJ%F9&;IX` zVgG;OFsmY{jLx{k*YI-3)0dhX7re01-}Yz1_Q$;q2zNiyN@$(4T}fU)rc`9&Hcb;9b=UDsq{ih*j`hd*b^^}t#{eLjyE`B$GoCn=aUwaUy zCb6th*etr@SB}}Xbs2{bg)Yt9$+of|MSg?i{WH$@ComI+^LoQi{mQs(|gr9knOE z#Sag>-*r{5wWFg$BX@}`!_|%k#)Mnu#RWOC;mofYQ9{f4(=v;MjnNmj<_qjhw=}A& z4dO$QZ(e2h_TahyS>d~RKo#MZ31ti++DvQ>o#xM!Dl{4^THdo?GeUL4%{yv8R;{u4 zlOwC_Ng&QbA5*qNi{ANb2`HC`0P$l0rsf;1Snu)F94X0Ve z8L=(7@&Q}-wYqhp*fy!_bnT0~;;vdN%Tz%oo|;g` zFh!AxjUjv<>$Du>S@)l9Ziqrj8(|lJtd6nxvfOS;?3|@#nkHPxNh?7#H14sx$^@55 z(M&>uf)n{XtQk%TF)}w?&Trg(qUYFkN#<-$)2r}fsx@cOv8xXe;}%b*d>!o$EfUvis`G}G+gprAo29FkacIy39vAKv&i zaZ1dS(A^?ZY$(naVtrmIc*E*IIkT{!;KVdfYlc%2jLgUH^S;sAqOieqO_`<%4~o-V zy7UVcW;(MReZKRmUh6qjt0ta0z4^kvR&K|Rjt-C6OJo^Rdm9)Ps+)K>G8$LCezK;& zFZqiTs1X1vQzaO^mzB&DJ|4%Hbj-Zys$Of`0Yu)-)ML80zuvx2&UV(i%d6MDDbk;1 z-f~5+xuc^a+pkVXtoz|bclo1#etzDVa&l5|-P@9sz6Qnx zFLK1$s)Z%qJgqa>l!UO@M40h-lHr}=^R^i|IW|upuMS&#=+@Tk!vBB2FW$IOuv=XJ zSf8x5LFK0>m*=^EjdzNo`0v%5|sqyKThAq@$Y0 zrOoqP6fTr~`h3+w&*FUZw6%ZL?{h6Gv$e7N`oramlDDnRwR?pM7ux>*&R)6_)G)Ox zSnw@{MaS9JX4>b9x3{*=pHb{5BX51J;oZAD4d0E78}{#)FTYo*UiYy(e!-qSf2PF0 ze)*D-!D8;Y$@0a!@2*(8_Go~Qny^qOOVG+IlK(9odz0JG^Hi5{y|q?eb@jmGh5g-K zT!Aq$D!yQ=0|OO{ivC;zw>W?QtrHO!_YCIme`@2W);_Oz(n%H{HQ^~~UbS187(wms zYwP3t=l^-4ZZP{Ss2N=F_v`h?m*@X8$uW~&wR-i%D_2@BT?%42@apR7k8d`ge{?f_ zzU|weVtO$hAzH1x(q ze@dQho`39EuQbS%MSJ$heEs@Wo&7-V?{9NI?UWbqJ+|ON&aAM|`Ck){_wPL9?bhfZ z9klYw$@%*_-MS&oAuW)hf*+LoZHi9mo<6O=UuJdq`eVPozJ7>kR^HfGYaLhdkk#hj zkHk zKTq^px*_wjb=GvFHM{nIdGLZeQKx22$lIV7%BTspMtlR4m)^&P-%i>snRMK|L(Q^) zISkfRSh};EQM&VRx|>m^)AM_dTfG&Pl$4~VZDC&U<~`HD9hbcKWYj6ID3kmm*vN4I zpl8wzXK?e#+R=>@tRBz!l|8Yg{=m$i2Xt*7a z3bZ=kI&J$-ftHSrjwh;b8CD$SkYM0C_kGcGhWi{!%%}$K;MUn-v6jQ^*u+r1xAN;j zu0pg7OjXj`WOd)STm=_XFSaa^W!Ngk$n0>_MtxS?+wk2yhfzwuY>$6C5_RSyMwO#0fq^-S3BoL?GttJMX9WMSME90`%Cwe%J6ZaZ9;i}hLepg?xLb00d>G^L9 za=U=rIWW`UGdK}j=x8&lbSR@vQ3kWLJn`dfw-+c_i{idN0dxP6Fsp>Zop zhG9PXS=M5Q)rvC7Hk7<~t?|!IIi8XN=L}HOZTc4G1#-#>%kT4~@8kWeTe#KRkOd_i z{B9JleQ`ssZF`Hf({gE)@K*QI->~@QS5Oa4NlA%$av8&wWELKVf17PyvYeCwVMBuYGYx64D)!?M=DCYx&?QrTPiux(kq3GB%eI%Vie>8!QsQEfUfjVicmMds zsn1)M2i%(d!o|hKc~Tiek6ZGtz04WkXC<^r)*es#5PIv1p8IyDMzbYBnj+%jy+36h zIxzhXFJJihucJyb=Z(K&7PVnZgdCPEX5QE+64P4d^RG!fb$ZIdBb5Q+Q(h<6;qXndz4WzG_~dzbj_x;u?8w?9vK)LVD(l&4tFdvN#0km=gYJ&|T5?>ZN~s@kf( zu>W>ZJ^QSq`}XC3qg-5EK>c77Czd(3_u4{~IuFME)hi3UJaf_XQuC|H zr*e0$*cJ8jM!5!e_|(Oj$pY>F@66tqJxKx@NI&!v<`#+Fx=?M=bg1}Jz4)n&wy4_g z_iEp6_^q?zXHuGZj#+f0xBTh@tK%)(HnGO-;L?%YYu76mUN8yNeZ6!<>%m;>@<-F6 z^A28F8T@hT`Wja2@^>~rvo}WY^xOT?cy@k%zpD4Ngv-l(z2c%a7o96w^z`k9Td5Lm zk5@aEynS>0(f!aBMYZ=Uq|g1#Q;hrj+il6DRZ|>XT$bEnxN+@TTlV_BVxOO#JuK{R z)A;IDR!sHVtr@wwtPE?SwsK`(Uw8EH_xt@IHJ~Q+Lg)52oB4uD2@EIZCpw4g&S&;% z?|0X~m>~52Z#27J#fp%t<=^^t>u~S$G|Sn!-*>N4hKs~~@TkuuJEnWva-~n5I`yG9 zzlNE?p!8M9o$B{{Eoy!kY}&Nx!vSXg2P>D)JM`!0XGR9=vNs1V``h`^9~;*{%{6h#ndCqE+R%Vs z(%&GiAJ=1kzef0BfrNGWyCd`e|4E-y|L^AvyIQNL=;(_#Zgi|#wQBNw`{HLlZ|?1_ zKJmG)90M%k8$MqVseQ1rE}vgBjwbJyQkAXHsXV(%Tc&Mh?QHqP z-#a@9ZKXHC+ z6!}+lH8AGeZLS|$U-y2r6cp@@+x1NG=sO0xKOYW9Y|FWMvi^I84alu{^W%L}yzUElbB8UDWgiMO61p>&qP=V;4{ z-}7f(0ks@>TrQ?F{Qq&>{*izEFX3%DH=Eq$Dp}U=`=uqO8^zLpym|M#UEJE?>yCVT zd;8-dZv6)zkINtTHNR`}J}H{?N^jZ?h@wMlK425cSfJAsr{ah8-)a{hI&9i3BQaqnNCG-HpgAX$kN@oe?^eC^6dF1Ju zoU}D>M=i5P{FJ3-nulk7FW@w*+bk;DUArnT(0TbIEYjYSZl1e0e$UqXk+#fW#?`fDn#yGdn~&c3KU?mI^06JC z*Rnz@wx!u}lgb#j+|J!_yLYpeQT{gOzwW|^>y>!5CE5A@ZexG`e}AEdo?g)IIfBK{ zW=*uc+BtXsn~>cizwBLITW7seoWx(<;Q}$lWY(cbxg)=J2t{!A=CcOYYz&`#KDh4H z!X%UZ*ZM1FlrwgAXw)u|W#|=SynUzPYNlg-am?Ew6J0sa#Mk?`CZD_eb5>8BPMNTm z-nm29+HYOeYd!Fz@5sh>k11bPH%v`i5_&*AqBgDMJ@4V@Q>_ZgS9dIw7wz6t@Z+pT zlUexvv;U6A?D=!EN`)07z|>8jW+ zfz_8=Ho4xs@wcN)Y=eN{#BCnd3{tvGY<>DHoN|dQHetuFeOJn#CVlVHWAklWW%l~S zO=|rmUOgejdWF|%TbH5){)rlH2iPXqpY6>spZ`~?J4EQl$IE^(-%_5Q?9a5U`Zz6W zVtF4s`@A>DBIBR#Ysk7-#6M96)H@U0z-UpUW>K8WvAxiG#>K2ETaN!G?0Zr#{qr%- z+Pr^m|F87wiC(K>IRyn@_HJQbAgG+cu;F>sg<^$8Ws5&n_Evc7 zZg?JD_-5*YTT3OTU90Z60vdO}!*NT&Cbi6D|LInR#ez2k9WTkaUyS*#DS!Pi+nRLU zvq$&E7aDd%eYxPgMDLqG=OecZY=H*H`}Q3^w)?Z@qX)vaGya{m*?QWed(zIl11u$P zp7FD~xGa&$WOrEXc3|h_J03@W@;vZpW#%k-C*^Q-Y2^GLrh8rVh3%go-4`00H9gUz z>ZM{3C|MMWZgAXDbjD+*FUMB32b$Yf#P)s3`hTH$)%lczHM9J#oqn~u)BIV||MmSz z3$z*2^{SO7f4StmL=NOAt%St}jAreed)s#iq*$y}x$q-P;HvkRy}n&X`JT*uf9&4z zA36$xf||~^7%!|4G+_9VJ?DjW)|vaiGfeV1ly=|Z4FC7?eb-CPs+n7L_<7f}?^N|^ zd%3dRLl)F!4L!g%$^2^9G~dLnrzhN-=U!#-)Z2N#d8Gg9_nWlVdL>=e=Pp3?c?V!bgMM?A8G%W9LV#{i2a<%-me0wa@C+QaWjeT z!#(e6dm8=Kr?Skt&z1Ox*ZKJVhJTC1H;Zb-zB$^cUUfz9sKJ+O&QmHuuIIRsP_)eY zVzI@lGRdH<8v-p?vn7i@w$Jn3rzu^kZ_K8qq-1!BufeUafl*=kt)_3eHNOwPl=&Gf zzv}#kACFmV-Sqc;G_xpPlX3f8Zj+H?zpcwm&|pVT1Me~OldT>K>&hfANxtw}@I6E1 z@j8h*hAY`$jNWmXf3v%{vP`mV!F9PsXA&W)e{%3!hc)(nD(MG;W#TS3FXFzly!LZ@ zU`VYR`@>@_)y*9p$9!%vUU(t6;U?FO%{edZOjeX>265&{9GJ3S{^o4^FH(CgYSJpY z-mHkezPwCx*K(NI`?zm3$heBJ@bfF=uX27HwB$muMw>^v%)icWg0=hR-9H#vn{Ad= zf47%;+9>Aj@{>b^=@S!C(Vbg1qd{{ z-dmv11s?0^=V@k)6qFWSoLI(?!csHG_IA&` z!^ip>gFP2~|G;o8-h9{d|2+bGpUi)J=;898SvyJ$o%jQ*K(i|5O=cSLTZL=gdu%zj zzwbQ|p4zURU$VZ=_z>Hgso`l0GnwZz`-m%l5m+i--QiNGoY8!0W`^VI-^cbbSG+hW z=kudv+O1R%&GIhUJ{jvOyC>C)J3BI-ERki9WR!nC`_Z9<<%RPZeN-Fg{S6eE^IGof zbt%=VkJD~n{1CPOR_cpMfmI!s+HGAHNia%>e5;V2;Gt-@piK2uMx8#3`#Mp+`!Axe zYFs>hW&O{KzpU5~lo)vN2mad7?jfR_z+C(!@?!Ok8NvzY*E!DlTAEi}q;@mZjHk`_ z#oEWR#Y#$2l9}udZU^dedy<^~~bq!_9lMrZ2R1 zd!^{b4=Oa=F0h4JTyN{EEexr3d+TMiIyGA&tevy${f#<5!%J*yOwE69eCz(l+bt7W z%QMOM22bJ|de^rEY?`pT-u;mM&9(}G7q+K%KXY+0v4OZ^viZ^8i1{DyeY>z*z3F7U z%-8k?j?!ZIiYuXR1`j0To<~%p&V(z~` zfo{>Y=9QiDXY*UG=sAb?TytJx19GxfLbI`@?+e?MITAk}$X5lHG`hbvGEuFcCA=|! zsd=v_c);(Zr!|8Slik$aH~IoL7bZP;aFEgG$Cn$A;}29^>^{l##IO9}jrg7F3xZ(A zCb8rc94*e}s7dHqoAB#ZGRJqd388x1JnnW}^gb$@XR_}UV}PdH8^s_oNC~v#&dfZO z-@SrwWh)o*@4WmbYx+V@Rq?ci`!%~*;z9Eo=?ZTdPH_APSy!X&^wv;?cLLM)y`l$? zoKNGPf4t*j_ffuz(>WRUPwj5-GDzZ|xCxXMr5&Z`T+M8Umf<-#zr)=oQ-<(dj!2zvwH9od3Kkw(FVX)^-m|P)_Mxz^N&Gq1@yCdzQR=7ufx0 za;Ih1n@o(qxzw)Y>Z{%U9UUG$TbLD;ANXs(t!wec_1ucqw3Vz&}w&;}oQ@7*WHrz@T5S+MV$=})2 z1V5-(dh2)`Ynb8pKW=@sS;wW;$ga8PPh9l=r=u>!{2Dw=fbnn>;a{7i_ZUNjH!M8Wp?6u-1Ul{^Jb7Uz=3?( z^Ba@4)7j<^BCnONWk=lapL6rn{lCdym_+7*dSZGhL6f|LQot z_@Kz6TFxiu3$*IqIV|8^8FK1Fb-~oSQ;mnuw*(d~G+KYdc+$hF74IW>VyoE&c#lq$ zVpwB$Ao``tzev4^9g77U@3(uz?0Tl?%6{P3F}FW|p4%VKUca}k@bR&Seb(;|goekm zGIShGirAVZTKE6={l$wGF`b^Se|%5n=QbO0LoEy5sMhY`Tk>?T!kwjMGcFbA z&$@WJP?fvpM}hC=Td5hBURJHxuh}5(Q6Y>hY8u zdvI#|qF-t9TH?P$R?JB`_}*CsJn6xl5FanUZs#+pn8Kr?6;CF*Z%98c7gPK7YQ^id z+YRdf)gX@q96RQw6SqeK)Ioo8a`IHIg9e|L?GTiI5wM%b)nu<>L)Yobbgj8Dm10XK zseO^{J{W#-`|F<2CAN*<-|;T^ewZ;NZenlb{!{JSzkIpxBpmqQhW6eY|5wMo`IE_? z_9vO|Sp32_X6`O7N+P@8Z_2%GcIEo@$L#h$81?piaJp0X`|Xae*P=h3(q6A1!u7C` zUG6~b_q*+^++q)&&CWk|X{qVrwZO^ZcYAQy!O(S z#Zy6j=ilG2Cw;Fkf4J2D`oS={qhDXEC!DjbCN>o4^;l2*MFD4 z{LxQoeq}K~&yT+!6+Y|Zmsw%ElPjXFszT@a$6H&o?G2yLv9K#!!+7`ZsrUa^-(avP zeHF6j<1y(1tA6{_qMa@UPfiF?$~e1H~A(@d8exGbNi2cevx(iv*h~sA8Wn0H!$G!v{%WaPzYUKt^WhEYFTYCC(>stP@1F>p?GsBJDA zA3RLj>V0u>zmG<^Se{~{Vr&cV&1NoYJWaW-`{fndGq=5KZ-979y~aC zzI3OH-%jNpg?Bi9ohe!#d7)Op>W{kDnwMuej_G>ck85jWFi)^O_{{9!E7Qfr+v~RR zdo6G5V~P74#LOi&)9{_y(>gajtK+TO^W@rN`AZjlJi8$2xoYLd-vXVJ>Xcs;SoB2a z?PPsEw_MIF|6Wbe3+9GPmx88DpT2nCKDnh!mohZKm%sI|Za9C&`d(+_)?0C$GMSx= zE9!p9^*n#IJNT(z`NGm-iJzON+%IzXUsUV4r>t7}+^)O*tE8vw_nsb_cr^5U@8Q+) z*-lGz8&8`565Y>sviUxwgKM61qv7mq^T+3`-ydO@uPONS_~%dXSl7E=gQb zw_JYKtXUs^Jnnyd*8G0Rt50o{jZd~7jC*@e@vW6W>&>jD%ja#3g1huf9-R39Y2D?S zDvuBRn%%eXkCuN_+kWAgpJFk;?k|n_-^#Y*ldi`87VVZ!jm;v`4E!<{4d>4J#Z*3> zy1eGkzTfW-e|dR%W68^)8E4aUqPOu(ojNsQa~f|}R#uyhym^zI@oec0);cfp{SHKD zN&e$jx%?^5p6l7N@1bR&DEPZ#&r7{G|JLiNJhs?5q0#T_t?>SZFOvAr{cC4pc==Y* zzf-ikF1l^MUgm*E_pKkyKQVXl*N*$`9wED)E!?8R6v!mr$diBgQreXoO=;E7Tfa+f zy`GnSWJ>g=z%s-4H7S?k|L+qt-aobb(ayZB-V@K={(q_K-{QcSpT|%2#BG1K)Gp_L zp!&LKd$}<2``3M3mR~xeBydW{f1>DJ_O06mRHm{hHL7quD`0-%qdn`f+ws1RTd5w) zOTAg{?pE4SQ<~Op&md91m0xvw<%d_{(nfW;N=AB_><;UgZ=TKms5Kf7+pP#3bqwsHHut4%-u2mb2JunD#;IryDxq8CGT-P%8C z>zC-&$twk3TF(?F70$?iKK2d6qg%U~Z+!3)`1yRv|9A4?_iKavdYL;tggP}69xA6Z zH1NjfyPREXd1B_1A4k@wZfY0%yG-VMY@x{0=~X{$7W94i4V~+db&!g?S7LrB_m(t+ z^)~(mFJ$&j4&zPJo?v>jZ{6No_nbU5tL0v^oKZWo$N4Ys)whvpbN>481m(l5=}sw>9zHVNCf_RCTXNx^K@7`#-{J6Q@+PKWuOCRLFm$TZJ$TGYwlYUcnPt{6UIN@{j#qyMbALpmK z9G}N{qgE`!c!T)I&T_`fD^KZtRjBE)kJ-C(598JP>FnhXFZdr?anX6YO;}^{`R_`X z=AC#^s%iV97B;@LWMbitvoCI~o_oD|>mEs`Z&SX!`|qAU&HQR-K(I*R+bzA1`qsX7 z-g3R?<^2iob!E)9b57njrMdU?Y|wyWiO(&@34$L^{78?wSnOe4@3*&o(TY7f#>@(f zr(ZUXX3Ss8Y|GxeQ_5m0%eMb(q@Vxak#OUD=<4=m_op>S$Y1!9rt$c?S}LftNMi9> zJ~i~I!;VmC<^(NH7pp}YGbYdMT`jqnXYCIq@mXe4#fm5T+N~$rpKUxe^YDL#MQQx8 zJ2nYK@4MbNrA+SQ(#sFND++;{GTP~rq&F@wV>@uJw%K5&&z(9Krq<%SFYYa$o9;eU zNZ#Y&ja_zg{wx3d?{sedvzGrl*57O|E-lk^S>gdt+N#_I&=&q_JCvJM8%hGsW=T?Rq~F4f)o-bGA|X|NnkT zRn~WbOnK!5sY$xGo9;P(vOj-yqiIvuX`}CaQEvLqn{P1AKHs%2>hT(um$pLA@BeX} z`?q|LV14f`+b3sC7e7hT-}a|^c9!V`iLKQgTlxhJcuxi=xxW|ung8@iRZz|U^)9FM zzAkvG`&h^5#;5m;A$K3NADwqETtMsXo83i%H$)swPP8+2dRRidqqqB5)4Qmh|JUr| zbD3zQ{Pz9-zZ!1(>Cfh8IsXnzwqcjxRh!tOV)%FV;_AokJhmy{X9bolWWQZ^WnGr6 ziJ8$CflTT9tUk`4e2W%6T^qOcw8i3WYyW=mJ@Muh|GwtWXS$wMY_z}F8zN-j`FrsK zqq2YQ>S=$A|E;K>zv;^b%kR$~F3s*&l=3+=Z|Ud%5mw6H{1eMQ#vU+PGUue<3?EO2 zN0mEEgISh}r#|>9Z;}+IKEY?^g<6k)$Jt$2`j1}x;3dGCp8oD<@Y0WZ3WA!xki@pL z=xj&ZqzMFdJ*50{t;nVrANgKTP3x2r9ZXIZJ-+0;XdZSr?mg<<~zHv8NSElN+q!bb{$7Ayx zPm5WGm#~!`2`I@+yZ68PRH@`%p3^VQ53_LpivD6A=+f+2@`fYv*{6x`e}6CB^mW4W z8(jZc)cMxG4BIUtrFgTdV@tcB!Q8-eEY?bIi~gzJi?(h!#BKEX-~IliFGfb^inMRt zzka{E@%x`u+&?$Z*-(A_U#{OVxmiE?D_7Jma(`=-(US5;@sgO^ffBcSTI_79G7O6k ziY(g8`6TQ!>w!(Hb{uK_7F%@BR&I5OP{yqb{hQzG1?(0%`qXFJ;fLE7KK#1i)>46q zZQyRSoNWA3-jj_}O3t%w&E-EByXf49<1bQEip~`+4o}x^d_C`dqwQ^;^$#{L-}LRm zmFUeI-|v@X&Td>|>vHQL$Bc6)QrOP$3IFOl`PHHR`#cSmTJ6eax> zrKVIUpA5Y9-^D(yuJ{coPO|w!bxt{%)YP=vp7gyYr}Q>~#pm#;SsAV#n;18g21~FO z%VbS=JX|1p^TN((`P&Ce-pI0?;hb{j{Ga0GGeb)EitmC%Ne>Pq)3=*k=A^qrj4H<+c*_Tm2^%pS-{a9v^P(1sJo~m-zzB$c@_owstxGYu&4I@;1`6C=F@n+)l z^Ew~mSW?AzUfpfCc)H+bE5;dTUff#$_HDweYb^$kXaDS)>;Fye)cj`-r&1i-&u#m7 zQTS$c$CK&nDm`x9HZI%EyK>4m<&{30b~3hKZVFJ~^q=@^^Ir#r#yy*_N;7nKafR;Q zCie7f$)e?D0+pAy?95o6dtD*uB+I9Cjqv?qRT+EFLQj4N4#4rChFYnFTXLUD(!@%F6sXv?7wDT{a>7YU$SX)z}<)a7f-)P3AesjAM*E}WVegRi#x&x zcW6(I{9D~5R&`=(c4Ph=X8Cn*gLjKq@#=mNIN2*`(Ccxn!gv?|s(S4O-y*;dI#V<4YvcGGDC+CNhF{l+r zy;PheuACq>am`0Rm#8cExt0}|v$Q^6x#d51Z!c#-|MrE|^S>oMaQ-iS?8<)gcOQa{ zt~ZB>aIgzyO?PsUQUu3S&YsjVm#C|Mi@lFE7DS6!PQD_`Ai1|E&aapGt@YLaj7cU+ zZ#h-iCTq-mANuD`)T@PuKHFG)SuVHr^n(72LHrZv><~2R?bv>#DTS+J^$xxzuMa4! z$@^;7VAm+Up)c+pkGfU)%KaS~Cg1EXHt%_9HsjKVyFT8wCNH_|O!l37&XBrzshpB! zp>l?l;qv=!^PB~*uDie2zGu$Lva9R-lfw;8_eSIg1xq_^&*xnHY;|SU^o3UUE2Ph; zaa(;a@^?^@ng<>w+E^?WVQ^#JJ*NrR%Gm-e-2XoRFR{00-c`;Uo6fxW86a&~aQ&ag zq>U;?YD=H9^R4|FaQ3X-m4AKtY!f@LRdqZ$9D2Yc>0XibufC&!%PxIYKFj&!Oh7rS z!GWvg$y04l_Rae%_{3-0f4z76{=am(Y8shoz5iO(ijqLFtoqyCDKZA$zwh3S*mBKT zMSdOAGQrcUOcvHIw-Wr2kk9q8g6ZbtpMRExiOzlPaJww`NQW&uTVQ2xVfg>LImXAC zytlFH&Taxl$E3-iVPTuHeJA~v9C~%fHut`8ld*drC~4i;p)K>`miB+A3CFfO7o^R; z{A-(}?zamqukWn*D=Xf7!l~QVr6`fbNBNm7R&dN!Z$MSu%7o75M|9q)O%nn=r|JZnV`b>Ly-ERUX zRbxL-*ut#PJNGQ>jlE$lw;QQOg_lTowR|6};$OkX2DlRaWbEwG{&0Q8eg+l` zxxGt*G~XRGWqYS29d(b(+nT{_zkJ)_>_+YLldiqa*xJ4-@56(`_UydY49g}RSQlTm z_}$LD9bRu{lr&4TivOArTJ|cw($+nRV~^~idGpH?Z}VFPeioHqR-@^DF>1>#XO;ev zvygVzQH#%OILzG}EoVwIH1`%R$w;Ph9rjY_FxCkl0ZgEb;#C^|m#IVKKk#r=D0P@oT;1e-`OK!hcrmJlemmx%tMY zzh*L(Rvn+^b)4>nq{^jU%;KMDvx94f+fnm_P48}l?T4Gw)brS#)GSL-_s94^94~Zy}ueB@gX|uLZTKPl#LK%N(+(Y+z zXBp1iJy%lCEUNXiss9`J@$Ek^??=BXQtFWI7Pb-#6<^QXJ5?`==~C|AkJ6D~TpX<@d1 z#M_`tYS-o;>Arq_@h30&XaBMrA3fjy#eC_Ch&Kn1SxbxBn@_w29#ic$oXOaIyd@zy zP3^$E!+|q$)gn}*OisOU?Wyl#mbHFzOTcPFGDFyOQ{|8P3;xQjig6GySJm3@0xSCk4SHv-zR?yCoW^z_;VwhW z(;{81^e;+k%WriyemkS5B6aTK^6J1E&u8z9viG)o>3vkZG)JJzAa!e=r!3dG(noKs zi!ax$W!H=o}Ew@ zapAtsd4U@{%NatiZ#wUG>Y*_+--6uL`;s2`e&2d1$td|A^GxYAdX-ybVs{_8|9{FI?yBzSH!ez39k+yX=b_uC~14`8HFo zYDLmSarf(ASI4#`*~MIL{!3FbLP$L+%FfLAFuzo zFnztue_mJDTf#Sf7MP2z(F#x5`aWalC*2(%jvuXfoPE}Vr!OZt4rW1ASgY_uhhOg* z%H(BQ&L$h)k!J5X5dG9)y-M`KEqld}Em6IhHGQGwy&&yLfeTiL{r-J$PhmFmvU~H` z|GRwMcXLY21GkGd;`es!jod2LF$LUHj+mo;on`akmW)}K=g<3_SaJNM`Cr2-!RPCX zl*EGSSsj-sHy*raQMB6lc6!yzg_m08O!h8&>${NO<|eOqpJecVzrWw5<3BVjxi}YI zSu6d|R~tNAuACqx+lt!6K_LL+>xE1&U_GwxLEFUSN_RNfKg^uDIA zpb2%i114%=(0OeFrzk?8wLB0jfN-ZGrvh>%nj-zo6YFTQ1Z z2UGJxjUqM6OJ%(qlfEqKd85UBYuoDY|C5xsT&f>w?A7>uasIq2jp|k0yz8^y>m|H> zrMF~F+LfBfeb4Qm*RAKMEiE*e|M1Q&ap)m2NG*7uV1`yJ^F&+*LY8mxB4H>K82dn3`+LF44ba!}pDHXQCt5lxZ$c zHDWH_H@!Xg{(Fg{&AGN;1y(u-3w{8t75a2qfBEApCSLEjtc50ixy$(Go^0FlrsiGq zo`zTpW(B_IOL%;mc}9h>dA(TZA5Hm13t1+8dh^3tl*h{}Fytee6d|7i*ja zm(Pvgul2TiD@utTWj_$0(b6toCsF>T;TLqcFn;*Nl&9din z_=(Ez?-%P;t>9q~|30V8u>%xh3p+V}fJWJu`OIumpI`Gyc>lJnt6FPfcZ+?$UvIx| z@3$z>n14H;tkkENo`$K3#mviGxF(*Pm>Xpk)vO!%WBJLS1$K+YpH^tSd0{)@=k?A~ z`JO{b&TjYreGII6z3@_H-j)3;dbc<2`K<4|VxMLBs#baDk5eK73M<97+5R}`=$~)v zGIMz+$B)nR|J%%~erH+pu6Vxn)5i+~IQVToFudFS-VU^P+WdZv@ie{IuI2OVs!kN+ zS*llUXx`s$w=Fd~aQog$pZMx!kEbm&b8#;Gp>7fK-}aNm-s8(23Cy}DcjiawciYP? zYd-e%#D}jxK98IIV!Zt1W%CyHWyj|qbMO{Qzx!%Qvn6d|9sS2Q-71C7hdMsiC7|^YroBwQDnSa$G_~ zub%8bJI7LZclmp{PcbV)xYQ;~t_*qA^S|cH=gOL2$@#6a6OY}PvG1GD@q>#EZ`=us zc=PJNoPW6thfIH5%1Zas$NzbjxmZ7v+@JUS8q2Zx=O=%Z`pd90<-U0HM(XnQQx$Kt zOq6GuEZlFY=W6RL`~H8Vgz>(q+>;~ws@#flwlO{p|M}ryzm@;bitzC7CpjY8tSWV$ z7aCc9nUic5J;6iczs7XUuNS#QHkG^#@|j~H`2WxO|2zygw&hBnv;7`ZQc|Lzug}l$ z;f(S52krKM6u0Hyue((Y+8PbnJRGpg=J3t$Dds|dInVq^+qdZX-w*Y>FKkk6`)QxJ zXJPq9gAdpBuRNUf{CN7CS)Cdi{(W}e_=Ep(V)1Et=B%UD_j8@^l~}9QUtC|Nx%}!T z)uy%ga}SrkKE5*c)C$$d2jazQJ6u;ES6uCMXsw{fEw2dj${MW`(*V`ZF`z+bi!_%;I=~4}SefQehzbE_O z-rmk%{OpY6r|cbE5tjS1*oFR{+RvzPN67u0iDftgvxOw%WTQ293%)x%V{2c0zBRCB z-Y)K*FE^Z7H`{W^$xbP~NUh~Im%jK7=`xOP5H=X(+6KeJ3Y}vwx zUk$j<7p{&yFJF%KR*f1RhH^i83eMP3gJf&2nt@D@@cX&PdciK`Q!&z(UBmetc zC+9vobne-=o|x4izuG-HAM$@W_ug*-FOSOQ{Qe%c^|Zy&YwWiTUUq-l`@URfji2zb zdMozb|5gh>p6)c=)@A18b(J0Onans=U+(e94%ztmYt0Ooi#tEvv9*wO(?3=8!)Jr# zT8^z+|7~6#-@4Dfq09AEzQMZda;MgAR2K6}++bdt&neaW@0bN!RPu~Jm6qXZ{YSUl z{+m&`UA}T3Q|vc^)IK${Cbi0^=Qk}}@3(XR4(sjh6Kw8GviOtsMb7WN(}WZm z4%u$y{As(+|9?5_<&WF3QTL}l-db?G?owZu_Fse3vRXY|ceCB@>6r7`++19y>9Sn) zEyIb>1E)Ssxh24}{AG$Eo-|$8whUdE1bZM6o9Ksnt3u zYrfgt>z(t`@5bIs|4L2LB1OY?e4g8p6a2q2?2_Ku^-OA|mELbVw6=&R+;3<85N7z# zRO!En_@Vu;r2YxBmb{ZV@TU6Hr7B*gJFH?AU2k95iMXCiDg9UPr`mX|)KRBmjnCBA z>4zfY&mX)OUzU8R^jE}#A}5_GPOqocxiEQtF_-(lTh(%gQ`l5V*3GsjR&G4JWdF># z@FlVgJGmkv#iSpc*rs%6UdwKWH~&vjRw{ z0`eWo%ia2{)vDerW=~_iTdlN;w}E&4w}aENoAT$>U6%X3#o$*p+l5$ZSHIqG^#ysi z#SX5EFHsQkI24@Ln%^Y6GhOz-wW^rm!t{C>UVV0f%>Km>Wcbr+t@+p0WzT5LBb-+tb^>`EGQ7XQ(@ zyR5vaK=kI)%D=PkRj+tuEC=e@an28X{_aYn`fr2B(!ciIjt&Xw(Dc^1@!v)KnQg?H zGRx#^JN3K1@4L-2`DBl#_a2GzT0Ncg8R>@`S<39Mb(Wf6Tb!{@TKrBtOLuJF|Mgqd zVm}MyUy$lpeWdT8-$!=O-~C4(%vLxqD|N~yeP?Rp|6JZj0` z@Zm9I{(6Q7dYc~n`OQ`McbC?Um^7PuSKC?lrzBsslm(6Hg&sI{dY#eY=M!#lJFr&= zZ`k4ULF8n;cl>#6i>o0Px7()qzO3wx_|2Xua(;R`!9l&fI%b?=BDi zy<09s=iIfw#Yv40(pFWs9?VX0ej7Zqa&lkLqlqHd?iTPod}gt(dMDT9e5t?Lzt-pp zpZo2o_{=t>@8Fuzvh&m07rz76bXN889e6eI+sVe@7nQ8bXYCQ5 zrTFRG_R|6@`lPS_-^$hKo4DuaTCvtuCTw{uzK8CJ zn(>r$@HV`acbTyCG5fa6ihD-VY)zahNAuL`oqk`d`=R#nVv~<{+zXfI-wyV~eZH|v zE+zEVgEJHDJp;TBmaqqTa$A+9emh^K^Go`1#mmj^hF`zsoQYX=$9PuQbk1=5=+4G@ z5_T0Uf;WD?bi)I*AS`U&`b&)|F>lWKemKn1b$;@5zOyr5*zODNzM9#z{6gEi=2OSo zN?$8$*u}5A`qb^k!b3OZj<7xM;{TGq>HkRuwR*>YttnnVGWh3~?&fd*AO23l;?leW zzx1yqevAFNB68B~Y1uLRb_vN^J5~IBeKyl_|2+dHwL~%N1hbg}Ts^fMZ)}Yi&R68# zS!(E1vs`TT^{#(bsp3s6uh|1i7RIcvH+JH`c}K0HYt4(`AN!vNIO!I;>=xG8?ef3d zz9G1J#UCU2TiabH1l0e3d%4R!D5Czdq0iqaJGEM~z4>W7Yf4oY{pI}g|Jyqon{%7) zN&hJSW;^*zO7^TY%_Y9kwW7jQ5=%3fK; zw}}1TXC|-ywOvz!{LB}K>{=FmN$0H9WS_LNYo7PdS>`)Cx5(c3vrXl``LAD<+WJhM zIs41>^mjY=f8X>p-rD$=)nrEHgh?&uRyjX-al|c8D9~X!>!;<~$t)9$mM1LZo45B; ztrBmk+Jr8~b6RyTPP4jCm2{6geeT4Sf75GLM5N9B(rw?Cyd`_0`Pq)euYVWaoszhA z|4)rmdK;6!oq8j4di$I5%NI6kcI@3Bc<#=LlP=r#H}y0a9}ZdJ{^G`wXA3I6wlBe_1+hEB+uP?t;MUHFv%sZS8m>~ zPTy}I&D|E1pL*|GrIOY0zG&$S9c>}Omtjoi`!Dxz;(j9fpXt@gg^w!3j?{eYai1(} zt!%LC%j92A`h(|2IWOg2^Wwdy;*EJnJ*OvI>{`cvH`;l}m+em@>vO#q^3AW?BJ%$3 z&Ay6=nZH%iq(FAoEiF7WdtzCgPgUAs#mPwyz#<+<)1oimG7zC`#+K0yWUjvmQ3wF{oNPZR!)!n z`0aXB(Pf6IlXBT72o|J8EPEyR-S^Gp^c9v|*K;b@UXtmJdbB&n(5m!%7>Rv)ONa@{`Yb$~NOu%I|HH1rJ}HP`r!LBlBXl!uHdR<^NAH zXPC|?ul}W8lX9r_{M6+8-Ev>PKYNhlB$>?OGsz%=aV;Z9^`-8{=C_ih9 zr`YG}!q?6p_cA=S`tP61l^#98N7Z%EZJTAB%sgg1A%asskw$*IYl$fhtEtO||^IcUivwPtirqZMP z51cR!`Ms3s%iH(vTJFuUv)|_(o@u@-L3i5q7r*z^>p3W~nd~{K@x#LHfQn@CM7f9G zzI2@|S7G~JlbXX&J7?CT-7y7Wk_iTxEn$6bRg*Z%jy+sA_4%Ektyd?kz0mM%Kf6l( zy=7bX7RWT&PV_GOKS<|emtK)s#n-v_nmnK8g6RK#rC>i z+I^CJ)}j-$7p@U>Fg;qnGii2~!~Kb!?rn4BT>`&J{`Aucjpm(Ncd9Q%NAcwe?;X~S z9p-J9K3^&Pa(kuQ>V+@=2Qd79m=Sh9@!roZ?3Gikrp!^kf8UzH;{C*=E%{uEj~*H* zC?`znxww{H=y5#Dqz>JF6`#+ayQe=Fvb?kUS+^k|n)hW- z!?iVz*WEA8Jyr1IYn;V71I=eFiYt4=oR`_?_T(|&GW)&x(zFPJ&SUn!HvIUg8Chv7 zShbn$a(+qDT}trx0Qx@OpHtCpFTqh9K|*l3#||C(j3{$=mn`7DW=lfIvrAeFu1 zSC+NZJ&s3bL?6xFR#46qSTCeHe^U9K8J@{k1D;-GwfL}mvB$BA7bjaBisdZXCx7JX z>($2Jm4q&5H17N#)bRHH`vaT)&zbi6S&zmIx#Rn$@A*(C;;J;M;o`#9od+6L*f-Y2 z%Sf#i+omO)wzlMut@dC3Ju_eNpWMcrVIp;d@r*+7KE}pt95Ken-`&`n->1Up4)?WB#?2H>O4le`%?I+xNiv`#l4$^wZh9SL|DC znDaYPUFq##4a@tDo8N1@Je2vIDO(wz%(J$hKk4<~+zTpOW|WsRp1S;F`?-*MkE8pi zy-L%)u6)O8-wVAvQmw1QB0vAv@HS~UxYT`t<=<9?dwT_a_6F=#b}mXicin0A^@HjM z<@#l3AD&Uk{P_OH5dH0oey*Q+n(MD`VM(pPfj#d<(%zUFt$f$$l(X00u6tvU=YwYX z^2cwdo}E#(kTrbm?7r63T>&=#znqgd*pPbd&H5*gS6yFi@I9XE!^OvOf1b+E7I=B` zEtB}oGKQU87Mpv@{_fb`)Y&@e#C+xT{&9~NZWEoeGpXps#o5oE9Xx3vz<>O>Dbrl< zthI{iy@K=lkBT*Y`8|nat9xL`@^^m!Dza7-Roc~j4O}8LGxPB8{;~&WL|G+!XBe-{ znJ607EB5wwk^RDT`^`JIhMOJ;mp}T+K0NtbcueNsnB}og4qUBY=W?==Z%<_?bJVfI ze`eu~lb3v2Wy(?BU%Q>NdXmSHNd`|fwapl=Y;Q6S_d0KUpj>0`?BMV3I7-ykzABJmUF_vS2-s;)+gpZVfH<)#ACd# ztbgS$I;FpD;nUcc+wNPwYd*fq;H%Q?N7j#m&u>?}GSTTe-^+Tp#}|M33C!MiyHC~M z|Dfb`FXOcbBXf^jeZA_!9=(RM{Au=?6Cl;>)6f}$-+N^ltG~M?%w5!@%j~!G+7G?k zZ-RL{zNuZlb!J<2I!mzrgnIUXwWcZhmz|ei^r`3z(wX73z_M#n>l=o??QN^p7cR0#QI&-LiAeEH*3*R@sRFXnB?+*p&Vy6b*_ z&kfEL7N1F*M6)M$^w;)fzu#xkC9d)LzT?`DvD-4gHe6fa_~6Em{D8l*oFdn58r{&! zT(PrFCz$ueP3w&wpKi3--8g-<^4t5B$D*CT8QJQ%vVB=*IXB5f>*Zz?` z`Q`N|!Hlkj-wVF6wOeg#+wgK8U(AQkK{uAAq)CLXW}mh{vU6Q-;+u*?)wk{&GpU_? z*&|`h)>q<|xPIDy58kKKmR@1_zwqN(N$xF27jDf;Joo&q`#G6VHksQ-8Inp z$?IqNho@@)d4#3zfARbj)A{0n&6&#bU8%$s=TR}rI){1@d%XD>{8`G0B;bN%k6lIizjJ1(`} zj6Au{JiFkhx9E=_>9HT)*ne!CC1~JzF5$|->lU9l7OXd9>kZ`IYxKK*Rp0wfXJeMQ zKkxh4E4DSadtpkWkz4zf80UA9X)ZRMhRZwWtg3bAs9Mdo%;x64y(ytqrcYLxI<3BV z=(@tcwEf>*?Vr_V7oGWQJ@eig$Cs{UZ~kwe*R$_j^mdblUu>k$?XlbaJHkqf$Kw4% zp`1Nyc|U$Q@a%j~gJ%)zXC96%eSb7hl^ZZ^b~@Lyf0oi{CvyMQyH}*&F`G&jN<<{ zTifivF5UV;M|O8u^~Lwz5C0r5`us5eRw|#*-ZNd>_g|ZI@%lt?4tW;I#rwg$_s_?D z>GeE)nWyeDxt;wOyX~p(Wo0?9zRC}kSLQpdG3QL5YUK9)YMgUl7ITgJrP=H^_FDIB z;+|r7Xg2T5YXLV?!%A|*Zof)-zrHG<`t_5yEz#Yj1%BoW3uEU$IWKE7?ZCqORfR9} zPpj6yE77Uiz!DhGR9$!XG{=l3>DuNB59KRaHoOlIX-(Jse&6Bj$5_1?ci6tjN&M|y z{duy2j_WDOsR7lzzLw`(rfolFxAMV`d(MUmO9P+Wd^mGb9`97$EerRnw?z~rrSJbW z`%%ZYBY|IUo~&iLKgT-u#`ay|pPtvgzLs|8f1IJtyF1I{#sA(fdK3S+(K9)b#phDd z6lp=;o;!>^CC4Y4J1-Pw(0u)AnY3VzT(A<*oDX}qV$-c1&u$CkN&i3LeEI%g z4cqw2<8G%q!}udk@nqDNCmDWN7hhrcG>P40@}IkVpX6{)x+t>ex7xO#S_}DqJGieU z*)EIM?Y|g&Tf{bNe++|jZKm72_-(3(I4nMwxE;v2sC1$1SK9i{7Y8EOR*IgITx&V? z*5lndr&nJ)7;XP+iOlHdY{0#U$Zv#*GjLy6)p4Q$@IIEY{C~<_AT3L7xTeqrpkHssh07@ zic@_Qud4lhKQ+g?p{HT9^ZN^mH)?*VI&)J&Ra_mOof(Tz2H5|1w}FZ}e|e(9P`Q>6uSY->60 z)>|-bp8eGE@$~vO$%_w8xVL|m>2ONFcpziG%yOX>t_O3JZ*A+@EADfH*vH9f$Kfa_8LHvX0+yzF6_TS^9ye zRnj_?&3u0;e@OYwcb>z-nnzF9>F)CSs#lf*?gx2N>nDX**G$i?GHy6>Lc7t{nrVYy z@U<0o*&G_S*V*nb+~~e|*0olHuSrjP-ld&($v$Cwe9kqk>-&yLBrRlFrn7PR-zYD$ z&)TdOe_kC4XI?((K)e0Nq>KICDmhcP?BJSVq_jV%eP5!P>{t2TgQosB<)2F4Pn~{l zMdPsrt>xc8*dF*MTj&$HVvVco8byWu$65-G-;j1{yM_1RfP=cAHZ-e35z`jeBapvc_S`^+_#UY2C^-dJ~qccE?dB9Ba=O*gEP zHx-;UQk3*kx@7!mOV{Hd&jr613#8S%%68lPnZ6BT|F-9W@9*D92UVp0Uupj=clwN0 z^v=oP!lcaNlrBdp|KBGf{R|sDm?|$%v1K~(L+>|N+Utv(|9_c!V55M*cU95t_A=QS z4aIGe_iUB(w_p5IZXZ+p#JXvcxA0oUUa<|KvoCyoc;mN$+5WffhW_?1&;6BMe@ox; z;%4peSlX-HtsBHQQA%t7jh^K9y?dE{A2@9{ z=O^peWAC3a$@&~++pTqS(Xy>yI>Px-b`toUGuR8Z=d}Y{B!% z7tbg4{hG`n-@W*g*UWVZXD3DepVH2(-}^sh!dv@SN4A{Ky&!pYw{{~WW4gb-pj7ha z$l1G_jrqdYe#(E=^M2F$1lAx=UZ=Uv&vh%meOVmhDD5DBTuLfo{v!7+`mGzDfBBhl z-0XLvrdh%7>4F_sCLCSMToP6Bt#{9XhZp|8Zv5+`TypB}=6IoC{d3L_B6$8DK9gx5 zxtgwjFQF@a>pK%(wZ7_=Kduzzd@f2k>b~61h#k*L70>+Iw9EZ}L!J8=dLV#|DsdyKtM82pVm5tWu1i~5*MQD z7cooZ_$w^+Kltd;Tep)r!MquFORCs6y!Vhu71ydhA7H7y`v%+lcMshenAaA%G<#0b zdvac$A?@z8{%z%M_Rf}%d>4c}F0P-ca@*QS`N5mQkXaMo$V-^@w>DHLG5!{*T;CBeoO^$_m%02oTgV`FbnUxyY8R!n{%?^z zIPJ!U-v;mMZ?-Qgo?^aLDY)zTEWuYk3I#b&o*t05ed&6VHSB25KgOO0&E_>ntA6zJ ziG-g&Xy9}1Um+Zw|@n$IPgO}fp_j1)5}c}2`{8yZ!G`yU`5!^ zCz}mBu1z?{5#e3>;#T3WGK078Sq&b~k5M{1yJT71qXP@W+v0hI_?KVUsM^5#T`YZP zXp#hD@G*;1ycaUq-58j81!tt}Qx`fJ?byruaGkKML%p+zYx?HN3rgH}vsq0RtMyg? zmm7VLcX8ItPvSi%V{fLf$@q2p!jHYn?YMlo3eG-wf9>R|^(UQQAA7jH`Q}$|{ZG$l zC2ZOJn)$>I<`Wzi)my9_7GGAlG12Av@xmSHetkddd)^lvPGB|hf-X7u!|dm4+6xUwbI!<~R5|`$F{Mtgpg1`IO4d45nn{G1`lBG`v<~lsT$Fk_)$&E}B2{qHbzj^6a zy;^LsZKph|@l@+YXQZS~Rkm?9F8=Frtnu7;(gp%(7C)_9tt{<{W((&Ytw(?)l1sV_Z-3rmg+KapU#ZNB@4mcaM#geYfK=pNYA7 z`s<~cGCi$TUu51qmpq`-*(AlCRBQN!J*nV)P#*gPzvuJPolZyZHL&7cWnZT=Tjpui zl3&3Yb>>HYl*?r06)bumAo6#g3KMVOnH_;UWe@nA^kO%d>hB^yl;WE%JfOq<72%qUc8X` zQx(@88yzFnve(dJ)t5~f{UXx`K`MZ|gWMV}Yd>aR=V zwzB?H{yU?r(^u^6#NVb32RS0XSgu~bVbSZq0&~0Lsh)f<~qxvQ$#r=_S1apv6I#A<44db6SubV9t9XP6tq z^YiojU$5KE=WqK}By3GYBd@gCfee$RmzS0**ZG@BJ?2p01f97U;jGeS>axh_s{E-_ z91Fr$AANOowP2@9!RNE)ppBI4_kNSonKp5|-~+b{DqXXGOwza}(3!+LafQ#<_D$QI zr{8{i40Mv2`n(FIXnjx6vC)SoycJv9)8HwxSRCxVHHt!0SN@xjW8DzL>a%FauRfQJ z5T9sPi5dtB3JRWV&v9pvWSp!4IzrLZWf90lf(C)JB=$~#TBhN40A#WX=gmF_BPO-L zo;Yw&fXq=&SY!e^H4_?^Jq?|pbCZQO$ug#}_$)#-A`^6gG01dN(D4_-pu^`?QmmOy zaLjN)HR2NZ>`aI&4VlzHLF##;jAeqLK_{vaOTi~m@=gT1lLd5|w9BLkZ#fi{6NHX} zj=F>!;U$a+jiVehIyyQ!R1|LuI=CHB=>i=<2{&Sj8cHaD!pGBTn{oqcfF9+Tp#eU> zG`j;Fn#u_vH%@994ah+kkdR^K8@F$t=2PIjxuZ}yXr&0~nClqxRuNFfi*Q!a3jV8j zvy8!3PT#sgCFSj{|Brg7w;-~+%5A;{e#=wYC$!vnYj#b5GZEs)Qz7*otHH@VL`_-t zfcd=&=IZb7M63G17iDM$c6;B!BuW;fn$i!pB&>4D@-kXujA+qi$f z{QTc{o*NWD^O3PANH}qSX^`ffy5Da(TAeEX{d^u66x8&3{eHQL9xkgO7s(V}TB0&5 za)yuE#^UFGM;;gIL~Y>!-GCAEbxZMczbjX-9(?@qz|$f{x&GvLcXn3%dbxaK&dp7U z^Ve_LA|lt%{`Bcn50y#Hd+$HK{`2S04N8;n16&{lN8W73d^6US3{hyZhz)tHtM4 zy;|J@hjPzN*Z0TozN>eQbA9W&yLsYC8)xL)I-LD9_WASY&g-x9+f94@ z?L$=M^}OxUpe^A`%og$`Xiaqr4!-=LUR+f4VAj^YA7ADA*;`v%o#jrSKJD!8e!Sy; z@%z2sAKa^cKhxtgsJO@w5uU{?(d!l%8fy9@{ro)J#dWogj*NTa-cLCGr_TOP-glR{ zc_4)>J}a0rjApj+q|dK@w{vk_ZE2|~=rq22_u?2FlqMdSrW?)n{3__+>%ea*PgN%T z{QUg!uGi}h2m9N$ifV@?9B$(^km`-tRg!52G0<<))ZTAfqIAR7$Mt%sG|kT6C;9d3 z*NV4WuP1Db*!1xesAw(`I=QV8bkMa@>VrUOkcie{KJax&r#3KqxVX41F?s{OJZU8l z_!MOqtkS%X7x&v8`t|j7@m-!EgKGjOXRfhsczu2S@vYb6y1n)HisIbx%DUMu zT$`%D=Vh48QmoJ0`E**0dAEpy1ka;%``?n^-rO{P@-afEO+2oG(R-TCL0D_^~ET7Q4dt=gB@_x~&Hi}EbGdg)Tr?)UqoOG`_4e7lvsG4JlKgZn4HvYRq> zbAx))>uYPlXJDL;|0fi-TD1Ibskn|9x847r=b!J&Tq1D1Pu4j%S9f0ZyPXG*fBN*P z=!Q~I6KGNK>hSe$-rn5L&d;}h@;x$GZHTi?6S5UzB%| zNffWFZIuYGqo?ow4|2W8BH=%pk zB%j_8XNG|I`2KGFeKWqRN?mN1uj^Rs-v8)fyL{WtoMVrUc5f_w9rovKe!X*S?A(;^ ziZ|U3Ja@B7PzkK`l zY*(a8@z$+dBQ_*3c6E2R+uXZ%uW#Bkv5@d^emk@Ldp2)wZPmW*t1?r1WgA1n?{9A_ zKA$x&JZJg*L6`Qr1K#?3TPkekHHS@#IXlagTiPtA;bsmSL(Abr8OtJ;;AK7&zu2jq z=1b5JadLKMUK95EL9$pkE5i>qzsKA1@7wK&(YsUs|1apmfeodv!)kutzW*#HXx4ac~YR%vbTRX%EISHETd{`wCa3oLA&Ot!d~F~?6`sFP*Ot;UVh zKcBakKbRnJ-tMyg0!cdm1=QXc<^|SoAMk#b?bT)4H2$~-toVL z-jrE`u5_B6w~Mpt?}3WBQ@$%k_AxLW04=O7dg^s%rZKx%H>dGHrd~`VTm@+MWUL_moyjA6X znK`<+7hH=>_bijYv?OTd zl#pI<`|{`Hr!wGMrzWz4E}!#Ay}<~&A}C0bH0hu6W^9_wzmKP literal 0 HcmV?d00001 diff --git a/doc/notes.typ b/doc/notes.typ index fa0b780..ecc9b8f 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -335,7 +335,7 @@ composition breaks this: 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 & | @@ -353,7 +353,7 @@ needed. This is done through the following composition: 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 & | @@ -620,11 +620,11 @@ simulation loop as follows: models has much in common with code generation for synchronous languages. ]) -In [Branicky, 1995b], dynamical systems are defined as follows: +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 + 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$, @@ -694,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. 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? #breakl(```ocaml diff --git a/doc/pres.typ b/doc/pres.typ index 88c12a8..f4b52bd 100644 --- a/doc/pres.typ +++ b/doc/pres.typ @@ -145,7 +145,7 @@ Solution: simulate assertions with their own solver. Nodes take an additional reset parameter `'p` for their reset function: #align(top)[ - #grid(columns: (3fr, 2fr), align: (left, left), + #grid(columns: (3fr, 2fr), align: (left, left), ```ocaml type ('p, 'a, 'b) dnode = DNode : { @@ -200,7 +200,7 @@ Hybrid nodes are a bit more complex: #slide(repeat: 3, 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: $(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 => [ #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: $(h: bb(R)_+) times (u: [0, h] -> bb(V)) -->^D (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`) // (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)) // sim -> bot @@ -583,7 +583,7 @@ When receiving a new value, store it // solver.step box((6.5, 1), (3, 1.5), content: `step`) - + // solver.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((5.5, 2), $"Signal"(beta)$) content((14, 1.5), $"Signal"(gamma)$) - + content((-4, -1), $alpha:$) content((-4, -2.5), $beta:$) 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. #linebreak() This calls for a "lazy" composition: we request rather than provide data. - + #align(center, cetz-canvas({ import cetz.draw: * 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((5.5, 2), $"Signal"(beta)$) content((14, 1.5), $"Signal"(gamma)$) - + content((-4, -1), $alpha:$) content((-4, -2.5), $beta:$) content((-4, -4), $gamma:$) diff --git a/doc/rep.typ b/doc/rep.typ new file mode 100644 index 0000000..fe4061f --- /dev/null +++ b/doc/rep.typ @@ -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] +) + +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) + += 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) diff --git a/exm/ball.ml b/exm/ball.ml index dcdd63c..72e804e 100644 --- a/exm/ball.ml +++ b/exm/ball.ml @@ -40,7 +40,7 @@ let step ({ zin; lx; _ } as s) zfalse = 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 } -let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode = +let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec = let yd = cmake csize in let zout = cmake zsize 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 state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in let reset _ _ = state in - HNode { state; fder; fzer; fout; step; reset; horizon; - jump; cset; cget; zset; csize; zsize } + { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; + csize; zsize } let errmsg = "Too many arguments for the model (needed: 0)" let init = function - | [] -> bouncing_ball () + | [] -> HNode (bouncing_ball ()) | _ -> raise (Invalid_argument errmsg) diff --git a/exm/dune b/exm/dune index 9eb098f..4eecb1e 100644 --- a/exm/dune +++ b/exm/dune @@ -1,3 +1,5 @@ (library - (name examples) - (libraries hsim solvers)) + (name examples) + (libraries hsim solvers)) + +(include_subdirs unqualified) diff --git a/exm/sin1x.ml b/exm/sin1x.ml new file mode 100644 index 0000000..a6f2447 --- /dev/null +++ b/exm/sin1x.ml @@ -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) diff --git a/exm/sin1x_der.ml b/exm/sin1x_der.ml new file mode 100644 index 0000000..e4bf116 --- /dev/null +++ b/exm/sin1x_der.ml @@ -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) diff --git a/exm/sincos.ml b/exm/sincos.ml index 17131c6..ab62149 100644 --- a/exm/sincos.ml +++ b/exm/sincos.ml @@ -34,11 +34,11 @@ let sinus_cosinus theta0 omega = HNode { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset; csize; zsize } -let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)" -let errmsg_few = "Too few arguments to model (needed: 2 floats)" -let errmsg_many = "Too many 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: [float, float])" +let errmsg_many = "Too many arguments to model (needed: [float, float])" let init = function - | [t0; om] -> + | [om; t0] -> let t0, om = try float_of_string t0, float_of_string om with Failure _ -> raise (Invalid_argument errmsg_invalid) in sinus_cosinus t0 om diff --git a/exm/zelus/ballz/ballz.ml b/exm/zelus/ballz/ballz.ml new file mode 100644 index 0000000..08df6f6 --- /dev/null +++ b/exm/zelus/ballz/ballz.ml @@ -0,0 +1,131 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let (+=) r v = r := !r + v + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ball = + { mutable major : bool; + mutable h : float; + mutable init : bool; + mutable z : (bool, float) zerocrossing; + mutable y' : float continuous; + mutable y : float continuous } + +let ball cstate = + + let ball_alloc _ = + cstate.cmax <- cstate.cmax + 2; + cstate.zmax <- cstate.zmax + 1; + { major = false; + h = 42.; + init = false; + z = { zin = false; zout = 1. }; + y' = { pos = 42.; der = 0. }; + y = { pos = 42.; der = 0. } + } in + + let ball_step self (_, ()) = + let cidx = cstate.cindex in let cpos = ref cidx in + let zidx = cstate.zindex in let zpos = ref zidx in + cstate.cindex <- cstate.cindex + 2; + cstate.zindex <- cstate.zindex + 1; + self.major <- cstate.major; + if cstate.major then + for i = cidx to 1 do Zls.set cstate.dvec i 0. done + else begin + self.y'.pos <- Zls.get cstate.cvec !cpos; cpos += 1; + self.y.pos <- Zls.get cstate.cvec !cpos; cpos += 1 + end; + let res0, res1 = + let encore = ref false in + if self.init then self.y'.pos <- y'0; + let last_y' = self.y'.pos in + if self.z.zin then begin + encore := true; self.y'.pos <- -0.8 *. last_y' + end; + self.h <- if !encore then 0. else infinity; + if self.init then self.y.pos <- y0; + self.init <- false; + self.z.zout <- -. self.y.pos; + self.y'.der <- -. g; + self.y.der <- self.y'.pos; + self.y.pos, self.y.der + in + cstate.horizon <- min cstate.horizon self.h; + cpos := cidx; + if cstate.major then begin + Zls.set cstate.cvec !cpos self.y'.pos; cpos += 1; + Zls.set cstate.cvec !cpos self.y.pos; cpos += 1; + self.z.zin <- false + end else begin + self.z.zin <- Zls.get_zin cstate.zinvec !zpos; zpos += 1; + zpos := zidx; + Zls.set cstate.zoutvec !zpos self.z.zout; zpos += 1; + Zls.set cstate.dvec !cpos self.y'.der; cpos += 1; + Zls.set cstate.dvec !cpos self.y.der; cpos += 1 + end; + Bigarray.(Array1.of_array Float64 c_layout [| res0; res1 |]) in + + let ball_reset self = self.init <- true in + + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } + +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable i_73 : 'f ; + mutable major_62 : 'e ; + mutable h_72 : 'd ; + mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } + +let main (cstate_80:Ztypes.cstate) = + let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball + cstate_80 in + let main_alloc _ = + cstate_80.cmax <- (+) cstate_80.cmax 1; + { major_62 = false ; + h_72 = 42. ; + i_70 = (false:bool) ; + h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; + i_73 = i_73_alloc () (* continuous *) } in + let main_step self ((time_61:float) , ()) = + ((let (cindex_81:int) = cstate_80.cindex in + let cpos_83 = ref (cindex_81:int) in + cstate_80.cindex <- (+) cstate_80.cindex 1 ; + self.major_62 <- cstate_80.major ; + (if cstate_80.major then + for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done + else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; + cpos_83 := (+) !cpos_83 1))) ; + (let (result_85) = + let h_71 = ref (infinity:float) in + (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; + (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in + self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; + h_71 := min !h_71 self.h_68 ; + self.h_72 <- !h_71 ; + self.i_70 <- false ; + self.t_63.der <- 1. ; + (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in + (begin match z_69 with + | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 + | _ -> () end) ; + Bigarray.(Array1.create Float64 c_layout 0))) in + cstate_80.horizon <- min cstate_80.horizon self.h_72 ; + cpos_83 := cindex_81 ; + (if cstate_80.major then + (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; + cpos_83 := (+) !cpos_83 1))) + else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; + cpos_83 := (+) !cpos_83 1)))) ; result_85))) in + let main_reset self = + ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): + unit) in + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak b/exm/zelus/ballz/ballz.ml.bak new file mode 100644 index 0000000..d248fc0 --- /dev/null +++ b/exm/zelus/ballz/ballz.ml.bak @@ -0,0 +1,137 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = + { mutable major_50 : 'g ; + mutable h_60 : 'f ; + mutable h_58 : 'e ; + mutable i_56 : 'd ; + mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } + +let ball (cstate_74:Ztypes.cstate) = + + let ball_alloc _ = + cstate_74.cmax <- (+) cstate_74.cmax 2 ; + cstate_74.zmax <- (+) cstate_74.zmax 1; + { major_50 = false ; + h_60 = 42. ; + h_58 = (42.:float) ; + i_56 = (false:bool) ; + x_55 = { zin = false; zout = 1. } ; + y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in + let ball_step self ((_time_49:float) , ((y0_48:float) , (y'0_47:float))) = + Printf.printf "STEP (%d)\n" cstate_74.cindex; + ((let (cindex_75:int) = cstate_74.cindex in + let cpos_77 = ref (cindex_75:int) in + let (zindex_76:int) = cstate_74.zindex in + let zpos_78 = ref (zindex_76:int) in + cstate_74.cindex <- (+) cstate_74.cindex 2 ; + cstate_74.zindex <- (+) cstate_74.zindex 1 ; + self.major_50 <- cstate_74.major ; + (if cstate_74.major then + for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done + else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1) ; + (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1))) ; + (let (result_79:float) = + let h_59 = ref (infinity:float) in + let encore_57 = ref (false:bool) in + (if self.i_56 then self.y'_52.pos <- y'0_47) ; + (let (l_54:float) = self.y'_52.pos in + (begin match self.x_55.zin with + | true -> + encore_57 := true ; + self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end) + ; + self.h_58 <- (if !encore_57 then 0. else infinity) ; + h_59 := min !h_59 self.h_58 ; + self.h_60 <- !h_59 ; + (if self.i_56 then self.y_51.pos <- y0_48) ; + self.i_56 <- false ; + self.x_55.zout <- (~-.) self.y_51.pos ; + self.y'_52.der <- (~-.) g ; + self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in + cstate_74.horizon <- min cstate_74.horizon self.h_60 ; + cpos_77 := cindex_75 ; + (if cstate_74.major then + (((Printf.printf "idx: %d\n" !cpos_77; + Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; + cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) + else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; + zpos_78 := (+) !zpos_78 1)) ; + zpos_78 := zindex_76 ; + ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; + zpos_78 := (+) !zpos_78 1)) ; + ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; + cpos_77 := (+) !cpos_77 1)))) ; result_79)):float) in + let ball_reset self = + (self.i_56 <- true:unit) in + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable main_ball : 'f ; + mutable main_major : 'e ; + mutable h_72 : 'd; mutable i_70 : 'c; mutable h : 'b; + mutable t_63 : 'a } + +let main cs = + let Node + { alloc = ball_alloc; + step = ball_step; + reset = ball_reset } = ball cs in + + let main_alloc _ = + cs.cmax <- cs.cmax + 1; + { main_major = false; + h_72 = 42.0; i_70 = false; h = 42.0; + t_63 = { pos = 42.; der = 0. }; + main_ball = ball_alloc () } in + + let main_step self (time, ()) = + let cindex = cs.cindex in + let cpos = ref cindex in + Printf.printf "main:cindex: %d\n" cs.cindex; + cs.cindex <- cs.cindex + 1; + self.main_major <- cs.major; + if cs.major then for i = cindex to 0 do Zls.set cs.dvec i 0. done + else begin self.t_63.pos <- Zls.get cs.cvec !cpos; cpos := !cpos + 1 end; + let result = + if self.i_70 then self.h <- time; + let z = self.main_major && (time >= self.h) in + if z then self.h <- self.h +. 0.01; + self.h_72 <- min infinity self.h; + self.i_70 <- false; + self.t_63.der <- 1.; + let y_64 = ball_step self.main_ball (time, (y0, y'0)) in + if z then begin + print_float self.t_63.pos; + print_string "\t"; + print_float y_64; + print_newline () + end; + Bigarray.(Array1.create Float64 c_layout 0) in + cs.horizon <- min cs.horizon self.h_72; + cpos := cindex; + if cs.major then Zls.set cs.cvec !cpos self.t_63.pos + else Zls.set cs.dvec !cpos self.t_63.der; + result in + + let main_reset self = + self.i_70 <- true; + self.t_63.pos <- 0.; + ball_reset self.main_ball in + + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak2 b/exm/zelus/ballz/ballz.ml.bak2 new file mode 100644 index 0000000..05842bb --- /dev/null +++ b/exm/zelus/ballz/ballz.ml.bak2 @@ -0,0 +1,135 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +let g = 9.81 + +let y0 = 50. + +let y'0 = 0. + +type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = + { mutable major_50 : 'g ; + mutable h_60 : 'f ; + mutable h_58 : 'e ; + mutable i_56 : 'd ; + mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } + +let ball (cstate_74:Ztypes.cstate) = + + let ball_alloc _ = + cstate_74.cmax <- (+) cstate_74.cmax 2 ; + cstate_74.zmax <- (+) cstate_74.zmax 1; + { major_50 = false ; + h_60 = 42. ; + h_58 = (42.:float) ; + i_56 = (false:bool) ; + x_55 = { zin = false; zout = 1. } ; + y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in + + let ball_step self ((_time_49:float) , ()) = + ((let (cindex_75:int) = cstate_74.cindex in + let cpos_77 = ref (cindex_75:int) in + let (zindex_76:int) = cstate_74.zindex in + let zpos_78 = ref (zindex_76:int) in + cstate_74.cindex <- (+) cstate_74.cindex 2 ; + cstate_74.zindex <- (+) cstate_74.zindex 1 ; + self.major_50 <- cstate_74.major ; + (if cstate_74.major then + for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done + else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1) ; + (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; + cpos_77 := (+) !cpos_77 1))) ; + (let (result_79:float) = + let h_59 = ref (infinity:float) in + let encore_57 = ref (false:bool) in + (if self.i_56 then self.y'_52.pos <- y'0) ; + (let (l_54:float) = self.y'_52.pos in + begin match self.x_55.zin with + | true -> + encore_57 := true ; + self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end; + self.h_58 <- (if !encore_57 then 0. else infinity) ; + h_59 := min !h_59 self.h_58 ; + self.h_60 <- !h_59 ; + (if self.i_56 then self.y_51.pos <- y0) ; + self.i_56 <- false ; + self.x_55.zout <- (~-.) self.y_51.pos ; + self.y'_52.der <- (~-.) g ; + self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in + cstate_74.horizon <- min cstate_74.horizon self.h_60 ; + cpos_77 := cindex_75 ; + (if cstate_74.major then + (((Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; + cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) + else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; + zpos_78 := (+) !zpos_78 1)) ; + zpos_78 := zindex_76 ; + ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; + zpos_78 := (+) !zpos_78 1)) ; + ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; + cpos_77 := (+) !cpos_77 1) ; + (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; + cpos_77 := (+) !cpos_77 1)))) ; + Bigarray.(Array1.of_array Float64 c_layout [| result_79 |])))) in + + let ball_reset self = + (self.i_56 <- true:unit) in + + Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } + +type ('f , 'e , 'd , 'c , 'b , 'a) _main = + { mutable i_73 : 'f ; + mutable major_62 : 'e ; + mutable h_72 : 'd ; + mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } + +let main (cstate_80:Ztypes.cstate) = + let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball + cstate_80 in + let main_alloc _ = + cstate_80.cmax <- (+) cstate_80.cmax 1; + { major_62 = false ; + h_72 = 42. ; + i_70 = (false:bool) ; + h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; + i_73 = i_73_alloc () (* continuous *) } in + let main_step self ((time_61:float) , ()) = + ((let (cindex_81:int) = cstate_80.cindex in + let cpos_83 = ref (cindex_81:int) in + cstate_80.cindex <- (+) cstate_80.cindex 1 ; + self.major_62 <- cstate_80.major ; + (if cstate_80.major then + for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done + else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; + cpos_83 := (+) !cpos_83 1))) ; + (let (result_85) = + let h_71 = ref (infinity:float) in + (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; + (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in + self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; + h_71 := min !h_71 self.h_68 ; + self.h_72 <- !h_71 ; + self.i_70 <- false ; + self.t_63.der <- 1. ; + (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in + (begin match z_69 with + | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 + | _ -> () end) ; + Bigarray.(Array1.create Float64 c_layout 0))) in + cstate_80.horizon <- min cstate_80.horizon self.h_72 ; + cpos_83 := cindex_81 ; + (if cstate_80.major then + (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; + cpos_83 := (+) !cpos_83 1))) + else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; + cpos_83 := (+) !cpos_83 1)))) ; result_85))) in + let main_reset self = + ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): + unit) in + Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.zci b/exm/zelus/ballz/ballz.zci new file mode 100644 index 0000000000000000000000000000000000000000..08fa125c8a8c6a515d7c97cca6557fe41356e53b GIT binary patch literal 6309 zcmZpfx@;cI57NoU}WTVT(DsZ3nwEZpQ8gKBfkS9Bah>TDJ*ZaA5e)km2Cq;N-yYA0d+N;Naj4GRa}V#08UBSQkuU;dGe9 z!ZBfk16TyCb%GPae`g1fN{0UnCb6(MOo9kIF#LA`34=^r_Q1(u!G*~SE=*bAqf}|I zpt^R$6c%F-M~G4f$0`S?rH+*gSXkg32A2gatO%YHOc8^#14x;J1ITO^4!DXMxB^EQ zs}9Csh=#L#;VeHm%OB1PfV1G{GQ`5g;v5#Rh=E)Mw%lO>3kS%L4F4S#uy8v#I82zZ zz)i6V%#~ZfA`F)kgR>+dESQ1yFa|>mj8y?+RD*r#2#OPgQz93zusbYZQG+RBh+4qH z2IoY;Ss`#%AeU2FIs^BSSD;Q|JN~6^8{ZtP5BW z9>eEZ77kEgEMSp`o9zi@fh}IZq5}y#C*?|Y190g1E?_Z%Dp|ndzko#qq+sF#SEWjI zhW`s#j25t1IV@m_U%(OtQMG`@$YB9X9E9VrfF%*6egX?O*swH+@B)_H1uO{*SfUrO zApD;K^HeR2!Qi@prF;PkDD>GN8XXp}aKdHN7O-$REMTc{aC8J`4N#oBC{}_kLYR@f zfJGD$X^Ef|12-W7g_nfFOGV+ugTt5;6nIdVK*Ng_nG0brV1c<7ss9v|$R%Q%}c9EHwyqYaAA^ zoCd3NSWv~_3@(RQ_AOxHb?|aK{|*;Ge1`=W4zfsI*tWpI;ljd+u;NjNrP6h7(xg@(e=y zIQFGaf+D zjVxCuIIvuT#1+d$2bS}Uj544|Vmar)a)yyn8Wb_I4F6e9vz%IRVFD{7;{ry8EsP9X zH%wt=V(tXS*v2xI#ehX0;`wb1{~cHvSsgAourkUmn0R4=gM)*k z1B(KPudra^1r}CN)y2xFc)@{X$%2U&96$;jSe7|BIyo<3aARcH2J*Y!0tPn*cNX1A ztc;~B+AORrdMu1A+8{H0Kzcy-vmCtOz{w*V4R z6C4@-JAr)R0E#{)hXs%}>I6`To?v7+#Uck$=itbY>408KUI*Fb3=0?;o-i^zXJzVu z8L~p4SE>;X@mBCq+a8?Cs&_b*Ol}_;T-vK3PAn?bE%u!|S4 z2rXdYhnrXqXH~;lRoHETCFBxlX$WrFD#P6XZR0vDV9|s`f(Nus-T*V0kzv{b7DJ>+ zn1UM2lTo=-QMofniv))SEK?B{H!NV8ie&M8u+c1YppFN18(A_yMIxwVU%)aA)Ww6y zEnryzV}M%bpbj}GIKe>)D-IYLZaXYsWtM^{1eITx85y2~Dgh={hYJoY#}-Vyz{;oq zssuJnVP&?Dc64xDz{qf(k+J*2`U4K&igQ6Ai`fE3hL4O4pI8}vCb1+S41a*+@H?QS zhVa&Vus@&?!^rR#mHQi&`v=5@x)JPtC>zp(KoWt~)~t-W6C7BX99fxjSeZ;%nT!`) zm;q^nGgm298h}gc-;4}@S()=Du~Z=(`~=Ct4-pOqjZ-i(+ymJP4?+-!;Xl;SZcyhF zZ0dhTMkZF~`bjMC@|%(28Iq}w&`rG$G6>C7P(=J#{wfWIQp7F(H8)!Zdn;i9UL4$(Kj~+ z5`Aw`qA!T${{luvUPeYfRwmU+tc+Zc;6SctUPImJ$jI;l>KA0`R|x6nh#UqIWcYu9 zgMoovaDfA;TB;UV=l-y&#eM-JqXZ+P3@h^c425te69o#|>793a*%*x2(;IM#^QJImk-N9kO zg=LWDJ&WN2Mn-K$MjaNhNvwBIdv% zk_c*-FltR=VaM>+cO8hgic5+z^U^{7`nKUPi!P``%gXdm4I=C?!GV!c3*s?mW`xHy zK_1fvdF(4I)0YJoCQVV~T)@Z}$;cSBVdDZ;7A=H+hW`xz5#F4D@W({tA{w6m`#|Lu zLQyZsvj}cKDt8jNJO$Uj=otrG7J$ML%s@Ec7g9L?1j!=||Bh<-4-gk6guzBZMHVnJ zdV+n+qVK?>o5Jva5{ocKFtZd*aA0M&+;Etsz+uBimV5`6Tt-G+P;j#3IIv_fGU|ZB zQIFw2E3*+RvmrQv$Sh!F%wc5AWjQs8l`&uvE29`F87D9@=5E-?QtH4`oDPmt&q=I| zLWoq#9HI+}&b-v}oXosbkiEeh4zqNDld&xCm@9ML9Q?axgv;_IdH*+ zNgx{=85x^dnan4#GM<5i9#TeU1P>j7dnJsF&;Ui2W&s5fM4ACq4nV>fSp_p(1;am3 ze;-o?6GFv*kO~kNR=uIH7cepwGcq=B*vO*fz@m@^4xRUtShyB2GA1xGHnXs@Ftad0 zqN|yOjfI7Uc>yD%0V88O3;QHi#&0aFEKHyxhIx_$3zGvY^CSmH2M3m$Frk|l99V8I zm;f5>0gaHevaC@n&P>ls%mEFGizYZafQCpI{>M8yII=P?xG;g`70XK&o=L2XJ(!^e z9$Nzk3&TrjsDbkLXK<2dWmyF62!oRu%V!ob7A_V}NMv*}GIp~vl|qaJ=~>_a>VF~S zEJjcR+R@PwQsi~BGF5>*&dSIFXKBM(olq92b~(Yw*bVA!vD{*1YH@G`RbQa~7Rx#0 w{#H4ts{+np@ctIdW>9~N -0.8 *. (last y') + and z = up(-. y) + +let hybrid main () = + let der t = 1.0 init 0.0 in + let y = ball (y0, y'0) in + let z = period(0.01) in + present z -> ( + print_float t; + print_string "\t"; + print_float y; + print_newline () + ); () diff --git a/exm/zelus/sin_1_x.ml b/exm/zelus/sin_1_x.ml new file mode 100644 index 0000000..d185ae5 --- /dev/null +++ b/exm/zelus/sin_1_x.ml @@ -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 } diff --git a/exm/zelus/sin_1_x.zci b/exm/zelus/sin_1_x.zci new file mode 100644 index 0000000000000000000000000000000000000000..42b4e831f36d4800c963051a35613c59ba1a365c GIT binary patch literal 9262 zcmZpfx@;ccN?L@rLmg7aY`!p{xZC8xA`0 z#tROP6E-a1U|>*ZaA5e)$ne3z!NJji;Xjhldj|&xCy>bw3nsBJPh2pGg=4`a7A}Vg z8yp-KOag12;K=ab$pHiy{!fI7fy`s91_^@9gDP}j`0oN0a#*nJfs@063lkg`1c$O% zKomGQGK4ueID+hh2{1Bz1Ubsd(c!}S0}cx=99-b8RAO*p(nKBxF!*4;{Xi=M}}B9D+Zg{%3i=?31KZ@v2_5s32YWj zND;z-k_%WO;4CjF3#@zri$5wi8JP=VFJMW5FrefDmJSGO0ZWg=0+t#G2SP4jX@oE! z8ea_MaPbN_s|L=h z#h#lO%HS&M;jAhGg@rUI-NIFupg6f0#Z9HiJVu5FWG+Jys1N{^gNTBHArHBnV8{ij zfcY;UZet;g#mKN3!E8Y=8xhPV1hWmnoPuC>BAD$6W(R`Vg#%^uY5|Kg zM95(Qiyef8aB(w&*@|FpLNI3_m{SqVX$a;l1QXVPb7W+ghY+2OV9r4>7a*7mNeK73 z2yOEb%=HN75(IM*g1H#MT#8_>Loin&n9C5% zT=|WEg*6s+8q|KWWZV3a8}*|mO?nE2*O&xk`K->*$xX>Di*NpgNQmT zVA%<0?SZrQFJL(Y=NyKx7O)&xz*4t>Wsk!GmXinwZ#+ySvc z?E+X{zX}onS3R&6E?8y(3l~%Z%6)*$g|I>O3s@Gy1%(QzoepX~urhKwEMWNwl>jSR z!1523yM6&nC@8~&TlOxBC1}L_|y8FkS3Ul*`4$}M2oy?`Zg0gL4V zR>n85P+(+uh+sZKFdrkBPaGDoGRcCf5U@)fklV5t<>+&Sz847QOVk_&bsea&4>l5^ z@D+mj8o_*nL%S@@AXX-0gz5(iSQOzkJ0rtekjodaGFc(Fr`*(vQj1G-;Gw#JMP&gi zWBURYt_7@&O(0{z@x%meCoN!Q@&$zgR02dX{D*bu7#Y4ot()M)@E=?@BPF4@Vut?? z7eG2394?$6G^%gBH!tA&uok_f;}7l`Yiy?*3C1(!h(S?GWSvaBw` zWN}P;!A&8EUNKBreS}^KP`?D^Z*V+eDm6zamBf^_K*&mC%7S}V5JyO1%32}x%3#V` zBV-jYWx?GWh{^JpvbG4likPzC79T{fGN!COLa!30ELakvS7iaq3`EAoRO$#;3o0WR zRWW5@V{(p+jA{;`LIxcCD9vGb+03F0G74OWf)g8<1&%^63tXCkS>P~(jSwQI6H zVl;w=b&bK`vHb!DM+O&G<~yv6+N_LP3u+l$82&r3GTMXcc}D#O6E8S`Dtb^o&-}z4 zRL@tU)bpOKj1w0y_%Zk+hTGwz>Yfb$C$TaeoW#ngeqq}J2Zsv_CoW*&QPS0$;IP2K zVZ&huhYcGSu!z9bGcvX_{D%okBMG-M{NJ#VMZtkZmXWcA;Xfl|8;Iy&_&?EMB4~`9 zmHFWUh9HIzR>tBq)H4n0j9GVR?r4C8URD?rYk%T8h4DUh`nglX@3W%7-@Si1@C1(LB z1%c-m7(Uo9V8~}EVr6NW#L8?miIs8Yf(sK^<}6?+V<=}yoy5xQG>MgQo5KQzJcf!5 z8(BmgScJks>4MpO5-a0ehXoATU~w}C7LyRLxbh@c#tCpQuq>G1z%mcwe3rQmEVCII zXTjXfGRuKw1|#DPkb_yaJFskFWSkCi<65|B%mfru3h|5+Ze+(#I6c!C4VPpCl$9a#1=GOmFe zw9kQM4WSavcpt^}E}8bqvR_|Ni*nb#D?YlM^U+L8T-eSQ*9r z!BWzbSQ)2GfE6{YjHaM)gO(7ijK&VEjE0PiJ7FQm%4p!g%BaW4xE&OBtc zjN3o~w*y4%0)<`*D`PS|CKj+Vb}QF`ig{MXdPsaMU}c68)VpC zkYQDZLong9nv9I_X-!5(_#6QvBYcXJkue9G3s4$4ypZ+*!uSjXGZRI( zFhr%p0v0hiOA^kKfwSb{EG0M#+^GULhrpd6Fbmv?2D8AOdr)J?0@~LGMM1;@77K?3 ztjyr90ophbxM7Z*yxuEmX@OGhV|eQNhLn;H!qV;nR>r?b;`{El-$FkgkWho=$8Iaj5OB`4hF*2S88G06E zs5C2+6exQ!H!NUeILF9vo|VaO5-aN(RwiLk_BzGLaACtnmNgD6s~7@61qSQFNuVi7 zMutO-3>QJg%svNJ#yt$)V7V!iSXx24_AxSC0?9paU}e0=-~pCvn8Y#x$@So>;2I;t zO;$FyNvy2ZlUU|Kb>7^tkp)x@GdP3wRZU`936j~$$Z%`JMph;c2UaGwYOqwuB$gV7 z1q{iI47WgSH|DLNW-61~f{7PcnYVxzXKa|l$`1a7M;My=a3s?ETXJz z(^*7Vgjqzv0rb;>Mbv?n`KN=U!vwf9Z5A*x{AOhM$I3Q;5-W@1BozN9OeQGIS2F}VIykbj@Lrg};?3g4 z!ZL}KMP(9;8jBAL8w(Q)BRGg#7BGSU6HCk_R@Q|qPM`pOz{tq7VFGxRJjsC}#9`tB z7QO{6oFH0Y0gJB#OBxFo3mXeFE7J@Z<2IDRA_-wQEC^;~WOA6WfR*hUE87(ZP`4id Dwh!nG literal 0 HcmV?d00001 diff --git a/exm/zelus/sin_1_x.zls b/exm/zelus/sin_1_x.zls new file mode 100644 index 0000000..48cc71e --- /dev/null +++ b/exm/zelus/sin_1_x.zls @@ -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 diff --git a/exm/zelus/sincos/sincosz.ml b/exm/zelus/sincos/sincosz.ml new file mode 100644 index 0000000..ebbd09d --- /dev/null +++ b/exm/zelus/sincos/sincosz.ml @@ -0,0 +1,45 @@ +(* The Zelus compiler, version 2.2-dev + (2025-06-16-15:24) *) +open Common +open Ztypes +open Solvers + +type ('c , 'b , 'a) _f = + { mutable major_11 : 'c ; mutable sin_13 : 'b ; mutable cos_12 : 'a } + +let f (cstate_14:Ztypes.cstate) = + + let f_alloc _ = + cstate_14.cmax <- (+) cstate_14.cmax 2; + { major_11 = false ; + sin_13 = { pos = 42.; der = 0. } ; cos_12 = { pos = 42.; der = 0. } } in + let f_step self ((_time_10:float) , ()) = + ((let (cindex_15:int) = cstate_14.cindex in + let cpos_17 = ref (cindex_15:int) in + cstate_14.cindex <- (+) cstate_14.cindex 2 ; + self.major_11 <- cstate_14.major ; + (if cstate_14.major then + for i_1 = cindex_15 to 1 do Zls.set cstate_14.dvec i_1 0. done + else ((self.sin_13.pos <- Zls.get cstate_14.cvec !cpos_17 ; + cpos_17 := (+) !cpos_17 1) ; + (self.cos_12.pos <- Zls.get cstate_14.cvec !cpos_17 ; + cpos_17 := (+) !cpos_17 1))) ; + (let (result_19) = + self.cos_12.der <- (~-.) self.sin_13.pos ; + self.sin_13.der <- self.cos_12.pos ; + Bigarray.(Array1.of_array Float64 c_layout + [| self.sin_13.pos; self.cos_12.pos |]) in + cpos_17 := cindex_15 ; + (if cstate_14.major then + (((Zls.set cstate_14.cvec !cpos_17 self.sin_13.pos ; + cpos_17 := (+) !cpos_17 1) ; + (Zls.set cstate_14.cvec !cpos_17 self.cos_12.pos ; + cpos_17 := (+) !cpos_17 1))) + else (((Zls.set cstate_14.dvec !cpos_17 self.sin_13.der ; + cpos_17 := (+) !cpos_17 1) ; + (Zls.set cstate_14.dvec !cpos_17 self.cos_12.der ; + cpos_17 := (+) !cpos_17 1)))) ; result_19))) in + + let f_reset self = + ((self.sin_13.pos <- 0. ; self.cos_12.pos <- 1.):unit) in + Node { alloc = f_alloc; step = f_step ; reset = f_reset } diff --git a/exm/zelus/sincos/sincosz.zci b/exm/zelus/sincos/sincosz.zci new file mode 100644 index 0000000000000000000000000000000000000000..e531551d2cd8905615a83b010c30eb8d1cc752bc GIT binary patch literal 1170 zcmZpfx@;c<14|tP1EVkl1B*KY1B=}S_2A6BW~zW^39h(-rTUk3+Akm)df3dnjV zM~4gR4>&HkaBzXEa&czfg$WB>l#}y|7nGE5n8IS?<_J?A2z8~SKa3F%W5mE1u?`DZ zSU?V8`0oG{^IO2e3u6bt7{SO2SvX+Ap$ib~5QhaUtY~HhEMVb)8yp3aU%(;?(;T^g zg$*tl4weLmn1f>&GA9C=6TN_ibpZ=M$Qp34JD`yW?>K;@U{)n8V3BfIz@p{g2#P)j z2Zsq07C5swGW-V_1S1!)$ShznS-`>x)0Mb@#l~R)iwRs8*hm%|h$fh-IEMu+_7EY5 z1uWpG28AYz>jD-oh)qr`eqa%(s!1$a6D~}e$in~y3u+l0So$WhSi@p29U3YVoEZK) zgG@$vazPd-9Xeb9i!QiukR|EDwgnCj7Zy&0C2-|BUA+ko3mhCa9CmQnuyFwkCp>MJ zG5iM!@WBO28UB|u{GaGB5tQ~=niepGFobWI!m`8NaS}`PB$hc23m5_!!Z&PW5pZDP z^I-Tti6sQ$#);sF0_CD?uuB;JJ7z+{10jPV zM+%}3_6rzd84_5oPGX6G*ze$&=dggq9jajggF8b4!+!@BSBDD@EK3$lyf6WjQyf^9 zUT|PpyI|r4P`Y$ exit 2 +let args = List.rev !modelargs +let () = ignore lift + let m = try match !model with | 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 "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 "ballz" -> lift Ballz.ball + | Some "sincosz" -> lift Sincosz.f | 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) 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 = if !sundials then let open StatefulSundials in @@ -57,7 +83,7 @@ let sim = 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)) + run_until_n (output !sample (run m s)) else let open StatefulRK45 in let c = if !inplace then InPlace.csolve else Functional.csolve in @@ -66,7 +92,7 @@ let sim = let s = Solver.solver_c c z in let open Sim.Sim(val st) in let n = if !accel then accelerate m s else run m (d_of_dc s) in - run_until_n (Output.print !sample n) + run_until_n (output !sample n) let () = sim !stop !steps ignore diff --git a/src/bin/output.ml b/src/bin/output.ml index c397aa9..e60ea97 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -1,7 +1,6 @@ open Hsim.Types open Hsim.Utils -open Common let print_entry y t = let n = Bigarray.Array1.dim y in @@ -13,6 +12,16 @@ let print_entry y t = Format.printf "\n"; flush stdout +let print_entry_h y t h = + 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 = @@ -20,8 +29,16 @@ let print_sample samples ({ h; u; _ }, now) = else if i = samples then print_entry (u h) (now +. h) else let t = float_of_int i *. step in (print_entry (u t) (now +. t); loop (i+1)) in - if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end - else begin Debug.print "C: "; loop 0 end + if h <= 0.0 then print_entry (u 0.0) now else loop 0 + +let print_sample_h samples ({ h; u; _ }, now) = + let step = h /. (float_of_int samples) in + 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 @@ -30,3 +47,18 @@ let print_limits { 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, ()) } + diff --git a/src/lib/common/debug.ml b/src/lib/common/debug.ml index e34b3cf..8db8349 100644 --- a/src/lib/common/debug.ml +++ b/src/lib/common/debug.ml @@ -2,3 +2,10 @@ let debug = ref false 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) diff --git a/src/lib/common/ztypes.ml b/src/lib/common/ztypes.ml new file mode 100644 index 0000000..82ee092 --- /dev/null +++ b/src/lib/common/ztypes.ml @@ -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? +*) + diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 87406bc..ac8ed4f 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -8,6 +8,7 @@ module Sim (S : SimState) = include S let step_discrete s step hor fder fzer cget csize zsize jump reset = + Common.Debug.print "SIMU :: DISCRETE :: start"; 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 o, ms = step ms (i.u now) in @@ -26,9 +27,11 @@ 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 + Common.Debug.print "SIMU :: DISCRETE :: end"; Utils.dot o, s let step_continuous s step cset fout zset = + Common.Debug.print "SIMU :: CONTINUOUS :: start"; 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 @@ -46,6 +49,7 @@ module Sim (S : SimState) = let s = set_running ~mode:Discrete ~now:h s in update (zset ms z) ss s, Discontinuous in let h = h -. now in + Common.Debug.print "SIMU :: CONTINUOUS :: end"; { h; u=fout; c }, s, { h; c; u=fms } (** Simulation of a model with any solver. *) @@ -55,7 +59,7 @@ module Sim (S : SimState) = : ('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 + 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 = @@ -97,7 +101,7 @@ module Sim (S : SimState) = 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 diff --git a/src/lib/hsim/types.ml b/src/lib/hsim/types.ml index 705f20d..23e7b52 100644 --- a/src/lib/hsim/types.ml +++ b/src/lib/hsim/types.ml @@ -69,7 +69,7 @@ type ('p, 'a, 'b) sim = ('p, 'a signal, 'b signal) dnode (** 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 +(** Consider a model without assertions as a model with an empty list of assertions. *) let a_of_h (HNode body) = HNodeA { body; assertions=[] } diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml index 434be45..3f5cdb6 100644 --- a/src/lib/solvers/illinois.ml +++ b/src/lib/solvers/illinois.ml @@ -3,7 +3,9 @@ (* part of the Zelus standard library. *) (* 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 @@ -121,7 +123,7 @@ type t = { let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c = s.t1 <- t; 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); s.bothf_valid <- false @@ -152,6 +154,7 @@ let num_roots { f0; _ } = Zls.length f0 (* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *) let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c = + Common.Debug.print "ZSOL :: Calling [step]"; (* swap f0 and f1; f0 takes the previous value of f1 *) s.f0 <- f1; s.t0 <- t1; @@ -162,7 +165,7 @@ let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c = g t c s.f1; s.bothf_valid <- true; - if !debug then + if debug () then (printf "z|---------- step(%.24e, %.24e)----------@." s.t0 s.t1; log_limits s.f0 s.f1) @@ -212,7 +215,7 @@ let find ({ g = g; bothf_valid = bothf_valid; dky t_right 0; (* c = dky_0(t_right); update state *) ignore (update_roots calc_zc f_left (get_f_right f_right') roots); - if !debug then + if debug () then (printf "z|---------- stall(%.24e, %.24e) {interval < %.24e !}--@." t_left t_right ttol; @@ -280,20 +283,20 @@ let find ({ g = g; bothf_valid = bothf_valid; match check_interval calc_zc f_left f_mid with | SearchLeft -> - if !debug then printf "z| (%.24e -- %.24e] %.24e@." + if debug () then printf "z| (%.24e -- %.24e] %.24e@." t_left t_mid t_right; let alpha = if i >= 1 then alpha *. 0.5 else alpha 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) | SearchRight -> - if !debug then printf "z| %.24e (%.24e -- %.24e]@." + if debug () then printf "z| %.24e (%.24e -- %.24e]@." t_left t_mid t_right; 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) | FoundMid -> - if !debug then printf "z| %.24e [%.24e] %.24e@." + if debug () then printf "z| %.24e [%.24e] %.24e@." t_left t_mid t_right; ignore (update_roots calc_zc f_left f_mid roots); let f_tmp = f_mid_from_f_right f_right' in @@ -303,7 +306,7 @@ let find ({ g = g; bothf_valid = bothf_valid; if not bothf_valid then (clear_roots roots; assert false) else begin - if !debug then + if debug () then printf "z|\nz|---------- find(%.24e, %.24e)----------@." t0 t1; match check_interval calc_zc f0 f1 with @@ -314,7 +317,7 @@ let find ({ g = g; bothf_valid = bothf_valid; end | 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); s.bothf_valid <- false; t1 diff --git a/src/lib/solvers/odexx.ml b/src/lib/solvers/odexx.ml index 31637aa..30252be 100644 --- a/src/lib/solvers/odexx.ml +++ b/src/lib/solvers/odexx.ml @@ -51,7 +51,9 @@ module GenericODE (Butcher : BUTCHER_TABLEAU) : STATE_ODE_SOLVER = struct (* {{{1 *) open Bigarray - let debug = ref false (* !Debug.debug *) + let debug () = + false + (* !Common.Debug.debug *) 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)" 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 = @@ -288,11 +290,11 @@ struct (* {{{1 *) let tnew = if finished then max_t else t +. h *. (mA maxK) in mapinto ynew (make_newval y 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 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 (Printf.sprintf "Error (%e) > relative tolerance (%e) at t=%e" diff --git a/src/lib/solvers/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml index d85658e..997d0ee 100644 --- a/src/lib/solvers/statefulRK45.ml +++ b/src/lib/solvers/statefulRK45.ml @@ -22,6 +22,8 @@ module Functional = { state; vec = init } in let step ({ state ; vec=v } as s) h = + Common.Debug.print "SOLVER STEP"; + Common.Debug.print_entry v; let y_nv = vec v in let h = step state h y_nv in let state = copy state in diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml index cc8c255..f6dc60b 100644 --- a/src/lib/solvers/statefulZ.ml +++ b/src/lib/solvers/statefulZ.ml @@ -15,7 +15,10 @@ module Functional = vec = zmake 0 } in 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 + Common.Debug.print "ZSolver Reset"; + Common.Debug.print_entry init; { state = initialize size fzer init; vec = if length vec = size then vec else zmake size } in From 6d92261afd1ea1f4b84ff865d038b0f325549077 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Mon, 23 Jun 2025 15:48:58 +0200 Subject: [PATCH 5/5] feat: somewhat compatible with zelus output --- exm/ball.ml | 4 +- exm/dune | 5 ++ exm/zelus/ballz/ballz.ml | 131 ------------------------------ exm/zelus/ballz/ballz.ml.bak | 137 -------------------------------- exm/zelus/ballz/ballz.ml.bak2 | 135 ------------------------------- exm/zelus/ballz/ballz.zci | Bin 6309 -> 0 bytes exm/zelus/ballz/ballz.zls | 11 ++- exm/zelus/ballz/dune | 6 ++ exm/zelus/sincos/dune | 6 ++ exm/zelus/sincos/sincosz.ml | 45 ----------- exm/zelus/sincos/sincosz.zci | Bin 1170 -> 0 bytes exm/ztypes.ml | 4 + src/bin/lift.ml | 1 - src/bin/main.ml | 36 ++++++--- src/lib/hsim/sim.ml | 42 +++++----- src/lib/hsim/state.ml | 54 ++++++++----- src/lib/solvers/illinois.ml | 1 - src/lib/solvers/statefulRK45.ml | 2 - src/lib/solvers/statefulZ.ml | 2 - 19 files changed, 107 insertions(+), 515 deletions(-) delete mode 100644 exm/zelus/ballz/ballz.ml delete mode 100644 exm/zelus/ballz/ballz.ml.bak delete mode 100644 exm/zelus/ballz/ballz.ml.bak2 delete mode 100644 exm/zelus/ballz/ballz.zci create mode 100644 exm/zelus/ballz/dune create mode 100644 exm/zelus/sincos/dune delete mode 100644 exm/zelus/sincos/sincosz.ml delete mode 100644 exm/zelus/sincos/sincosz.zci create mode 100644 exm/ztypes.ml diff --git a/exm/ball.ml b/exm/ball.ml index 72e804e..d78109d 100644 --- a/exm/ball.ml +++ b/exm/ball.ml @@ -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; yd 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 horizon _ = max_float let cget s = s.lx @@ -38,7 +38,7 @@ let zset s zin = { s with zin } let step ({ zin; lx; _ } as s) zfalse = 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 } + of_array [| s.lx.{1}; s.lx.{0} |], { zin=zfalse; lx; i=false } let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec = let yd = cmake csize in diff --git a/exm/dune b/exm/dune index 4eecb1e..aa40c69 100644 --- a/exm/dune +++ b/exm/dune @@ -1,3 +1,8 @@ +(env + (dev + (flags + (:standard -w -a)))) + (library (name examples) (libraries hsim solvers)) diff --git a/exm/zelus/ballz/ballz.ml b/exm/zelus/ballz/ballz.ml deleted file mode 100644 index 08df6f6..0000000 --- a/exm/zelus/ballz/ballz.ml +++ /dev/null @@ -1,131 +0,0 @@ -(* The Zelus compiler, version 2.2-dev - (2025-06-16-15:24) *) -open Common -open Ztypes -open Solvers - -let (+=) r v = r := !r + v - -let g = 9.81 - -let y0 = 50. - -let y'0 = 0. - -type ball = - { mutable major : bool; - mutable h : float; - mutable init : bool; - mutable z : (bool, float) zerocrossing; - mutable y' : float continuous; - mutable y : float continuous } - -let ball cstate = - - let ball_alloc _ = - cstate.cmax <- cstate.cmax + 2; - cstate.zmax <- cstate.zmax + 1; - { major = false; - h = 42.; - init = false; - z = { zin = false; zout = 1. }; - y' = { pos = 42.; der = 0. }; - y = { pos = 42.; der = 0. } - } in - - let ball_step self (_, ()) = - let cidx = cstate.cindex in let cpos = ref cidx in - let zidx = cstate.zindex in let zpos = ref zidx in - cstate.cindex <- cstate.cindex + 2; - cstate.zindex <- cstate.zindex + 1; - self.major <- cstate.major; - if cstate.major then - for i = cidx to 1 do Zls.set cstate.dvec i 0. done - else begin - self.y'.pos <- Zls.get cstate.cvec !cpos; cpos += 1; - self.y.pos <- Zls.get cstate.cvec !cpos; cpos += 1 - end; - let res0, res1 = - let encore = ref false in - if self.init then self.y'.pos <- y'0; - let last_y' = self.y'.pos in - if self.z.zin then begin - encore := true; self.y'.pos <- -0.8 *. last_y' - end; - self.h <- if !encore then 0. else infinity; - if self.init then self.y.pos <- y0; - self.init <- false; - self.z.zout <- -. self.y.pos; - self.y'.der <- -. g; - self.y.der <- self.y'.pos; - self.y.pos, self.y.der - in - cstate.horizon <- min cstate.horizon self.h; - cpos := cidx; - if cstate.major then begin - Zls.set cstate.cvec !cpos self.y'.pos; cpos += 1; - Zls.set cstate.cvec !cpos self.y.pos; cpos += 1; - self.z.zin <- false - end else begin - self.z.zin <- Zls.get_zin cstate.zinvec !zpos; zpos += 1; - zpos := zidx; - Zls.set cstate.zoutvec !zpos self.z.zout; zpos += 1; - Zls.set cstate.dvec !cpos self.y'.der; cpos += 1; - Zls.set cstate.dvec !cpos self.y.der; cpos += 1 - end; - Bigarray.(Array1.of_array Float64 c_layout [| res0; res1 |]) in - - let ball_reset self = self.init <- true in - - Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } - -type ('f , 'e , 'd , 'c , 'b , 'a) _main = - { mutable i_73 : 'f ; - mutable major_62 : 'e ; - mutable h_72 : 'd ; - mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } - -let main (cstate_80:Ztypes.cstate) = - let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball - cstate_80 in - let main_alloc _ = - cstate_80.cmax <- (+) cstate_80.cmax 1; - { major_62 = false ; - h_72 = 42. ; - i_70 = (false:bool) ; - h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; - i_73 = i_73_alloc () (* continuous *) } in - let main_step self ((time_61:float) , ()) = - ((let (cindex_81:int) = cstate_80.cindex in - let cpos_83 = ref (cindex_81:int) in - cstate_80.cindex <- (+) cstate_80.cindex 1 ; - self.major_62 <- cstate_80.major ; - (if cstate_80.major then - for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done - else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; - cpos_83 := (+) !cpos_83 1))) ; - (let (result_85) = - let h_71 = ref (infinity:float) in - (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; - (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in - self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; - h_71 := min !h_71 self.h_68 ; - self.h_72 <- !h_71 ; - self.i_70 <- false ; - self.t_63.der <- 1. ; - (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in - (begin match z_69 with - | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 - | _ -> () end) ; - Bigarray.(Array1.create Float64 c_layout 0))) in - cstate_80.horizon <- min cstate_80.horizon self.h_72 ; - cpos_83 := cindex_81 ; - (if cstate_80.major then - (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; - cpos_83 := (+) !cpos_83 1))) - else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; - cpos_83 := (+) !cpos_83 1)))) ; result_85))) in - let main_reset self = - ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): - unit) in - Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak b/exm/zelus/ballz/ballz.ml.bak deleted file mode 100644 index d248fc0..0000000 --- a/exm/zelus/ballz/ballz.ml.bak +++ /dev/null @@ -1,137 +0,0 @@ -(* The Zelus compiler, version 2.2-dev - (2025-06-16-15:24) *) -open Common -open Ztypes -open Solvers - -let g = 9.81 - -let y0 = 50. - -let y'0 = 0. - -type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = - { mutable major_50 : 'g ; - mutable h_60 : 'f ; - mutable h_58 : 'e ; - mutable i_56 : 'd ; - mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } - -let ball (cstate_74:Ztypes.cstate) = - - let ball_alloc _ = - cstate_74.cmax <- (+) cstate_74.cmax 2 ; - cstate_74.zmax <- (+) cstate_74.zmax 1; - { major_50 = false ; - h_60 = 42. ; - h_58 = (42.:float) ; - i_56 = (false:bool) ; - x_55 = { zin = false; zout = 1. } ; - y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in - let ball_step self ((_time_49:float) , ((y0_48:float) , (y'0_47:float))) = - Printf.printf "STEP (%d)\n" cstate_74.cindex; - ((let (cindex_75:int) = cstate_74.cindex in - let cpos_77 = ref (cindex_75:int) in - let (zindex_76:int) = cstate_74.zindex in - let zpos_78 = ref (zindex_76:int) in - cstate_74.cindex <- (+) cstate_74.cindex 2 ; - cstate_74.zindex <- (+) cstate_74.zindex 1 ; - self.major_50 <- cstate_74.major ; - (if cstate_74.major then - for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done - else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; - cpos_77 := (+) !cpos_77 1) ; - (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; - cpos_77 := (+) !cpos_77 1))) ; - (let (result_79:float) = - let h_59 = ref (infinity:float) in - let encore_57 = ref (false:bool) in - (if self.i_56 then self.y'_52.pos <- y'0_47) ; - (let (l_54:float) = self.y'_52.pos in - (begin match self.x_55.zin with - | true -> - encore_57 := true ; - self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end) - ; - self.h_58 <- (if !encore_57 then 0. else infinity) ; - h_59 := min !h_59 self.h_58 ; - self.h_60 <- !h_59 ; - (if self.i_56 then self.y_51.pos <- y0_48) ; - self.i_56 <- false ; - self.x_55.zout <- (~-.) self.y_51.pos ; - self.y'_52.der <- (~-.) g ; - self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in - cstate_74.horizon <- min cstate_74.horizon self.h_60 ; - cpos_77 := cindex_75 ; - (if cstate_74.major then - (((Printf.printf "idx: %d\n" !cpos_77; - Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; - cpos_77 := (+) !cpos_77 1) ; - (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; - cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) - else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; - zpos_78 := (+) !zpos_78 1)) ; - zpos_78 := zindex_76 ; - ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; - zpos_78 := (+) !zpos_78 1)) ; - ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; - cpos_77 := (+) !cpos_77 1) ; - (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; - cpos_77 := (+) !cpos_77 1)))) ; result_79)):float) in - let ball_reset self = - (self.i_56 <- true:unit) in - Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } -type ('f , 'e , 'd , 'c , 'b , 'a) _main = - { mutable main_ball : 'f ; - mutable main_major : 'e ; - mutable h_72 : 'd; mutable i_70 : 'c; mutable h : 'b; - mutable t_63 : 'a } - -let main cs = - let Node - { alloc = ball_alloc; - step = ball_step; - reset = ball_reset } = ball cs in - - let main_alloc _ = - cs.cmax <- cs.cmax + 1; - { main_major = false; - h_72 = 42.0; i_70 = false; h = 42.0; - t_63 = { pos = 42.; der = 0. }; - main_ball = ball_alloc () } in - - let main_step self (time, ()) = - let cindex = cs.cindex in - let cpos = ref cindex in - Printf.printf "main:cindex: %d\n" cs.cindex; - cs.cindex <- cs.cindex + 1; - self.main_major <- cs.major; - if cs.major then for i = cindex to 0 do Zls.set cs.dvec i 0. done - else begin self.t_63.pos <- Zls.get cs.cvec !cpos; cpos := !cpos + 1 end; - let result = - if self.i_70 then self.h <- time; - let z = self.main_major && (time >= self.h) in - if z then self.h <- self.h +. 0.01; - self.h_72 <- min infinity self.h; - self.i_70 <- false; - self.t_63.der <- 1.; - let y_64 = ball_step self.main_ball (time, (y0, y'0)) in - if z then begin - print_float self.t_63.pos; - print_string "\t"; - print_float y_64; - print_newline () - end; - Bigarray.(Array1.create Float64 c_layout 0) in - cs.horizon <- min cs.horizon self.h_72; - cpos := cindex; - if cs.major then Zls.set cs.cvec !cpos self.t_63.pos - else Zls.set cs.dvec !cpos self.t_63.der; - result in - - let main_reset self = - self.i_70 <- true; - self.t_63.pos <- 0.; - ball_reset self.main_ball in - - Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.ml.bak2 b/exm/zelus/ballz/ballz.ml.bak2 deleted file mode 100644 index 05842bb..0000000 --- a/exm/zelus/ballz/ballz.ml.bak2 +++ /dev/null @@ -1,135 +0,0 @@ -(* The Zelus compiler, version 2.2-dev - (2025-06-16-15:24) *) -open Common -open Ztypes -open Solvers - -let g = 9.81 - -let y0 = 50. - -let y'0 = 0. - -type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball = - { mutable major_50 : 'g ; - mutable h_60 : 'f ; - mutable h_58 : 'e ; - mutable i_56 : 'd ; - mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a } - -let ball (cstate_74:Ztypes.cstate) = - - let ball_alloc _ = - cstate_74.cmax <- (+) cstate_74.cmax 2 ; - cstate_74.zmax <- (+) cstate_74.zmax 1; - { major_50 = false ; - h_60 = 42. ; - h_58 = (42.:float) ; - i_56 = (false:bool) ; - x_55 = { zin = false; zout = 1. } ; - y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in - - let ball_step self ((_time_49:float) , ()) = - ((let (cindex_75:int) = cstate_74.cindex in - let cpos_77 = ref (cindex_75:int) in - let (zindex_76:int) = cstate_74.zindex in - let zpos_78 = ref (zindex_76:int) in - cstate_74.cindex <- (+) cstate_74.cindex 2 ; - cstate_74.zindex <- (+) cstate_74.zindex 1 ; - self.major_50 <- cstate_74.major ; - (if cstate_74.major then - for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done - else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ; - cpos_77 := (+) !cpos_77 1) ; - (self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ; - cpos_77 := (+) !cpos_77 1))) ; - (let (result_79:float) = - let h_59 = ref (infinity:float) in - let encore_57 = ref (false:bool) in - (if self.i_56 then self.y'_52.pos <- y'0) ; - (let (l_54:float) = self.y'_52.pos in - begin match self.x_55.zin with - | true -> - encore_57 := true ; - self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end; - self.h_58 <- (if !encore_57 then 0. else infinity) ; - h_59 := min !h_59 self.h_58 ; - self.h_60 <- !h_59 ; - (if self.i_56 then self.y_51.pos <- y0) ; - self.i_56 <- false ; - self.x_55.zout <- (~-.) self.y_51.pos ; - self.y'_52.der <- (~-.) g ; - self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in - cstate_74.horizon <- min cstate_74.horizon self.h_60 ; - cpos_77 := cindex_75 ; - (if cstate_74.major then - (((Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ; - cpos_77 := (+) !cpos_77 1) ; - (Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ; - cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false))) - else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ; - zpos_78 := (+) !zpos_78 1)) ; - zpos_78 := zindex_76 ; - ((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ; - zpos_78 := (+) !zpos_78 1)) ; - ((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ; - cpos_77 := (+) !cpos_77 1) ; - (Zls.set cstate_74.dvec !cpos_77 self.y_51.der ; - cpos_77 := (+) !cpos_77 1)))) ; - Bigarray.(Array1.of_array Float64 c_layout [| result_79 |])))) in - - let ball_reset self = - (self.i_56 <- true:unit) in - - Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset } - -type ('f , 'e , 'd , 'c , 'b , 'a) _main = - { mutable i_73 : 'f ; - mutable major_62 : 'e ; - mutable h_72 : 'd ; - mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a } - -let main (cstate_80:Ztypes.cstate) = - let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball - cstate_80 in - let main_alloc _ = - cstate_80.cmax <- (+) cstate_80.cmax 1; - { major_62 = false ; - h_72 = 42. ; - i_70 = (false:bool) ; - h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. }; - i_73 = i_73_alloc () (* continuous *) } in - let main_step self ((time_61:float) , ()) = - ((let (cindex_81:int) = cstate_80.cindex in - let cpos_83 = ref (cindex_81:int) in - cstate_80.cindex <- (+) cstate_80.cindex 1 ; - self.major_62 <- cstate_80.major ; - (if cstate_80.major then - for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done - else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ; - cpos_83 := (+) !cpos_83 1))) ; - (let (result_85) = - let h_71 = ref (infinity:float) in - (if self.i_70 then self.h_68 <- (+.) time_61 0.) ; - (let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in - self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ; - h_71 := min !h_71 self.h_68 ; - self.h_72 <- !h_71 ; - self.i_70 <- false ; - self.t_63.der <- 1. ; - (let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in - (begin match z_69 with - | true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64 - | _ -> () end) ; - Bigarray.(Array1.create Float64 c_layout 0))) in - cstate_80.horizon <- min cstate_80.horizon self.h_72 ; - cpos_83 := cindex_81 ; - (if cstate_80.major then - (((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ; - cpos_83 := (+) !cpos_83 1))) - else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ; - cpos_83 := (+) !cpos_83 1)))) ; result_85))) in - let main_reset self = - ((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ): - unit) in - Node { alloc = main_alloc; step = main_step ; reset = main_reset } diff --git a/exm/zelus/ballz/ballz.zci b/exm/zelus/ballz/ballz.zci deleted file mode 100644 index 08fa125c8a8c6a515d7c97cca6557fe41356e53b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6309 zcmZpfx@;cI57NoU}WTVT(DsZ3nwEZpQ8gKBfkS9Bah>TDJ*ZaA5e)km2Cq;N-yYA0d+N;Naj4GRa}V#08UBSQkuU;dGe9 z!ZBfk16TyCb%GPae`g1fN{0UnCb6(MOo9kIF#LA`34=^r_Q1(u!G*~SE=*bAqf}|I zpt^R$6c%F-M~G4f$0`S?rH+*gSXkg32A2gatO%YHOc8^#14x;J1ITO^4!DXMxB^EQ zs}9Csh=#L#;VeHm%OB1PfV1G{GQ`5g;v5#Rh=E)Mw%lO>3kS%L4F4S#uy8v#I82zZ zz)i6V%#~ZfA`F)kgR>+dESQ1yFa|>mj8y?+RD*r#2#OPgQz93zusbYZQG+RBh+4qH z2IoY;Ss`#%AeU2FIs^BSSD;Q|JN~6^8{ZtP5BW z9>eEZ77kEgEMSp`o9zi@fh}IZq5}y#C*?|Y190g1E?_Z%Dp|ndzko#qq+sF#SEWjI zhW`s#j25t1IV@m_U%(OtQMG`@$YB9X9E9VrfF%*6egX?O*swH+@B)_H1uO{*SfUrO zApD;K^HeR2!Qi@prF;PkDD>GN8XXp}aKdHN7O-$REMTc{aC8J`4N#oBC{}_kLYR@f zfJGD$X^Ef|12-W7g_nfFOGV+ugTt5;6nIdVK*Ng_nG0brV1c<7ss9v|$R%Q%}c9EHwyqYaAA^ zoCd3NSWv~_3@(RQ_AOxHb?|aK{|*;Ge1`=W4zfsI*tWpI;ljd+u;NjNrP6h7(xg@(e=y zIQFGaf+D zjVxCuIIvuT#1+d$2bS}Uj544|Vmar)a)yyn8Wb_I4F6e9vz%IRVFD{7;{ry8EsP9X zH%wt=V(tXS*v2xI#ehX0;`wb1{~cHvSsgAourkUmn0R4=gM)*k z1B(KPudra^1r}CN)y2xFc)@{X$%2U&96$;jSe7|BIyo<3aARcH2J*Y!0tPn*cNX1A ztc;~B+AORrdMu1A+8{H0Kzcy-vmCtOz{w*V4R z6C4@-JAr)R0E#{)hXs%}>I6`To?v7+#Uck$=itbY>408KUI*Fb3=0?;o-i^zXJzVu z8L~p4SE>;X@mBCq+a8?Cs&_b*Ol}_;T-vK3PAn?bE%u!|S4 z2rXdYhnrXqXH~;lRoHETCFBxlX$WrFD#P6XZR0vDV9|s`f(Nus-T*V0kzv{b7DJ>+ zn1UM2lTo=-QMofniv))SEK?B{H!NV8ie&M8u+c1YppFN18(A_yMIxwVU%)aA)Ww6y zEnryzV}M%bpbj}GIKe>)D-IYLZaXYsWtM^{1eITx85y2~Dgh={hYJoY#}-Vyz{;oq zssuJnVP&?Dc64xDz{qf(k+J*2`U4K&igQ6Ai`fE3hL4O4pI8}vCb1+S41a*+@H?QS zhVa&Vus@&?!^rR#mHQi&`v=5@x)JPtC>zp(KoWt~)~t-W6C7BX99fxjSeZ;%nT!`) zm;q^nGgm298h}gc-;4}@S()=Du~Z=(`~=Ct4-pOqjZ-i(+ymJP4?+-!;Xl;SZcyhF zZ0dhTMkZF~`bjMC@|%(28Iq}w&`rG$G6>C7P(=J#{wfWIQp7F(H8)!Zdn;i9UL4$(Kj~+ z5`Aw`qA!T${{luvUPeYfRwmU+tc+Zc;6SctUPImJ$jI;l>KA0`R|x6nh#UqIWcYu9 zgMoovaDfA;TB;UV=l-y&#eM-JqXZ+P3@h^c425te69o#|>793a*%*x2(;IM#^QJImk-N9kO zg=LWDJ&WN2Mn-K$MjaNhNvwBIdv% zk_c*-FltR=VaM>+cO8hgic5+z^U^{7`nKUPi!P``%gXdm4I=C?!GV!c3*s?mW`xHy zK_1fvdF(4I)0YJoCQVV~T)@Z}$;cSBVdDZ;7A=H+hW`xz5#F4D@W({tA{w6m`#|Lu zLQyZsvj}cKDt8jNJO$Uj=otrG7J$ML%s@Ec7g9L?1j!=||Bh<-4-gk6guzBZMHVnJ zdV+n+qVK?>o5Jva5{ocKFtZd*aA0M&+;Etsz+uBimV5`6Tt-G+P;j#3IIv_fGU|ZB zQIFw2E3*+RvmrQv$Sh!F%wc5AWjQs8l`&uvE29`F87D9@=5E-?QtH4`oDPmt&q=I| zLWoq#9HI+}&b-v}oXosbkiEeh4zqNDld&xCm@9ML9Q?axgv;_IdH*+ zNgx{=85x^dnan4#GM<5i9#TeU1P>j7dnJsF&;Ui2W&s5fM4ACq4nV>fSp_p(1;am3 ze;-o?6GFv*kO~kNR=uIH7cepwGcq=B*vO*fz@m@^4xRUtShyB2GA1xGHnXs@Ftad0 zqN|yOjfI7Uc>yD%0V88O3;QHi#&0aFEKHyxhIx_$3zGvY^CSmH2M3m$Frk|l99V8I zm;f5>0gaHevaC@n&P>ls%mEFGizYZafQCpI{>M8yII=P?xG;g`70XK&o=L2XJ(!^e z9$Nzk3&TrjsDbkLXK<2dWmyF62!oRu%V!ob7A_V}NMv*}GIp~vl|qaJ=~>_a>VF~S zEJjcR+R@PwQsi~BGF5>*&dSIFXKBM(olq92b~(Yw*bVA!vD{*1YH@G`RbQa~7Rx#0 w{#H4ts{+np@ctIdW>9~N -0.8 *. (last y') and z = up(-. y) let hybrid main () = let der t = 1.0 init 0.0 in - let y = ball (y0, y'0) in - let z = period(0.01) in - present z -> ( + 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 () ); () diff --git a/exm/zelus/ballz/dune b/exm/zelus/ballz/dune new file mode 100644 index 0000000..81087c1 --- /dev/null +++ b/exm/zelus/ballz/dune @@ -0,0 +1,6 @@ +(rule + (targets ballz.ml ballz.zci) + (deps + (:zl ballz.zls)) + (action + (run zeluc %{zl}))) diff --git a/exm/zelus/sincos/dune b/exm/zelus/sincos/dune new file mode 100644 index 0000000..ef360bb --- /dev/null +++ b/exm/zelus/sincos/dune @@ -0,0 +1,6 @@ +(rule + (targets sincosz.ml sincosz.zci) + (deps + (:zl sincosz.zls)) + (action + (run zeluc %{zl}))) diff --git a/exm/zelus/sincos/sincosz.ml b/exm/zelus/sincos/sincosz.ml deleted file mode 100644 index ebbd09d..0000000 --- a/exm/zelus/sincos/sincosz.ml +++ /dev/null @@ -1,45 +0,0 @@ -(* The Zelus compiler, version 2.2-dev - (2025-06-16-15:24) *) -open Common -open Ztypes -open Solvers - -type ('c , 'b , 'a) _f = - { mutable major_11 : 'c ; mutable sin_13 : 'b ; mutable cos_12 : 'a } - -let f (cstate_14:Ztypes.cstate) = - - let f_alloc _ = - cstate_14.cmax <- (+) cstate_14.cmax 2; - { major_11 = false ; - sin_13 = { pos = 42.; der = 0. } ; cos_12 = { pos = 42.; der = 0. } } in - let f_step self ((_time_10:float) , ()) = - ((let (cindex_15:int) = cstate_14.cindex in - let cpos_17 = ref (cindex_15:int) in - cstate_14.cindex <- (+) cstate_14.cindex 2 ; - self.major_11 <- cstate_14.major ; - (if cstate_14.major then - for i_1 = cindex_15 to 1 do Zls.set cstate_14.dvec i_1 0. done - else ((self.sin_13.pos <- Zls.get cstate_14.cvec !cpos_17 ; - cpos_17 := (+) !cpos_17 1) ; - (self.cos_12.pos <- Zls.get cstate_14.cvec !cpos_17 ; - cpos_17 := (+) !cpos_17 1))) ; - (let (result_19) = - self.cos_12.der <- (~-.) self.sin_13.pos ; - self.sin_13.der <- self.cos_12.pos ; - Bigarray.(Array1.of_array Float64 c_layout - [| self.sin_13.pos; self.cos_12.pos |]) in - cpos_17 := cindex_15 ; - (if cstate_14.major then - (((Zls.set cstate_14.cvec !cpos_17 self.sin_13.pos ; - cpos_17 := (+) !cpos_17 1) ; - (Zls.set cstate_14.cvec !cpos_17 self.cos_12.pos ; - cpos_17 := (+) !cpos_17 1))) - else (((Zls.set cstate_14.dvec !cpos_17 self.sin_13.der ; - cpos_17 := (+) !cpos_17 1) ; - (Zls.set cstate_14.dvec !cpos_17 self.cos_12.der ; - cpos_17 := (+) !cpos_17 1)))) ; result_19))) in - - let f_reset self = - ((self.sin_13.pos <- 0. ; self.cos_12.pos <- 1.):unit) in - Node { alloc = f_alloc; step = f_step ; reset = f_reset } diff --git a/exm/zelus/sincos/sincosz.zci b/exm/zelus/sincos/sincosz.zci deleted file mode 100644 index e531551d2cd8905615a83b010c30eb8d1cc752bc..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1170 zcmZpfx@;c<14|tP1EVkl1B*KY1B=}S_2A6BW~zW^39h(-rTUk3+Akm)df3dnjV zM~4gR4>&HkaBzXEa&czfg$WB>l#}y|7nGE5n8IS?<_J?A2z8~SKa3F%W5mE1u?`DZ zSU?V8`0oG{^IO2e3u6bt7{SO2SvX+Ap$ib~5QhaUtY~HhEMVb)8yp3aU%(;?(;T^g zg$*tl4weLmn1f>&GA9C=6TN_ibpZ=M$Qp34JD`yW?>K;@U{)n8V3BfIz@p{g2#P)j z2Zsq07C5swGW-V_1S1!)$ShznS-`>x)0Mb@#l~R)iwRs8*hm%|h$fh-IEMu+_7EY5 z1uWpG28AYz>jD-oh)qr`eqa%(s!1$a6D~}e$in~y3u+l0So$WhSi@p29U3YVoEZK) zgG@$vazPd-9Xeb9i!QiukR|EDwgnCj7Zy&0C2-|BUA+ko3mhCa9CmQnuyFwkCp>MJ zG5iM!@WBO28UB|u{GaGB5tQ~=niepGFobWI!m`8NaS}`PB$hc23m5_!!Z&PW5pZDP z^I-Tti6sQ$#);sF0_CD?uuB;JJ7z+{10jPV zM+%}3_6rzd84_5oPGX6G*ze$&=dggq9jajggF8b4!+!@BSBDD@EK3$lyf6WjQyf^9 zUT|PpyI|r4P`Y$ exit 2 let args = List.rev !modelargs let () = ignore lift +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 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 "ballz" -> lift Ballz.ball - | Some "sincosz" -> lift Sincosz.f - | Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2 + 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) diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index ac8ed4f..d6dd628 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -7,12 +7,13 @@ module Sim (S : SimState) = struct include S - let step_discrete s step hor fder fzer cget csize zsize jump reset = - Common.Debug.print "SIMU :: DISCRETE :: start"; - let ms, ss = get_mstate s, get_sstate s in + 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 @@ -27,30 +28,25 @@ 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 - Common.Debug.print "SIMU :: DISCRETE :: end"; Utils.dot o, s - let step_continuous s step cset fout zset = - Common.Debug.print "SIMU :: CONTINUOUS :: start"; + 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 (h, f, z), ss = step ss stop 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 -> - let s, c = if h >= stop - then set_running ~mode:Discrete ~now:h s, Discontinuous - else set_running ~now: h s, Continuous in - update ms ss s, c - | Some z -> - let s = set_running ~mode:Discrete ~now:h s in - update (zset ms z) ss s, Discontinuous in + 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 - Common.Debug.print "SIMU :: CONTINUOUS :: end"; - { h; u=fout; c }, s, { h; c; u=fms } + { h; u=fout; c }, update ms ss (set_zin z s), { h; c; u=fms } (** Simulation of a model with any solver. *) let run @@ -59,11 +55,11 @@ module Sim (S : SimState) = : ('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 + 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.zset in + let o, s, _ = step_continuous st s.step m.cset m.fout m.horizon in Some o, s in let step st = function @@ -95,8 +91,8 @@ module Sim (S : SimState) = 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 + 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 @@ -104,7 +100,7 @@ module Sim (S : SimState) = 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 + 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 = @@ -140,10 +136,10 @@ module Sim (S : SimState) = = 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.csize m.zsize m.jump s.reset in + 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.zset in + let o, st, _ = step_continuous st s.step m.cset m.fout m.horizon in o, st in let rec step st = function diff --git a/src/lib/hsim/state.ml b/src/lib/hsim/state.ml index 676e12a..29fd8b2 100644 --- a/src/lib/hsim/state.ml +++ b/src/lib/hsim/state.ml @@ -15,59 +15,65 @@ module type SimState = - Idle: waiting for input; - Running: currently integrating; in this case, we have access to the step mode, current input, timestamp and stop time. *) - type ('a, 'ms, 'ss) state + type ('a, 'ms, 'ss, 'zin) 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. *) - 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. ⚠ 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. ⚠ 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. ⚠ 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. ⚠ 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. *) - 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 ? *) - val is_running : ('a, 'ms, 'ss) state -> bool + val is_running : ('a, 'ms, 'ss, 'zin) state -> bool (** 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. *) - 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. *) - 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. *) val set_running : ?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. *) - 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 module type SimStateCopy = sig 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 module FunctionalSimState : SimState = @@ -88,15 +94,17 @@ module FunctionalSimState : SimState = (** Internal state of the simulation node: model state, solver state and current simulation status. *) - type ('a, 'ms, 'ss) state = + type ('a, 'ms, 'ss, 'zin) state = { status : 'a status; (** Current simulation status. *) mstate : 'ms; (** Model state. *) - sstate : 'ss } (** Solver state. *) + sstate : 'ss; (** Solver state. *) + zin : 'zin option; } (** Last zero-crossing vector *) exception Not_running let get_mstate state = state.mstate let get_sstate state = state.sstate + let get_zin state = state.zin let is_running state = 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_sstate sstate state = { state with sstate } + let set_zin zin state = { state with zin } let update mstate sstate state = { state with mstate; sstate } @@ -132,7 +141,7 @@ module FunctionalSimState : SimState = let get_stop s = 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 module InPlaceSimState : SimState = @@ -146,15 +155,17 @@ module InPlaceSimState : SimState = mutable stop : time; } -> 'a status - type ('a, 'ms, 'ss) state = + type ('a, 'ms, 'ss, 'zin) state = { mutable status : 'a status; mutable mstate : 'ms; - mutable sstate : 'ss } + mutable sstate : 'ss; + mutable zin : 'zin option } exception Not_running let get_mstate state = state.mstate let get_sstate state = state.sstate + let get_zin state = state.zin let is_running state = 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_sstate sstate state = state.sstate <- sstate; state + let set_zin zin state = state.zin <- zin; state let update mstate sstate state = state.mstate <- mstate; state.sstate <- sstate; state @@ -192,6 +204,6 @@ module InPlaceSimState : SimState = let get_stop s = 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 diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml index 3f5cdb6..7e280ce 100644 --- a/src/lib/solvers/illinois.ml +++ b/src/lib/solvers/illinois.ml @@ -154,7 +154,6 @@ let num_roots { f0; _ } = Zls.length f0 (* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *) let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c = - Common.Debug.print "ZSOL :: Calling [step]"; (* swap f0 and f1; f0 takes the previous value of f1 *) s.f0 <- f1; s.t0 <- t1; diff --git a/src/lib/solvers/statefulRK45.ml b/src/lib/solvers/statefulRK45.ml index 997d0ee..d85658e 100644 --- a/src/lib/solvers/statefulRK45.ml +++ b/src/lib/solvers/statefulRK45.ml @@ -22,8 +22,6 @@ module Functional = { state; vec = init } in let step ({ state ; vec=v } as s) h = - Common.Debug.print "SOLVER STEP"; - Common.Debug.print_entry v; let y_nv = vec v in let h = step state h y_nv in let state = copy state in diff --git a/src/lib/solvers/statefulZ.ml b/src/lib/solvers/statefulZ.ml index f6dc60b..b642a08 100644 --- a/src/lib/solvers/statefulZ.ml +++ b/src/lib/solvers/statefulZ.ml @@ -17,8 +17,6 @@ module Functional = let reset { fzer; init; size } { vec; _ } = let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in - Common.Debug.print "ZSolver Reset"; - Common.Debug.print_entry init; { state = initialize size fzer init; vec = if length vec = size then vec else zmake size } in