feat: correct greedy/lazy and inplace/functional, split into multiple inputs

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).

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

src/bin/main.ml

@@ -3,20 +3,29 @@ open Hsim
 open Examples
 open Common
 
-let sample = ref 10
+let sample  = ref 10
-let stop   = ref 30.0
+let stop    = ref 30.0
-let greedy = ref false
+let greedy  = ref false
+let inplace = ref false
+let steps   = ref 1
 
-let doc_sample = "n \tSample count [10]"
+let gt0i v i = v := if i <= 0   then 1   else i
-let doc_stop   = "n \tStop time [10.0]"
+let gt0f v f = v := if f <= 0.0 then 1.0 else f
-let doc_debug  = "\tPrint debug information"
+
-let doc_greedy = "\tUse greedy simulation"
+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;
+  "-sample",  Arg.Int (gt0i sample), doc_sample;
-  "-stop", Arg.Set_float stop, doc_stop;
+  "-stop",    Arg.Float (gt0f stop), doc_stop;
-  "-debug", Arg.Set Debug.debug, doc_debug;
+  "-debug",   Arg.Set Debug.debug,   doc_debug;
-  "-greedy", Arg.Set greedy, doc_greedy;
+  "-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
+  if !inplace then
-    let open Sim.GreedySim(State.FunctionalSimState) in
+    let module S = State.InPlaceSimState in
-    run_until model solver !stop (Output.print !sample)
+    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
+    let module S = State.FunctionalSimState in
-    run_until model (Solver.solver_from_c solver) !stop (Output.print !sample)
+    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)

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
+    else (Format.printf "\t% .10e" y.{i}; loop (i+1)) in
-  Printf.printf "% .10e" t;
+  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
+  if length <= 0.0 then begin Debug.print "D: "; print_entry start (u 0.0) end
-    begin Debug.print "D: "; print_entry start (u 0.0) end
+  else begin Debug.print "C: "; loop 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)
+

src/lib/common/debug.ml

@@ -1,6 +1,7 @@
 
 let debug = ref false
 
-let print (s: string) =
+let fmt () = if !debug then Format.std_formatter
-  if !debug then begin Format.printf "%s" s; flush stdout end else ()
+             else Format.make_formatter (fun _ _ _ -> ()) (fun () -> ())
 
+let print = Format.fprintf (fmt ()) "%s"

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

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,

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

src/lib/hsim/statefulRK45.ml

@@ -38,7 +38,7 @@ module Functional =
 
 module InPlace =
   struct
-    
+
     type ('state, 'vec) state =
       { mutable state: 'state; mutable vec : 'vec }
 

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;

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;

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 (* }}} *)