hms/hsim
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat: remove greedy sim, acceleration based on continuity, and composition

Browse source
This commit is contained in:
Henri Saudubray 2025-06-04 16:26:37 +02:00
parent 1a4f950324
commit 65918ab59b
Signed by: hms
GPG key ID: 7065F57ED8856128
5 changed files with 195 additions and 177 deletions
Download patch file Download diff file Expand all files Collapse all files

19
src/bin/main.ml
View file

@ -7,7 +7,7 @@ open Types
let sample = ref 10 let sample = ref 10
let stop = ref 30.0 let stop = ref 30.0
let greedy = ref false let accel = ref false
let inplace = ref false let inplace = ref false
let sundials = ref false let sundials = ref false
let steps = ref 1 let steps = ref 1
@ -26,8 +26,8 @@ let opts = [
"-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)"; "-sample", Arg.Int (gt0i sample), "n \tSample count (default=10)";
"-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)"; "-stop", Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
"-debug", Arg.Set Debug.debug, "\tPrint debug information"; "-debug", Arg.Set Debug.debug, "\tPrint debug information";
"-greedy", Arg.Set greedy, "\tUse greedy simulation"; "-accelerate", Arg.Set accel, "\tConcatenate continuous functions";
"-sundials", Arg.Set sundials, "\tUse sundials (not compatible with greedy)"; "-sundials", Arg.Set sundials, "\tUse sundials (does not support acceleration)";
"-inplace", Arg.Set inplace, "\tUse imperative solvers"; "-inplace", Arg.Set inplace, "\tUse imperative solvers";
"-steps", Arg.Int (gt0i steps), "n \tSplit into [n] steps (default=1)"; "-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 = let sim =
if !sundials then if !sundials then
if !greedy then
(Format.eprintf "Sundials does not support greedy simulation\n"; exit 2)
else
let open StatefulSundials in let open StatefulSundials in
let c = if !inplace then InPlace.csolve else Functional.csolve in let c = if !inplace then InPlace.csolve else Functional.csolve in
let s = Solver.solver c (d_of_dc z) 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 Sim.Sim(val st) in
run_until_n (Output.print !sample (run m s))
else else
let open StatefulRK45 in let open StatefulRK45 in
let c = if !inplace then InPlace.csolve else Functional.csolve in let c = if !inplace then InPlace.csolve else Functional.csolve in
let s = Solver.solver_c c z in let s = Solver.solver_c c z in
if !greedy then let open Sim.GreedySim(val st) in run_until_n m s let open Sim.Sim(val st) in
else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s) 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

14
src/bin/output.ml
View file

@ -1,8 +1,9 @@
open Hsim.Types open Hsim.Types
open Hsim.Utils
open Common open Common
let print_entry t y = let print_entry y t =
let n = Bigarray.Array1.dim y in let n = Bigarray.Array1.dim y in
let rec loop i = let rec loop i =
if i = n then () if i = n then ()
@ -12,17 +13,20 @@ let print_entry t y =
Format.printf "\n"; Format.printf "\n";
flush stdout flush stdout
let print samples { h; u } = let print_sample samples ({ h; u; _ }, now) =
let step = h /. (float_of_int samples) in let step = h /. (float_of_int samples) in
let rec loop i = let rec loop i =
if i > samples then () 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 else let t = float_of_int i *. step in
(print_entry t (u t); loop (i+1)) in (print_entry (u t) (now +. t); loop (i+1)) in
if h <= 0.0 then begin Debug.print "D: "; print_entry 0.0 (u 0.0) end if h <= 0.0 then begin Debug.print "D: "; print_entry (u 0.0) now end
else begin Debug.print "C: "; loop 0 end else begin Debug.print "C: "; loop 0 end
let print_limits { h; _ } = let print_limits { h; _ } =
if h <= 0.0 then Format.printf "D: % .10e\n" 0.0 if h <= 0.0 then Format.printf "D: % .10e\n" 0.0
else Format.printf "C: % .10e\t% .10e\n" 0.0 h 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, ((), ())) }

256
src/lib/hsim/sim.ml
View file

