From b037dacccf3c6b579d67f909c4a310d2dca3b9cd Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Fri, 25 Apr 2025 16:41:41 +0200 Subject: [PATCH] feat: greedy simulation --- doc/notes.typ | 107 ++++++++++++++++++++++++++++++++++++++++ src/bin/main.ml | 16 ++++-- src/bin/output.ml | 9 +++- src/lib/common/debug.ml | 6 +++ src/lib/hsim/sim.ml | 24 ++++++--- src/lib/hsim/utils.ml | 19 +++++-- 6 files changed, 164 insertions(+), 17 deletions(-) create mode 100644 src/lib/common/debug.ml diff --git a/doc/notes.typ b/doc/notes.typ index b552d9a..0d40e28 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -300,6 +300,113 @@ Two possible options for the simulation reset: makes this impossible. We thus need reset parameters for both the model and solver. +=== Steps + +The _lazy simulation_'s step function is as follows: + +```ocaml +let step s i = + let ms, ss = S.get_mstate s, S.get_sstate s in + match i, S.is_running s with + | Some i, _ -> + let mode, now, stop = Discrete, 0.0, i.length in + None, S.set_running ~mode ~input:i ~now ~stop s + | None, false -> None, s + | None, true -> + let i, now, stop = S.get_input s, S.get_now s, S.get_stop s in + match S.get_mode s with + | Discrete -> + let o, ms = model.step ms (i.u now) in + let s = + let h = model.horizon ms in + if h <= 0.0 then S.set_mstate ms s + else if now >= stop then S.set_idle s + else if model.jump ms then + let init = model.cget ms in + let fder t = model.fder ms (Utils.offset i now t) in + let fzer t = model.fzer ms (Utils.offset i now t) in + let ivp = { fder; stop = stop -. now; init } in + let zc = { init ; fzer; size = model.zsize } in + let ss = solver.reset (ivp, zc) ss in + let i = { start=i.start +. now; length=i.length -. now; + u=Utils.offset i now } in + let mode, stop, now = Continuous, i.length, 0.0 in + S.update ms ss (S.set_running ~mode ~input:i ~stop ~now s) + else S.set_running ~mode:Continuous s in + Some { start = i.start+. now; length = 0.0; u = fun _ -> o }, s + | Continuous -> + let (h, f, z), ss = solver.step ss stop in + let ms = model.cset ms (f h) in + let start = i.start +. now in + let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in + let out = { start; length=h -. now; u=fout } in + let s = match z with + | None -> + let s = if h >= stop + then S.set_running ~mode:Discrete ~now:h s + else S.set_running ~now:h s in + S.update ms ss s + | Some z -> + let s = S.set_running ~mode:Discrete ~now:h s in + S.update (model.zset ms z) ss s in + Some out, s in +``` + +The _greedy simulation_'s step function is as follows: + +```ocaml +let rec step s i = + let ms, ss = S.get_mstate s, S.get_sstate s in + if not (S.is_running s) then + let mode, now, stop = Discrete, 0.0, i.length in + step (S.set_running ~mode ~input:i ~now ~stop s) i + else let now, stop = S.get_now s, S.get_stop s in + match S.get_mode s with + | Discrete -> + let o, ms = model.step ms (i.u now) in + let h = model.horizon ms in + let rest, s = + if h <= 0.0 then step (S.set_mstate ms s) i + else if now >= stop then [], s + else if model.jump ms then + let init = model.cget ms in + let fder t = model.fder ms (Utils.offset i now t) in + let fzer t = model.fzer ms (Utils.offset i now t) in + let ivp = { fder; stop = stop -. now; init } in + let zc = { init; fzer; size = model.zsize } in + let ss = solver.reset (ivp, zc) ss in + let i = { start=i.start +. now; length=i.length -. now; + u=Utils.offset i now } in + let mode, stop, now = Continuous, i.length, 0.0 in + let s = S.set_running ~mode ~input:i ~stop ~now s in + step (S.update ms ss s) i + else step (S.set_running ~mode:Continuous s) i in + { start = i.start +. now; length = 0.0; u = fun _ -> o }::rest, s + | Continuous -> + let (h, f, z), ss = solver.step ss stop in + let ss = solver.copy ss in + let ms = model.cset ms (f h) in + let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in + let out = { start = i.start +. now; length = h -. now; u = fout } in + match z with + | None -> + if h >= stop then + let s = S.set_running ~mode:Discrete ~now:h s in + let rest, s = step (S.update ms ss s) i in + out::rest, s + else + let s = S.set_running ~now:h s in + let rest, s = step (S.update ms ss s) i in + (match rest with + | [] -> [out], s + | f::rest -> Utils.compose [out;f] :: rest, s) + | Some z -> + let s = S.set_running ~mode:Discrete ~now:h s in + let ms = model.zset ms z in + let rest, s = step (S.update ms ss s) i in + out::rest, s in +``` + == Mathematical model #link("https://zelus.di.ens.fr/cc2015/fullpaper.pdf")[[CC'15]] defines the diff --git a/src/bin/main.ml b/src/bin/main.ml index 13d5094..8487520 100644 --- a/src/bin/main.ml +++ b/src/bin/main.ml @@ -1,19 +1,22 @@ open Hsim open Examples +open Common let sample = ref 10 let stop = ref 30.0 -let debug = ref false +let greedy = ref false let doc_sample = "n \tSample count [10]" let doc_stop = "n \tStop time [10.0]" let doc_debug = "\tPrint debug information" +let doc_greedy = "\tUse greedy simulation" let opts = [ "-sample", Arg.Set_int sample, doc_sample; "-stop", Arg.Set_float stop, doc_stop; - "-debug", Arg.Set debug, doc_debug + "-debug", Arg.Set Debug.debug, doc_debug; + "-greedy", Arg.Set greedy, doc_greedy; ] let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS]\nOptions are:" @@ -27,6 +30,11 @@ let () = let zsolver = StatefulZ.Functional.zsolve in let solver = Solver.solver_c csolver zsolver in let model = Ball.bouncing_ball () in - let open Sim.LazySim(State.FunctionalSimState) in - run_until model (Solver.solver_from_c solver) !stop (Output.print !sample) + if !greedy then + let open Sim.GreedySim(State.FunctionalSimState) in + run_until model solver !stop (Output.print !sample) + else + let open Sim.LazySim(State.FunctionalSimState) in + run_until model (Solver.solver_from_c solver) !stop (Output.print !sample) + diff --git a/src/bin/output.ml b/src/bin/output.ml index c28e8c1..d977541 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -1,5 +1,6 @@ open Hsim.Types +open Common let print_entry t y = let n = Bigarray.Array1.dim y in @@ -17,5 +18,11 @@ let print samples { start; length; u } = if i > samples then () else let t = float_of_int i *. step in (print_entry (start +. t) (u t); loop (i+1)) in - loop 0 + if length <= 0.0 then + begin Debug.print "D: "; print_entry start (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) diff --git a/src/lib/common/debug.ml b/src/lib/common/debug.ml new file mode 100644 index 0000000..8c2c3c4 --- /dev/null +++ b/src/lib/common/debug.ml @@ -0,0 +1,6 @@ + +let debug = ref false + +let print (s: string) = + if !debug then begin Format.printf "%s" s; flush stdout end else () + diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 00594db..4344693 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -41,7 +41,7 @@ module LazySim (S : SimState) = let mode, stop, now = Continuous, i.length, 0.0 in S.update ms ss (S.set_running ~mode ~input:i ~stop ~now s) else S.set_running ~mode:Continuous s in - Some { start = i.start+. now; length = 0.0; u = fun _ -> o }, s + Some { start = i.start +. now; length = 0.0; u = fun _ -> o }, s | Continuous -> let (h, f, z), ss = solver.step ss stop in let ms = model.cset ms (f h) in @@ -120,31 +120,29 @@ module GreedySim (S : SimState) = let s = S.set_running ~mode ~input:i ~stop ~now s in step (S.update ms ss s) i else step (S.set_running ~mode:Continuous s) i in - let start = i.start +. now in - { start; length = 0.0; u = fun _ -> o }::rest, s + { start = i.start +. now; length = 0.0; u = fun _ -> o }::rest, s | Continuous -> let (h, f, z), ss = solver.step ss stop in (* Copy the state to allow [f] to remain independent from further modifications. *) let ss = solver.copy ss in let ms = model.cset ms (f h) in - let h' = i.start +. h in let fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in let out = { start = i.start +. now; length = h -. now; u = fout } in match z with | None -> if h >= stop then - let s = S.set_running ~mode:Discrete ~now:h' s in + let s = S.set_running ~mode:Discrete ~now:h s in let rest, s = step (S.update ms ss s) i in out::rest, s else - let s = S.set_running ~now:h' s in + let s = S.set_running ~now:h s in let rest, s = step (S.update ms ss s) i in (match rest with | [] -> [out], s | f::rest -> Utils.compose [out;f] :: rest, s) | Some z -> - let s = S.set_running ~mode:Discrete ~now:h' s in + let s = S.set_running ~mode:Discrete ~now:h s in let ms = model.zset ms z in let rest, s = step (S.update ms ss s) i in out::rest, s in @@ -156,4 +154,16 @@ module GreedySim (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 o, _ = sim.step sim.state input in + List.iter use o + + (** Run the model autonomously until [length], or until the model stops + answering. *) + let run_until model solver length = + run_on model solver { start = 0.0; length; u = fun _ -> () } + end diff --git a/src/lib/hsim/utils.ml b/src/lib/hsim/utils.ml index 8db8a52..e4da67f 100644 --- a/src/lib/hsim/utils.ml +++ b/src/lib/hsim/utils.ml @@ -1,15 +1,24 @@ open Types -(** Offset the [input] function by [now]. *) +(** Offset the [input.u] function by [now]. *) let offset (input : 'a value) (now : time) : time -> 'a = fun t -> input.u ((now -. input.start) +. t) +(** +Concatenate functions. [ +^ ^ +| ---, | ---, +| ___ `--- = | _ `--- +| --' | --' ++--------------> +-------------->] +*) let rec compose = function - | [] -> assert false + | [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [f] -> f - | { start=sl; u=ul; _ } :: l -> + | { start; u; _ } :: l -> let { start=sr; length=lr; u=ur } = compose l in - let length = sr +. lr -. sl in - { start=sl; length; u=fun t -> if t <= sr then ur t else ul t } + 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) }