feat: runtime as library

src/lib/hsim/sim.ml

@@ -9,7 +9,9 @@ module Sim (S : SimState) =
 
     let step_discrete
       s step hor fder fzer cget zset csize zsize jump reset reinit
-    = let ms, ss, zin = get_mstate s, get_sstate s, get_zin s in
+    = let ms, ss = get_mstate s, get_sstate s in
+      let zin, last = get_zin s, get_last s in
+      (match last with Some { h; u; _ } -> ignore (u h) | None -> ());
       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 now (i.u now) in
@@ -28,10 +30,12 @@ module Sim (S : SimState) =
           let mode, stop, now = Continuous, i.h, 0.0 in
           update ms ss (set_running ~mode ~input ~stop ~now s)
         end else set_running ~mode:Continuous s in
-      Utils.dot o, (set_zin None s)
+      let o = Utils.dot o in
+      o, (set_last (Some o) (set_zin None s))
 
     let step_continuous s step cset fout hor =
-      let ms, ss = get_mstate s, get_sstate s in
+      let ms, ss, last = get_mstate s, get_sstate s, get_last s in
+      (match last with None -> () | Some { h; u; _ } -> ignore (u h));
       let i, now, stop = get_input s, get_now s, get_stop s in
       let stop = min stop (hor ms) in
       let (h, f, z), ss = step ss (min stop (hor ms)) in
@@ -47,7 +51,8 @@ module Sim (S : SimState) =
           else set_running ~now:h s, Continuous
       | Some _ -> set_running ~mode:Discrete ~now:h s, Discontinuous in
       let h = h -. now in
-      { h; u=fout; c }, update ms ss (set_zin z s), { h; c; u=fms }
+      let o = { h; u=fout; c } in
+      o, update ms ss (set_last (Some o) (set_zin z s)), { h; c; u=fms }
 
     (** Simulation of a model with any solver. *)
     let run
@@ -171,7 +176,7 @@ module Sim (S : SimState) =
         model stops answering. *)
     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 state = match out with None, s -> s | Some o, s -> use o; s in
       let rec loop state =
         let o, state = n.step state None in
         match o with None -> () | Some o -> use o; loop state in

src/lib/hsim/state.ml

