feat: greedy simulation

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

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

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)

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 ()
+

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 = i.start +. now; length = 0.0; u = fun _ -> o }::rest, s
-            { start; 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

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
+      let sw = sr -. start in
-      { start=sl; length; u=fun t -> if t <= sr then ur t else ul t }
+      let length = sw +. lr  in
+      { start; length; u=fun t -> if t < sw then u t else ur (t -. sw) }