feat (solver/rk45): state copy is handled by simulation

src/bin/main.ml

@@ -1,5 +1,6 @@
 
 open Hsim
+open Types
 open Examples
 open Common
 
@@ -12,46 +13,31 @@ let steps   = ref 1
 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.Int (gt0i sample), doc_sample;
+  "-sample",  Arg.Int (gt0i sample), "n \tSample count (default=10)";
-  "-stop",    Arg.Float (gt0f stop), doc_stop;
+  "-stop",    Arg.Float (gt0f stop), "n \tStop time (default=10.0)";
-  "-debug",   Arg.Set Debug.debug,   doc_debug;
+  "-debug",   Arg.Set Debug.debug,   "\tPrint debug information";
-  "-greedy",  Arg.Set greedy,        doc_greedy;
+  "-greedy",  Arg.Set greedy,        "\tUse greedy simulation";
-  "-inplace", Arg.Set inplace,       doc_inplace;
+  "-inplace", Arg.Set inplace,       "\tUse greedy simulation";
-  "-steps",   Arg.Int (gt0i steps),  doc_steps;
+  "-steps",   Arg.Int (gt0i steps),  "n \tSplit into [n] steps (default=1)";
 ]
 
 let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS]\nOptions are:"
 
-let () =
+let () = try Arg.parse (Arg.align opts) (fun _ -> ()) errmsg with _ -> exit 2
-  try Arg.parse (Arg.align opts) (fun _ -> ()) errmsg
+
-  with _ -> exit 2
+let output = Output.print !sample
+
+let c = StatefulRK45.(if !inplace then InPlace.csolve else Functional.csolve)
+let z = StatefulZ.(Functional.zsolve)
+let s = Solver.solver_c c z
+let m = Ball.bouncing_ball ()
+
+let state = if !inplace then (module State.InPlaceSimState    : State.SimState)
+                        else (module State.FunctionalSimState : State.SimState)
+let sim = if !greedy
+  then let open Sim.GreedySim(val state) in run_until_n m s
+  else let open Sim.LazySim(val state) in run_until_n m (d_of_dc s)
+
+let () = sim !stop !steps output
 
-let () =
-  let csolver = StatefulRK45.Functional.csolve in
-  let zsolver = StatefulZ.Functional.zsolve in
-  let solver = Solver.solver_c csolver zsolver in
-  let model = Ball.bouncing_ball () in
-  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 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)

src/lib/hsim/dune

@@ -1,4 +1,3 @@
 (library
  (name hsim)
- (libraries common solvers)
+ (libraries common solvers))
- (modules_without_implementation types))

src/lib/hsim/sim.ml

@@ -6,7 +6,7 @@ open State
 module LazySim (S : SimState) =
   struct
 
-    (** "Lazy" simulation of a model with an appropriate solver. *)
+    (** "Lazy" simulation of a model with any solver. *)
     let run
       (HNode model : ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
       (DNode solver : ('y, 'yder, 'zin, 'zout) solver)
@@ -63,6 +63,7 @@ module LazySim (S : SimState) =
         let ms = model.reset pm (S.get_mstate s) in
         let ss = solver.reset ps (S.get_sstate s) in
         S.update ms ss (S.set_idle s) in
+
       DNode { state; step; reset }
 
     (** Run the model on the given input until the end of the input or until the
@@ -79,7 +80,7 @@ module LazySim (S : SimState) =
       loop (DNode { sim with state })
 
     (** Run the model on multiple inputs. *)
-    let run_on_multiple model solver inputs use =
+    let run_on_n 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);
@@ -97,20 +98,19 @@ module LazySim (S : SimState) =
     let run_until model solver length =
       run_on model solver { start = 0.0; length; u = fun _ -> () }
 
-    let run_until_multiple model solver length steps =
+    (** Run the model autonomously until [length], split in multiple [steps]. *)
+    let run_until_n 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 = stop -. start; u = fun _ -> () }) in
-        { start; length; u = fun _ -> () }) in
+      run_on_n model solver inputs
-      run_on_multiple model solver inputs
-
   end
 
 module GreedySim (S : SimState) =
   struct
-
+    (** "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)
@@ -146,8 +146,6 @@ module GreedySim (S : SimState) =
             { start = i.start +. now; 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 fout t = model.fout ms (i.u (now +. t)) (f (now +. t)) in
@@ -184,8 +182,8 @@ module GreedySim (S : SimState) =
       let o, _ = sim.step sim.state input in
       List.iter use o
 
-    (** Run the model on the given input list. *)
+    (** Run the model on multiple inputs. *)
-    let run_on_multiple model solver inputs use =
+    let run_on_n 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);
@@ -199,13 +197,11 @@ module GreedySim (S : SimState) =
       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 run_until_n 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 = stop -. start; u = fun _ -> () }) in
-        { start; length; u = fun _ -> () }) in
+      run_on_n model solver inputs
-      run_on_multiple model solver inputs
-
   end

src/lib/hsim/statefulRK45.ml

@@ -60,14 +60,11 @@ module InPlace =
       let step ({ state; vec=v } as s) h =
         let y_nv = vec v in
         let h = step state h y_nv in
-        let state = copy state in
         let get_dky t = get_dky state y_nv t 0; unvec y_nv in
         (h, get_dky), s in
 
       let copy { state; vec } =
-        let vec' = Zls.cmake 0 in
+        { state = copy state; vec = Zls.copy vec } in
-        Zls.blit vec vec';
-        { state; vec = vec' } in
 
       DNodeC { state; reset; step; copy }
   end

src/lib/hsim/types.ml

@@ -64,3 +64,8 @@ type ('p, 'a, 'b) lazy_sim =
     subsystem steps as needed to reach the input's horizon. *)
 type ('p, 'a, 'b) greedy_sim =
   ('p, 'a value, 'b value list) dnode
+
+(** Utils *)
+
+let d_of_dc (DNodeC { state; step; reset; _ }) =
+  DNode { state; step; reset }