@@ -15,65 +15,71 @@ 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, 'zin) state
+    type ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Get the model state. *)
-    val get_mstate : ('a, 'ms, 'ss, 'zin) state -> 'ms
+    val get_mstate : ('a, 'b, 'ms, 'ss, 'zin) state -> 'ms
 
     (** Get the solver state. *)
-    val get_sstate : ('a, 'ms, 'ss, 'zin) state -> 'ss
+    val get_sstate : ('a, 'b, 'ms, 'ss, 'zin) state -> 'ss
 
     (** Get the last zero-crossing value. *)
-    val get_zin : ('a, 'ms, 'ss, 'zin) state -> 'zin option
+    val get_zin : ('a, 'b, 'ms, 'ss, 'zin) state -> 'zin option
+
+    (** Get the last produced value. *)
+    val get_last : ('a, 'b, 'ms, 'ss, 'zin) state -> 'b signal
 
     (** Get the current step mode.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_mode : ('a, 'ms, 'ss, 'zin) state -> mode
+    val get_mode : ('a, 'b, 'ms, 'ss, 'zin) state -> mode
 
     (** Get the current input.
         ⚠ Should only be called when running (see [is_running]). *)
-    val get_input : ('a, 'ms, 'ss, 'zin) state -> 'a value
+    val get_input : ('a, 'b, '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, 'zin) state -> time
+    val get_now : ('a, 'b, '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, 'zin) state -> time
+    val get_stop : ('a, 'b, 'ms, 'ss, 'zin) state -> time
 
     (** Build an initial state. *)
-    val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state
+    val get_init : 'ms -> 'ss -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Is the simulation running or idle ? *)
-    val is_running : ('a, 'ms, 'ss, 'zin) state -> bool
+    val is_running : ('a, 'b, 'ms, 'ss, 'zin) state -> bool
 
     (** Update the model state. *)
-    val set_mstate : 'ms -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val set_mstate : 'ms -> ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Update the solver state. *)
-    val set_sstate : 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val set_sstate : 'ss -> ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Update the zero-crossing value. *)
-    val set_zin : 'zin option -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val set_zin : 'zin option -> ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
+
+    (** Update the last produced value. *)
+    val set_last : 'b signal -> ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Update both the solver and model states. *)
-    val update : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val update : 'ms -> 'ss -> ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Update the status to running. *)
     val set_running :
       ?mode:mode -> ?input:'a value -> ?now:time -> ?stop:time ->
-        ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+        ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
 
     (** Update the status to idle. *)
-    val set_idle : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val set_idle : ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
   end
 
 module type SimStateCopy =
   sig
     include SimState
 
-    val copy : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
+    val copy : ('a, 'b, 'ms, 'ss, 'zin) state -> ('a, 'b, 'ms, 'ss, 'zin) state
   end
 
 module FunctionalSimState : SimState =
@@ -94,17 +100,19 @@ module FunctionalSimState : SimState =
 
     (** Internal state of the simulation node: model state, solver state and
         current simulation status. *)
-    type ('a, 'ms, 'ss, 'zin) state =
+    type ('a, 'b, 'ms, 'ss, 'zin) state =
       { status : 'a status;                  (** Current simulation status.   *)
         mstate : 'ms;                        (** Model state.                 *)
         sstate : 'ss;                        (** Solver state.                *)
-        zin    : 'zin option; }              (** Last zero-crossing vector    *)
+        zin    : 'zin option;                (** Last zero-crossing vector.   *)
+        last   : 'b signal }                 (** Last produced value.         *)
 
     exception Not_running
 
     let get_mstate state = state.mstate
     let get_sstate state = state.sstate
     let get_zin state = state.zin
+    let get_last state = state.last
 
     let is_running state =
       match state.status with Running _ -> true | Idle -> false
@@ -129,6 +137,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 set_last last state = { state with last }
 
     let update mstate sstate state = { state with mstate; sstate }
 
@@ -141,7 +150,8 @@ 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; zin = None }
+    let get_init mstate sstate =
+      { status=Idle; mstate; sstate; zin=None; last=None }
   end
 
 module InPlaceSimState : SimState =
@@ -155,17 +165,19 @@ module InPlaceSimState : SimState =
             mutable stop  : time;
           } -> 'a status
 
-    type ('a, 'ms, 'ss, 'zin) state =
+    type ('a, 'b, 'ms, 'ss, 'zin) state =
       { mutable status : 'a status;
         mutable mstate : 'ms;
         mutable sstate : 'ss;
-        mutable zin    : 'zin option }
+        mutable zin    : 'zin option;
+        mutable last   : 'b signal }
 
     exception Not_running
 
     let get_mstate state = state.mstate
     let get_sstate state = state.sstate
     let get_zin state = state.zin
+    let get_last state = state.last
 
     let is_running state =
       match state.status with Running _ -> true | Idle -> false
@@ -191,6 +203,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 set_last last state = state.last <- last; state
 
     let update mstate sstate state =
       state.mstate <- mstate; state.sstate <- sstate; state
@@ -204,6 +217,7 @@ 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; zin=None }
+    let get_init mstate sstate =
+      { status=Idle; mstate; sstate; zin=None; last=None }
   end
 

src/lib/hsim/utils.ml

@@ -27,6 +27,9 @@ let sample { h; u; _ } n =
     (t, u t) :: step (i+1) in
   if h <= 0.0 then [(0.0, u 0.0)] else step 0
 
+let sample_tracked (o, t) n =
+  List.map (fun (h, v) -> h +. t, v) @@ sample o n
+
 (** Compose two nodes together. *)
 let compose (DNode m) (DNode n) =
   let state = m.state, n.state in
@@ -81,11 +84,20 @@ let map f =
 
 let ignore _ n =
   let state = () in
-  let step () = function
-    | None -> None, ()
-    | Some _ -> Some (), () 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
+  let DNode n = compose n @@ DNode { state; step; reset } in
   DNode { n with reset=fun p -> n.reset (p, ()) }
 
+let do_and_reset (DNode m) (DNode n) f =
+  let state = m.state, n.state in
+  let step (ms, ns) i =
+    let o, ms = m.step ms i in
+    let v, ns = n.step ns o in
+    begin match v with Some v -> f v; | None -> () end;
+    begin match o with Some { h; u; _ } -> Stdlib.ignore (u h) | None -> () end;
+    v, (ms, ns) in
+  let reset (ms, ns) (mp, np) =
+    m.reset ms mp, n.reset ns np in
+  DNode { state; step; reset }
+

src/lib/std/lift.ml

@@ -26,37 +26,25 @@ let lift f =
 
   (* the function that compute the derivatives *)
   let fder { state; time; _ } offset input y =
-    cstate.major <- false;
-    cstate.zinvec <- no_roots_in;
-    cstate.zoutvec <- no_roots_out;
-    cstate.cvec <- y;
-    cstate.dvec <- ignore_der;
-    cstate.cindex <- 0;
-    cstate.zindex <- 0;
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0;
     ignore (f_step state (time +. offset, input));
     cstate.dvec in
 
   (* the function that compute the zero-crossings *)
   let fzer { state; time; _ } offset input y =
-    cstate.major <- false;
-    cstate.zinvec <- no_roots_in;
-    cstate.dvec <- ignore_der;
-    cstate.zoutvec <- no_roots_out;
-    cstate.cvec <- y;
-    cstate.cindex <- 0;
-    cstate.zindex <- 0;
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0;
     ignore (f_step state (time +. offset, input));
     cstate.zoutvec in
 
   (* the function which compute the output during integration *)
   let fout { state; time; _ } offset input y =
-    cstate.major <- false;
-    cstate.zoutvec <- no_roots_out;
-    cstate.dvec <- ignore_der;
-    cstate.zinvec <- no_roots_in;
-    cstate.cvec <- y;
-    cstate.cindex <- 0;
-    cstate.zindex <- 0;
+    cstate.major <- false; cstate.cvec <- y; cstate.dvec <- ignore_der;
+    cstate.zinvec <- no_roots_in; cstate.zoutvec <- no_roots_out;
+    cstate.cindex <- 0; cstate.zindex <- 0;
     f_step state (time +. offset, input) in
 
   (* the function which compute a discrete step *)
@@ -77,7 +65,6 @@ let lift f =
 
   (* horizon *)
   let horizon { time; _ } =
-    (* Printf.printf "\tCalling horizon :: cstate.horizon=%.10e\ttime=%.10e\n" cstate.horizon time; *)
     cstate.horizon -. time in
 
   let jump _ = true in

src/lib/std/output.ml

@@ -0,0 +1,53 @@
+
+open Hsim.Types
+open Hsim.Utils
+
+let print_entry y t =
+  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;
+  loop 0;
+  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 =
+    if i > samples then ()
+    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 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
+  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, ((), ())) }
+
+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, ((), ())) }

src/lib/std/runtime.ml

@@ -0,0 +1,31 @@
+
+open Hsim.Types
+
+let sample = ref 0
+let stop   = ref 10.0
+
+let options = [
+  "-sample", Arg.Set_int sample, "\tSampling frequency (default=0)";
+  "-stop", Arg.Set_float stop, "\tStop time (default=10.0)";
+  "-debug", Arg.Set Common.Debug.debug, "\tShow debug information";
+]
+
+let arg s =
+  Format.eprintf "Unexpected argument: %s\n" s; exit 1
+
+let usage = ""
+
+let go
+  (input : time -> 'a)
+  (model : Ztypes.cstate -> (time * 'a, 'b) Ztypes.node)
+  (output : (time * 'b) -> unit)
+= Arg.parse options arg usage;
+  let input = { h=(!stop); c=Discontinuous; u=input } in
+  let output o = List.iter output @@ Hsim.Utils.sample_tracked o !sample in
+  let model = Lift.lift model in
+  let open Hsim in
+  let solver = Solver.solver_c Solvers.StatefulRK45.InPlace.csolve
+                               Solvers.StatefulZ.InPlace.zsolve in
+  let open Sim.Sim(State.InPlaceSimState) in
+  let sim = Hsim.Utils.(compose (run model (d_of_dc solver)) track) in
+  run_on sim input output