@ -3,60 +3,70 @@ open Types
open Solver open Solver
open State open State
module LazySim (S : SimState) = module Sim (S : SimState) =
struct struct
include S 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 let run
(HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode) (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
(DNode solver : ('y, 'yder, 'zin, 'zout) solver) (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
: ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
= let state = get_init model.state solver.state in = 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 step_continuous s =
let ms, ss = get_mstate s, get_sstate s in step_continuous s solver.step model.cset model.fout model.zset in
match i, is_running s with
| Some i, _ -> let step s = function
| Some i ->
let mode, now, stop = Discrete, 0.0, i.h in let mode, now, stop = Discrete, 0.0, i.h in
None, set_running ~mode ~input:i ~now ~stop s step_discrete (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 -> | None ->
let s = if h >= stop if is_running s then match get_mode s with
then set_running ~mode:Discrete ~now:h s | Discrete -> step_discrete s
else set_running ~now:h s in | Continuous -> step_continuous s
update ms ss s else None, s in
| Some z ->
let s = set_running ~mode:Discrete ~now:h s in
update (model.zset ms z) ss s in
Some out, s in
let reset (pm, ps) s = let reset (pm, ps) s =
let ms = model.reset pm (get_mstate s) in let ms = model.reset pm (get_mstate s) in
@ -65,128 +75,70 @@ module LazySim (S : SimState) =
DNode { state; step; reset } DNode { state; step; reset }
(** Run the model on the given input until the end of the input or until the let accelerate
model stops answering. *) (HNode m : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
let run_on model solver input use = (DNodeC s : ('y, 'yder, 'zin, 'zout) solver_c)
let DNode sim = run model solver in : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
let out = sim.step sim.state (Some input) in = let state = get_init m.state s.state in
let state = match out with None, s -> s | _ -> assert false in let step_discrete st =
let rec loop state = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize
let o, state = sim.step state None in m.jump s.reset in
match o with None -> () | Some o -> use o; loop state in let step_continuous st =
loop state step_continuous st s.step m.cset m.fout m.zset in
(** Run the model on multiple inputs. *) let rec step st = function
let run_on_n model solver inputs use = | Some i ->
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
(** 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 let mode, now, stop = Discrete, 0.0, i.h in
step (set_running ~mode ~input:i ~now ~stop s) i step_discrete (set_running ~mode ~input:i ~now ~stop st)
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 -> | None ->
if h >= stop then if is_running st then match get_mode st with
let s = set_running ~mode:Discrete ~now:h s in | Discrete -> step_discrete st
let rest, s = step (update ms ss s) i in | Continuous ->
out::rest, s let o, st = step_continuous st in
else match o with
let s = set_running ~now:h s in | None -> None, st
let rest, s = step (update ms ss s) i in | Some { c=Discontinuous; _ } -> o, st
(match rest with | Some ({ c=Continuous; _ } as o) ->
| [] -> [out], s let o', st = step st None in
| f::rest -> Utils.concat [out;f] :: rest, s) match o' with
| Some z -> | None -> assert false
let s = set_running ~mode:Discrete ~now:h s in | Some o' -> Some (Utils.concat [o;o']), st
let ms = model.zset ms z in else None, st in
let rest, s = step (update ms ss s) i in
out::rest, s in
let reset (mp, sp) s = let reset (pm, ps) st =
let ms = model.reset mp (get_mstate s) in let ms = m.reset pm (get_mstate st) in
let ss = solver.reset sp (get_sstate s) in let ss = s.reset ps (get_sstate st) in
update ms ss (set_idle s) in update ms ss (set_idle st) in
DNode { state; step; reset } DNode { state; step; reset }
(** Run the model on the given input until the end of the input or until the (** Run the model on the given input until the end of the input or until the
model stops answering. *) model stops answering. *)
let run_on model solver input use = let run_on (DNode n) input use =
let DNode sim = run model solver in let out = n.step n.state (Some input) in
let o, _ = sim.step sim.state input in let state = match out with None, s -> s | _ -> assert false in
List.iter use o 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. *) (** Run the model on multiple inputs. *)
let run_on_n model solver inputs use = let run_on_n (DNode n) inputs use =
let DNode sim = run model solver in ignore @@ List.fold_left (fun state i ->
let o, _ = List.fold_left (fun (acc, state) i -> let o, state = n.step state (Some i) in
let o, state = sim.step state i in begin match o with None -> () | Some o -> use o end;
o::acc, state) ([], sim.state) inputs in let rec loop state =
List.iter use (List.concat (List.rev o)) 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 (** Run the model autonomously until [h], or until the model stops
answering. *) answering. *)
let run_until model solver h = let run_until n h = run_on n { h; c=Discontinuous; u = fun _ -> () }
run_on model solver { h; 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 end

4
src/lib/hsim/types.ml
View file

@ -1,10 +1,12 @@
type time = float type time = float
type continuity = Continuous | Discontinuous
(** Input and output values are functions defined on intervals. *) (** Input and output values are functions defined on intervals. *)
type 'a value = type 'a value =
{ h : time; { 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 (** A time signal is a sequence of possibly absent α-values
[{ h; u }] where: [{ h; u }] where:

67
src/lib/hsim/utils.ml
View file

@ -16,7 +16,68 @@ Concatenate functions. [
let rec concat = function let rec concat = function
| [] -> raise (Invalid_argument "Cannot concatenate an empty value list") | [] -> raise (Invalid_argument "Cannot concatenate an empty value list")
| [f] -> f | [f] -> f
| { u; h } :: l -> | { c=Discontinuous; _ } :: _ ->
let { h=hr; u=ur } = concat l in raise (Invalid_argument "Cannot concatenate discontinuous functions")
{ h=h+.hr; u=fun t -> if t <= h then u t else ur (t -. h) } | { 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 }