From 5bce9e5b01683410ae92d484ab6b092dba838e07 Mon Sep 17 00:00:00 2001 From: Henri Saudubray Date: Mon, 28 Apr 2025 15:13:15 +0200 Subject: [PATCH] feat: correct greedy/lazy and inplace/functional, split into multiple inputs --- README.md | 2 +- doc/notes.typ | 4 +-- src/bin/main.ml | 53 ++++++++++++++++++++++++------------ src/bin/output.ml | 14 +++++----- src/lib/common/debug.ml | 5 ++-- src/lib/hsim/sim.ml | 44 +++++++++++++++++++++++++++++- src/lib/hsim/solver.ml | 4 +-- src/lib/hsim/state.ml | 16 +++-------- src/lib/hsim/statefulRK45.ml | 2 +- src/lib/hsim/statefulZ.ml | 2 +- src/lib/solvers/illinois.ml | 24 ++++++++-------- src/lib/solvers/odexx.ml | 12 ++++---- 12 files changed, 117 insertions(+), 65 deletions(-) diff --git a/README.md b/README.md index f94eac7..ff14ee9 100644 --- a/README.md +++ b/README.md @@ -2,5 +2,5 @@ A hybrid system simulation semantics. -Implemented with, and heavily inspired by work from, Marc Pouzet and Timothy +Implemented with, and heavily inspired by work from, Marc Pouzet and Timothy Bourke (PARKAS, Inria, École Normale Supérieure). diff --git a/doc/notes.typ b/doc/notes.typ index 0d40e28..9a3ba54 100644 --- a/doc/notes.typ +++ b/doc/notes.typ @@ -297,7 +297,7 @@ Two possible options for the simulation reset: `fder : 'a -> 'y -> 'yder` and `fzer : 'a -> 'y -> 'zout`, and we thus need a way to obtain a value of type `'a`. This is usually done through `input.u : time -> 'a`, but we have no input available during the reset, which - makes this impossible. We thus need reset parameters for both the model and + makes this impossible. We thus need reset parameters for both the model and solver. === Steps @@ -375,7 +375,7 @@ let rec step s i = 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; + 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 diff --git a/src/bin/main.ml b/src/bin/main.ml index 8487520..3984e49 100644 --- a/src/bin/main.ml +++ b/src/bin/main.ml @@ -3,20 +3,29 @@ open Hsim open Examples open Common -let sample = ref 10 -let stop = ref 30.0 -let greedy = ref false +let sample = ref 10 +let stop = ref 30.0 +let greedy = ref false +let inplace = ref false +let steps = ref 1 -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 gt0i v i = v := if i <= 0 then 1 else i +let gt0f v f = v := if f <= 0.0 then 1.0 else f + +let doc_sample = "n \tSample count (default=10)" +let doc_stop = "n \tStop time (default=10.0)" +let doc_debug = "\tPrint debug information" +let doc_greedy = "\tUse greedy simulation" +let doc_inplace = "\tUse greedy simulation" +let doc_steps = "n \tSplit simulation into [n] steps (default=1)" let opts = [ - "-sample", Arg.Set_int sample, doc_sample; - "-stop", Arg.Set_float stop, doc_stop; - "-debug", Arg.Set Debug.debug, doc_debug; - "-greedy", Arg.Set greedy, doc_greedy; + "-sample", Arg.Int (gt0i sample), doc_sample; + "-stop", Arg.Float (gt0f stop), doc_stop; + "-debug", Arg.Set Debug.debug, doc_debug; + "-greedy", Arg.Set greedy, doc_greedy; + "-inplace", Arg.Set inplace, doc_inplace; + "-steps", Arg.Int (gt0i steps), doc_steps; ] let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS]\nOptions are:" @@ -30,11 +39,19 @@ let () = let zsolver = StatefulZ.Functional.zsolve in let solver = Solver.solver_c csolver zsolver in let model = Ball.bouncing_ball () in - if !greedy then - let open Sim.GreedySim(State.FunctionalSimState) in - run_until model solver !stop (Output.print !sample) + if !inplace then + let module S = State.InPlaceSimState in + if !greedy then + let open Sim.GreedySim(S) in + run_until_multiple model solver !stop !steps (Output.print !sample) + else + let open Sim.LazySim(S) in + run_until_multiple model (Solver.solver_from_c solver) !stop !steps (Output.print !sample) else - let open Sim.LazySim(State.FunctionalSimState) in - run_until model (Solver.solver_from_c solver) !stop (Output.print !sample) - - + let module S = State.FunctionalSimState in + if !greedy then + let open Sim.GreedySim(S) in + run_until_multiple model solver !stop !steps (Output.print !sample) + else + let open Sim.LazySim(S) in + run_until_multiple model (Solver.solver_from_c solver) !stop !steps (Output.print !sample) diff --git a/src/bin/output.ml b/src/bin/output.ml index d977541..715a739 100644 --- a/src/bin/output.ml +++ b/src/bin/output.ml @@ -6,23 +6,23 @@ let print_entry t y = let n = Bigarray.Array1.dim y in let rec loop i = if i = n then () - else (Printf.printf "\t% .10e" y.{i}; loop (i+1)) in - Printf.printf "% .10e" t; + else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in + Format.printf "% .10e" t; loop 0; - Printf.printf "\n"; + Format.printf "\n"; flush stdout let print samples { start; length; u } = let step = length /. (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 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 - else - begin Debug.print "C: "; loop 0 end + 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 index 8c2c3c4..fdd37e7 100644 --- a/src/lib/common/debug.ml +++ b/src/lib/common/debug.ml @@ -1,6 +1,7 @@ let debug = ref false -let print (s: string) = - if !debug then begin Format.printf "%s" s; flush stdout end else () +let fmt () = if !debug then Format.std_formatter + else Format.make_formatter (fun _ _ _ -> ()) (fun () -> ()) +let print = Format.fprintf (fmt ()) "%s" diff --git a/src/lib/hsim/sim.ml b/src/lib/hsim/sim.ml index 4344693..5fd0c7a 100644 --- a/src/lib/hsim/sim.ml +++ b/src/lib/hsim/sim.ml @@ -78,11 +78,34 @@ module LazySim (S : SimState) = | Some o -> use o; loop (DNode { s with state }) in loop (DNode { sim with state }) + (** Run the model on multiple inputs. *) + let run_on_multiple model solver inputs use = + ignore @@ List.fold_left (fun (DNode sim) i -> + Common.Debug.print + (Format.sprintf "New input: %.10e\t%.10e\n" i.start i.length); + let state = match sim.step sim.state (Some i) with + | None, s -> s | _ -> assert false in + let rec loop (DNode s) = + let o, state = s.step s.state None in + match o with + | None -> DNode { s with state } + | Some o -> use o; loop (DNode { s with state }) in + loop (DNode { sim with state })) (run model solver) inputs + (** 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 _ -> () } + let run_until_multiple 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 + let length = stop -. start in + { start; length; u = fun _ -> () }) in + run_on_multiple model solver inputs + end module GreedySim (S : SimState) = @@ -106,7 +129,7 @@ module GreedySim (S : SimState) = 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 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 @@ -161,9 +184,28 @@ module GreedySim (S : SimState) = let o, _ = sim.step sim.state input in List.iter use o + (** Run the model on the given input list. *) + let run_on_multiple model solver inputs use = + let o, _ = List.fold_left (fun (acc, DNode sim) i -> + Common.Debug.print + (Format.sprintf "new input: %.10e\t%.10e\n" i.start i.length); + let o, state = sim.step sim.state i in + o::acc, DNode { sim with state }) ([], run model solver) inputs in + List.iter use (List.concat (List.rev o)) + (** Run the model autonomously until [length], or until the model stops answering. *) let run_until model solver length = run_on model solver { start = 0.0; length; u = fun _ -> () } + (** Run the model autonomously until [length], split in multiple [steps]. *) + let run_until_multiple 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 + let length = stop -. start in + { start; length; u = fun _ -> () }) in + run_on_multiple model solver inputs + end diff --git a/src/lib/hsim/solver.ml b/src/lib/hsim/solver.ml index 4bafabb..c8be4a7 100644 --- a/src/lib/hsim/solver.ml +++ b/src/lib/hsim/solver.ml @@ -32,7 +32,7 @@ type ('y, 'yder) csolver_c = type ('y, 'zin, 'zout) zsolver = (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode -(** A zero-crossing solver can optionally provide a state copy method, in which +(** A zero-crossing solver can optionally provide a state copy method, in which case greedy simulation is possible. *) type ('y, 'zin, 'zout) zsolver_c = (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode_c @@ -46,7 +46,7 @@ type ('y, 'yder, 'zin, 'zout) solver = time, time * (time -> 'y) * 'zin option) dnode -(** A solver can optionally provide a state copy method, in which case greedy +(** A solver can optionally provide a state copy method, in which case greedy simulation is possible. *) type ('y, 'yder, 'zin, 'zout) solver_c = (('y, 'yder) ivp * ('y, 'zout) zc, diff --git a/src/lib/hsim/state.ml b/src/lib/hsim/state.ml index dd643ec..676e12a 100644 --- a/src/lib/hsim/state.ml +++ b/src/lib/hsim/state.ml @@ -39,8 +39,6 @@ module type SimState = ⚠ Should only be called when running (see [is_running]). *) val get_stop : ('a, 'ms, 'ss) state -> time - val get_t0 : ('a, 'ms, 'ss) state -> time - (** Build an initial state. *) val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss) state @@ -86,7 +84,6 @@ module FunctionalSimState : SimState = input : 'a value; (** Function to integrate. *) now : time; (** Current time of integration. *) stop : time; (** How long until we stop. *) - t0 : time; (** Initial start time. *) } -> 'a status (** Internal state of the simulation node: model state, solver state and @@ -111,15 +108,15 @@ module FunctionalSimState : SimState = | Idle -> begin match mode, input, now, stop with | Some mode, Some input, Some now, Some stop -> - { state with status = Running { mode; input; now; stop; t0 = input.start } } + { state with status = Running { mode; input; now; stop } } | _ -> raise (Invalid_argument "") end - | Running { mode=m; input=i; now=n; stop=s; t0 } -> + | Running { mode=m; input=i; now=n; stop=s } -> let mode = Option.value mode ~default:m in let input = Option.value input ~default:i in let now = Option.value now ~default:n in let stop = Option.value stop ~default:s in - { state with status = Running { mode; input; now; stop; t0 } } + { state with status = Running { mode; input; now; stop } } let set_mstate mstate state = { state with mstate } let set_sstate sstate state = { state with sstate } @@ -134,8 +131,6 @@ module FunctionalSimState : SimState = match s.status with Running r -> r.now | Idle -> raise Not_running let get_stop s = match s.status with Running r -> r.stop | Idle -> raise Not_running - let get_t0 s = - match s.status with Running r -> r.t0 | Idle -> raise Not_running let get_init mstate sstate = { status = Idle; mstate; sstate } end @@ -149,7 +144,6 @@ module InPlaceSimState : SimState = mutable input : 'a value; mutable now : time; mutable stop : time; - mutable t0 : time; } -> 'a status type ('a, 'ms, 'ss) state = @@ -172,7 +166,7 @@ module InPlaceSimState : SimState = | Idle -> begin match mode, input, now, stop with | Some mode, Some input, Some now, Some stop -> - state.status <- Running { mode; input; now; stop; t0 = input.start }; + state.status <- Running { mode; input; now; stop }; state | _ -> raise (Invalid_argument "") end @@ -197,8 +191,6 @@ module InPlaceSimState : SimState = match s.status with Running r -> r.now | Idle -> raise Not_running let get_stop s = match s.status with Running r -> r.stop | Idle -> raise Not_running - let get_t0 s = - match s.status with Running r -> r.t0 | Idle -> raise Not_running let get_init mstate sstate = { status = Idle; mstate; sstate } end diff --git a/src/lib/hsim/statefulRK45.ml b/src/lib/hsim/statefulRK45.ml index 05ec6b4..880a41d 100644 --- a/src/lib/hsim/statefulRK45.ml +++ b/src/lib/hsim/statefulRK45.ml @@ -38,7 +38,7 @@ module Functional = module InPlace = struct - + type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec } diff --git a/src/lib/hsim/statefulZ.ml b/src/lib/hsim/statefulZ.ml index 2f8ed8d..cf8787c 100644 --- a/src/lib/hsim/statefulZ.ml +++ b/src/lib/hsim/statefulZ.ml @@ -11,7 +11,7 @@ module Functional = let zsolve : (Zls.carray, Zls.zarray, Zls.carray) zsolver_c = let state = { state = Illinois.initialize 0 (fun _ _ _ -> ()) (Zls.cmake 0); - vec = Zls.zmake 0 } in + vec = Zls.zmake 0 } in let reset { fzer; init; size } { vec; _ } = let fzer t cvec zout = let zout' = fzer t cvec in Zls.blit zout' zout in { state = Illinois.initialize size fzer init; diff --git a/src/lib/solvers/illinois.ml b/src/lib/solvers/illinois.ml index 5fd2dc9..2919029 100644 --- a/src/lib/solvers/illinois.ml +++ b/src/lib/solvers/illinois.ml @@ -54,7 +54,7 @@ let make_check_root rdir f0 f1 = let check = get_check_root rdir in (fun i -> check f0.{i} f1.{i}) **) - + (* update roots and returns true if there was at least one root *) (* between f0 and f1 for one component of index [i in [0..length f0 - 1]] *) (* update [roots] *) @@ -85,7 +85,7 @@ type zcfn = float -> Zls.carray -> Zls.carray -> unit (* type of a session with the solver *) (* zx = g(t, c) yields the values of system zero-crossing expressions - + f0/t0 are the zero-crossing expression values at the last mesh point f1/t1 are the zero-crossing expression values at the next mesh point @@ -129,16 +129,16 @@ let initialize_only nroots g = { g = g; bothf_valid = false; - + f0 = Zls.cmake nroots; t0 = 0.0; - + f1 = Zls.cmake nroots; t1 = 0.0; - + fta = Zls.cmake nroots; ftb = Zls.cmake nroots; - + calc_zc = get_check_root Up; } @@ -224,8 +224,8 @@ let find ({ g = g; bothf_valid = bothf_valid; (* Searches between (t_left, f_left) and (t_right, f_right) to find the leftmost (t_mid, f_mid): - - | + + | | f_right | | f_mid @@ -323,7 +323,7 @@ let find ({ g = g; bothf_valid = bothf_valid; | SearchLeft -> begin let (t, v, f0', fta', ftb') = seek (t0, f0, fta, t1, None, 1.0, 0) in - + s.t0 <- t; s.f0 <- f0'; s.bothf_valid <- v; @@ -337,12 +337,12 @@ let find ({ g = g; bothf_valid = bothf_valid; (* the main function of this module *) (* locate a root *) let find s (dky, c) roots = find s (dky, c) roots - + (* is there a root? [has_root s: bool] is true is there is a change in sign *) (* for one component [i in [0..length f0 - 1]] beetwen [f0.(i)] and [f1.(i)] *) let has_roots { bothf_valid = bothf_valid; t0; f0; t1; f1; calc_zc = calc_zc } = bothf_valid && (check_interval calc_zc f0 f1 <> SearchRight) - + let takeoff { bothf_valid = bothf_valid; f0; f1 } = bothf_valid && (takeoff f0 f1) @@ -351,7 +351,7 @@ let takeoff { bothf_valid = bothf_valid; f0; f1 } = (* with function [find] when the signal is taking off from [0.0] to a *) (* strictly positive value *) let find_takeoff ({ f0; f1 } as s) roots = - let calc_zc x0 x1 = + let calc_zc x0 x1 = if (x0 = 0.0) && (x1 > 0.0) then 1l else 0l in let b = update_roots calc_zc f0 f1 roots in if b then begin s.t1 <- s.t0; s.f1 <- s.f0; s.ftb <- s.fta end; diff --git a/src/lib/solvers/odexx.ml b/src/lib/solvers/odexx.ml index 0928df5..3eca888 100644 --- a/src/lib/solvers/odexx.ml +++ b/src/lib/solvers/odexx.ml @@ -40,10 +40,10 @@ sig (* {{{ *) size(b) = ns x ns (but only the lower strictly triangular entries) size(e) = ns - size(bi) = ns x po + size(bi) = ns x po (where po is the order of the interpolating polynomial) *) - + end (* }}} *) @@ -322,7 +322,7 @@ struct (* {{{1 *) s.time <- nextt; s.h <- nexth; s.time - + let get_dky { last_time = t; time = t'; h = h; yold = y; k = k } yi ti kd = if kd > 0 then @@ -334,7 +334,7 @@ struct (* {{{1 *) failwith (Printf.sprintf "get_dky: requested time %.24e is out of range\n\ [%.24e,...,%.24e]" - ti t t'); + ti t t'); let h = t' -. t in let th = (ti -. t) /. h in @@ -358,9 +358,9 @@ struct (* {{{1 *) let copy ({ last_time; time; h; yold; k } as s) = { s with last_time; time; h; yold = Zls.copy yold; k = Zls.copy_matrix k } - let blit { last_time = l1; time = t1; h = h1; yold = yhold1; k = k1 } + let blit { last_time = l1; time = t1; h = h1; yold = yhold1; k = k1 } ({ last_time; time; h; yold; k } as s2) = - s2.last_time <- l1; s2.time <- t1; + s2.last_time <- l1; s2.time <- t1; Zls.blit yhold1 yold; Zls.blit_matrix k1 k end (* }}} *)