feat: support for model-requested horizons

exm/ball.ml

@@ -29,7 +29,7 @@ let fder y yd =
   else begin yd.{0} <- 0.0; yd.{1} <- 0.0; yd.{2} <- 0.0; yd.{3} <- 0.0 end;
   yd
 let fzer y zo  = zo.{0} <- -. y.{1}; zo
-let fout _ _ y = of_array [| y.{1}; y.{0} |]
+let fout _ _ _ y = of_array [| y.{1}; y.{0} |]
 let jump _     = true
 let horizon _  = max_float
 let cget s     = s.lx
@@ -44,9 +44,9 @@ let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec =
   let yd         = cmake csize in
   let zout       = cmake zsize in
   let zfalse     = zmake 1 in
-  let fder _ _ y = fder y yd in
+  let fder _ _ _ y = fder y yd in
-  let fzer _ _ y = fzer y zout in
+  let fzer _ _ _  y = fzer y zout in
-  let step s _   = step s zfalse in
+  let step s _ _   = step s zfalse in
   let state      = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
   let reset _ _  = state in
   { state; fder; fzer; fout; step; reset; horizon; jump; cset; cget; zset;

exm/dune

@@ -5,6 +5,9 @@
 
 (library
  (name examples)
- (libraries hsim solvers))
+ (libraries hsim solvers std))
+
+; (executable
+;  (name ballz_main))
 
 (include_subdirs unqualified)

exm/sin1x.ml

@@ -58,10 +58,10 @@ let sin_1_x () =
   let zfalse = zmake 1 in
   let yd     = cmake 1 in
   let zout   = cmake 1 in
-  let fder _ _ y = fder 0.0 y yd; yd in
+  let fder _ _ _ y = fder 0.0 y yd; yd in
-  let fzer _ _ y = fzer 0.0 y zout; zout in
+  let fzer _ _ _ y = fzer 0.0 y zout; zout in
-  let fout s _ y = fout s y in
+  let fout s _ _ y = fout s y in
-  let step s _ = step s zfalse in
+  let step s _ _ = step s zfalse in
   let state = { s_x = of_array [| 0.0 |]; zin = zfalse } in
   { state; fder; fzer; fout; step; reset; horizon;
     cset; cget; zset; csize; zsize; jump }

exm/sin1x_der.ml

@@ -58,10 +58,10 @@ let sin_1_x z d =
   let zfalse = zmake 1 in
   let yd     = cmake 2 in
   let zout   = cmake 1 in
-  let fder _ _ y = fder d y yd; yd in
+  let fder _ _ _ y = fder d y yd; yd in
-  let fzer _ _ y = fzer z y zout; zout in
+  let fzer _ _ _ y = fzer z y zout; zout in
-  let fout s _ y = fout s y in
+  let fout s _ _ y = fout s y in
-  let step s _ = step s zfalse in
+  let step s _  _ = step s zfalse in
   let state = { s_x = of_array [| 0.0; y0 |]; zin = zfalse } in
   { state; fder; fzer; fout; step; horizon;
     cset; cget; zset; csize; zsize; reset; jump }

exm/sincos.ml

@@ -11,7 +11,7 @@ let zsize = 1
 
 let fder y yd omega =
   yd.{0} <- omega *. y.{1}; yd.{1} <- -.omega *. y.{0}; yd.{2} <- 1.0; yd
-let fout _ _ y = of_array [| y.{0}; y.{1}; y.{2} |]
+let fout _ _ _ y = of_array [| y.{0}; y.{1}; y.{2} |]
 let step { si; sx } sin0 cos0 =
   let sx = if si then of_array [| sin0; cos0; 0.0 |] else sx in
   of_array [| sx.{0}; sx.{1}; sx.{2} |], { sx; si=false }
@@ -26,9 +26,9 @@ let sinus_cosinus theta0 omega =
   let cos0 = Float.cos theta0 in
   let yd   = cmake csize in
   let zout = cmake zsize in
-  let fder _ _ y = fder y yd omega in
+  let fder _ _ _ y = fder y yd omega in
-  let fzer _ _ _ = zout in
+  let fzer _ _ _ _ = zout in
-  let step s _   = step s sin0 cos0 in
+  let step s _ _   = step s sin0 cos0 in
   let state      = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in
   let reset _ _  = state in
   HNode { state; fder; fzer; fout; step; reset; horizon;

exm/sqrt.ml

@@ -66,10 +66,10 @@ let sqrt () =
   let reset _ _ = s_init in
   let jump _ = true in
   HNode { state = s_init;
-    fder = (fun s _ y -> fder s y yd; yd);
+    fder = (fun s _ _ y -> fder s y yd; yd);
-    fzer = (fun s _ y -> fzero s y zout; zout);
+    fzer = (fun s _ _ y -> fzero s y zout; zout);
-    fout = (fun s _ y -> fout s y);
+    fout = (fun s _ _ y -> fout s y);
-    step = (fun s a -> fstep s a);
+    step = (fun s _ a -> fstep s a);
     horizon = (fun s -> if s.s_encore then 0.0 else max_float);
     cset; cget; zset; zsize; csize; reset; jump }
 

exm/vdp.ml

@@ -18,8 +18,8 @@ let fder y yd =
   yd.{0} <- y.{1};
   yd.{1} <- (mu *. (1.0 -. (y.{0} *. y.{0})) *. y.{1}) -. y.{0};
   yd
-let fout _ _ y = of_array [| y.{0}; y.{1} |]
+let fout _ _ _ y = of_array [| y.{0}; y.{1} |]
-let step { i; lx } _ =
+let step { i; lx } _ _ =
   let lx = if i then of_array [| x0; y0 |] else lx in
   of_array [| lx.{0}; lx.{1} |], { lx; i=false }
 let cget s    = s.lx
@@ -31,8 +31,8 @@ let horizon _ = max_float
 let van_der_pol () : (_, _, carray, carray, carray, zarray, carray) hnode =
   let yd         = cmake csize in
   let zout       = cmake zsize in
-  let fder _ _ y = fder y yd in
+  let fder _ _ _ y = fder y yd in
-  let fzer _ _ _ = zout in
+  let fzer _ _ _ _ = zout in
   let state      = { lx=of_array [| x0; y0 |]; i=true } in
   let reset _ _  = state in
   HNode { state; fder; fzer; fout; step; reset; horizon;

exm/zelus/ballz/ballz.zls

@@ -9,8 +9,7 @@ let hybrid ball (y0, y'0) = (y, y', z) where
 
 let hybrid main () =
   let der t = 1.0 init 0.0 in
-  let rec der p = 1.0 init -0.01 reset s -> -0.01
+  let s = period(0.01) in
-      and s = up(p) in 
   let (y, y', z) = ball (y0, y'0) in
   present z | s -> (
     print_float t;

exm/zelus/ballz/dune

@@ -1,6 +1,6 @@
 (rule
- (targets ballz.ml ballz.zci)
+ (targets ballz.ml ballz.zci ballz_main.ml)
  (deps
   (:zl ballz.zls))
  (action
-  (run zeluc %{zl})))
+  (run zeluc -s main -o ballz_main %{zl})))

exm/zelus/count/count.zls

@@ -0,0 +1,4 @@
+
+let node count () =
+  let rec n = 0 -> (pre n + 1) in
+  print_int n; print_newline ()

exm/zelus/count/dune

@@ -0,0 +1,6 @@
+(rule
+ (targets count.ml count.zci)
+ (deps
+  (:zl count.zls))
+ (action
+  (run zeluc %{zl})))

exm/zelus/sincos/sincosz.zls

@@ -1,4 +1,17 @@
 
-let hybrid f () = (sin, cos) where
+let hybrid g () = (sin, cos) where
   rec der sin = cos init 0.0
   and der cos = -. sin init 1.0
+
+let hybrid f () =
+  let der t = 1.0 init 0.0 in
+  let sin, cos = g () in
+  let z = period (0.01) in
+  present z -> (
+    print_float t;
+    print_string "\t";
+    print_float sin;
+    print_string "\t";
+    print_float cos;
+    print_newline ()
+  ); ()

exm/ztypes.ml

@@ -1,4 +1,21 @@
-include Common
+include Std
 include Ztypes
 include Solvers
 
+module type IGNORE = sig end
+module Defaultsolver : IGNORE = struct end
+
+module Zlsrun = struct
+  module Make (S : IGNORE) = struct
+    let go s =
+      let s = Lift.lift_hsim s in
+      let open Hsim in
+      let state = (module State.InPlaceSimState : State.SimState) in
+      let solver =
+        Solver.solver (StatefulSundials.InPlace.csolve)
+        (Types.d_of_dc StatefulZ.InPlace.zsolve) in
+      let open Sim.Sim(val state) in
+      ()
+      (* run_until_n (Utils.ignore 0 (run s solver)) 30. 1 ignore *)
+  end
+end

src/bin/main.ml

@@ -4,7 +4,7 @@ open Solvers
 open Examples
 open Common
 open Types
-open Lift
+open Std.Lift
 
 let sample   = ref 1
 let stop     = ref 10.0
@@ -53,12 +53,11 @@ let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
 let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
 
 let args = List.rev !modelargs
-let () = ignore lift
 
 let wrap_zelus (HNode m) =
   let ret = Bigarray.(Array1.create Float64 c_layout 0) in
-  let fout s a y = ignore (m.fout s a y); ret in
+  let fout s t a y = ignore (m.fout s t a y); ret in
-  let step s () = let _, s = m.step s () in ret, s in
+  let step s t () = let _, s = m.step s t () in ret, s in
   HNode { m with fout; step }
 
 let m =
@@ -67,7 +66,9 @@ let m =
       match !model with
       | None           -> Format.eprintf "Missing model\n"; exit 2
       | Some "ballz"   -> wrap_zelus (lift Ballz.main)
+      | Some "ballzm"  -> wrap_zelus (lift_hsim Ballz_main.main)
       | Some "sincosz" -> wrap_zelus (lift Sincosz.f)
+      (* | Some "count"   -> wrap_zelus (lift Count.count) *)
       | Some s         -> Format.eprintf "Unknown model: %s\n" s; exit 2
     else
       match !model with
@@ -85,7 +86,7 @@ let st = if !inplace then (module State.InPlaceSimState    : State.SimState)
                      else (module State.FunctionalSimState : State.SimState)
 
 let output =
-  if !no_print then Output.ignore
+  if !no_print then Hsim.Utils.ignore
   else if !speed then Output.print_h
   else Output.print (* Output.ignore *)
 

src/bin/output.ml

@@ -51,14 +51,3 @@ let print samples n =
 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, ((), ())) }
-
-let ignore _ n =
-  let state = () 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
-  DNode { n with reset=fun p -> n.reset (p, ()) }
-

src/lib/common/utils.ml

@@ -1,9 +0,0 @@
-
-exception TODO
-
-let pair = fun a b -> a, b
-
-let uncurry = fun f (a, b) -> f a b
-
-let (@.) = fun f g x -> f @@ g x
-

src/lib/hsim/sim.ml

@@ -11,34 +11,34 @@ module Sim (S : SimState) =
       let ms, ss, zin = get_mstate s, get_sstate s, get_zin s in
       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 (i.u now) in
+      let o, ms = step ms now (i.u now) in
       let s =
-        let s = set_zin None s in
         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 fder t = fder ms t (Utils.offset i.u now t) in
-          let fzer t = fzer ms (Utils.offset i.u now t) in
+          let fzer t = fzer ms t (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 input = { 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)
+          update ms ss (set_running ~mode ~input ~stop ~now s)
         end else set_running ~mode:Continuous s in
-      Utils.dot o, s
+      Utils.dot o, (set_zin None s)
 
     let step_continuous s step cset fout hor =
       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 stop = min stop (hor ms) in
       let (h, f, z), ss = step ss (min stop (hor ms)) in
+      let h = min h (min stop (hor ms)) in
       let ms = cset ms (f h) in
       let fy t = f (now +. t) in
       let fms t = cset ms (fy t) in
-      let fout t = fout ms (i.u (now +. t)) (fy t) in
+      let fout t = fout ms t (i.u (now +. t)) (fy t) in
       let s, c = match z with
       | None ->
           if h >= stop

src/lib/hsim/solver.ml

@@ -60,8 +60,10 @@ let solver (DNode csolver : ('y, 'yder) csolver)
     : ('y, 'yder, 'zin, 'zout) solver =
   let state = csolver.state, zsolver.state in
   let step (cstate, zstate) h =
-    let (h, f), cstate = csolver.step cstate h in
+    let (h', f), cstate = csolver.step cstate h in
-    let (h, z), zstate = zsolver.step zstate (h, f) in
+    let h = min h h' in
+    let (h', z), zstate = zsolver.step zstate (h, f) in
+    let h = min h h' in
     (h, f, z), (cstate, zstate) in
   let reset (ivp, zc) (cstate, zstate) =
     csolver.reset ivp cstate, zsolver.reset zc zstate in
@@ -73,9 +75,10 @@ let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
     : ('y, 'yder, 'zin, 'zout) solver_c =
   let state = csolver.state, zsolver.state in
   let step (cstate, zstate) h =
-    let (h, f), cstate = csolver.step cstate h in
+    let (h', f), cstate = csolver.step cstate h in
-    let (h, z), zstate = zsolver.step zstate (h, f) in
+    let h = min h h' in
-    (h, f, z), (cstate, zstate) in
+    let (h', z), zstate = zsolver.step zstate (h, f) in
+    (h', f, z), (cstate, zstate) in
   let reset (ivp, zc) (cstate, zstate) =
     csolver.reset ivp cstate, zsolver.reset zc zstate in
   let copy (cstate, zstate) =

src/lib/hsim/types.ml

@@ -35,10 +35,10 @@ type ('p, 'a, 'b) dnode_c =
 
 type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hrec =
   { state: 's;
-    step : 's -> 'a -> 'b * 's;        (** Discrete step function.            *)
+    step : 's -> time -> 'a -> 'b * 's;     (** Step function.                *)
-    fder : 's -> 'a -> 'y -> 'yder;    (** Continuous derivative function.    *)
+    fder : 's -> time -> 'a -> 'y -> 'yder; (** Derivative function.          *)
-    fout : 's -> 'a -> 'y -> 'b;       (** Continuous output function.        *)
+    fout : 's -> time -> 'a -> 'y -> 'b;    (** Output function.              *)
-    fzer : 's -> 'a -> 'y -> 'zout;    (** Continuous zero-crossing function. *)
+    fzer : 's -> time -> 'a -> 'y -> 'zout; (** Zero-crossing function.       *)
     reset : 'p -> 's -> 's;                 (** Reset function.               *)
     horizon : 's -> time;                   (** Next integration horizon.     *)
     jump : 's -> bool;                      (** Discontinuity flag.           *)

src/lib/hsim/utils.ml

@@ -78,3 +78,14 @@ let map f =
   let step () = function None -> None, () | Some v -> Some (f v), () in
   let reset () () = () in
   DNode { state; step; reset }
+
+let ignore _ n =
+  let state = () 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
+  DNode { n with reset=fun p -> n.reset (p, ()) }
+

src/lib/solvers/odexx.ml

@@ -52,8 +52,8 @@ struct (* {{{1 *)
   open Bigarray
 
   let debug () =
-    false
+    (* false *)
-    (* !Common.Debug.debug *)
+    !Common.Debug.debug
 
   let pow = 1.0 /. float(Butcher.order)
 

src/lib/std/dune

@@ -0,0 +1,3 @@
+(library
+ (name std)
+ (libraries hsim solvers))

src/lib/std/lift.ml

@@ -1,7 +1,7 @@
 
 open Hsim.Types
 open Solvers.Zls
-open Common.Ztypes
+open Ztypes
 
 type ('s, 'a) state =
   { mutable state : 's; mutable input : 'a option; mutable time : time }
@@ -25,7 +25,7 @@ let lift f =
   let no_time = -1.0 in
 
   (* the function that compute the derivatives *)
-  let fder { state; _ } input y =
+  let fder { state; time; _ } offset input y =
     cstate.major <- false;
     cstate.zinvec <- no_roots_in;
     cstate.zoutvec <- no_roots_out;
@@ -33,11 +33,11 @@ let lift f =
     cstate.dvec <- ignore_der;
     cstate.cindex <- 0;
     cstate.zindex <- 0;
-    ignore (f_step state (no_time, input));
+    ignore (f_step state (time +. offset, input));
     cstate.dvec in
 
   (* the function that compute the zero-crossings *)
-  let fzer { state; _ } input y =
+  let fzer { state; time; _ } offset input y =
     cstate.major <- false;
     cstate.zinvec <- no_roots_in;
     cstate.dvec <- ignore_der;
@@ -45,11 +45,11 @@ let lift f =
     cstate.cvec <- y;
     cstate.cindex <- 0;
     cstate.zindex <- 0;
-    ignore (f_step state (no_time, input));
+    ignore (f_step state (time +. offset, input));
     cstate.zoutvec in
 
   (* the function which compute the output during integration *)
-  let fout { state; _ } input y =
+  let fout { state; time; _ } offset input y =
     cstate.major <- false;
     cstate.zoutvec <- no_roots_out;
     cstate.dvec <- ignore_der;
@@ -57,10 +57,12 @@ let lift f =
     cstate.cvec <- y;
     cstate.cindex <- 0;
     cstate.zindex <- 0;
-    f_step state (no_time, input) in
+    f_step state (time +. offset, input) in
 
   (* the function which compute a discrete step *)
-  let step ({ state; time; _ } as st) input =
+  let step ({ state; time; _ } as st) offset input =
+    st.input <- Some input;
+    st.time <- time +. offset;
     cstate.major <- true;
     cstate.horizon <- infinity;
     cstate.zinvec <- no_roots_in;
@@ -68,13 +70,15 @@ let lift f =
     cstate.dvec <- ignore_der;
     cstate.cindex <- 0;
     cstate.zindex <- 0;
-    let o = f_step state (time, input) in
+    let o = f_step state (st.time, input) in
-    o, { st with state; input=Some input } in
+    o, st in
 
   let reset _ ({ state; _ } as st) = f_reset state; st in
 
   (* horizon *)
-  let horizon _ = cstate.horizon in
+  let horizon { time; _ } =
+    (* Printf.printf "\tCalling horizon :: cstate.horizon=%.10e\ttime=%.10e\n" cstate.horizon time; *)
+    cstate.horizon -. time in
 
   let jump _ = true in
 
@@ -116,3 +120,36 @@ let lift f =
     { state; fder; fzer; step; fout; reset;
       horizon; cset; cget; zset; zsize; csize; jump }
 
+let lift_hsim n =
+  let Hsim {
+    alloc; step; reset; derivative; crossings; maxsize; horizon; _
+  } = n in
+  let s = alloc () in
+  let state = { state = s; input = None; time = 0.0 } in
+  let csize, zsize = maxsize s in
+  let no_roots_in = zmake zsize in
+  let no_roots_out = cmake zsize in
+  let ignore_der = cmake csize in
+  let cstates = cmake csize in
+  let no_time = -1.0 in
+  reset s;
+  let fder { state; time; _ } offset () y =
+    derivative state y ignore_der no_roots_in no_roots_out (time +. offset);
+    ignore_der in
+  let fzer { state; time; _ } offset () y =
+    crossings state y no_roots_in no_roots_out (time +. offset); no_roots_out in
+  let fout _ _ () _ = () in
+  let step { state; time; _ } offset () =
+    step state cstates ignore_der no_roots_in (time +. offset),
+    { state; time=time +. offset; input=Some () } in
+  let reset _ ({ state; _ } as st) = reset state; st in
+  let horizon { state; time; _ } = horizon state -. time in
+  let jump _ = true in
+  let cset ({ state; _ } as st) _ =
+    derivative state cstates ignore_der no_roots_in no_roots_out no_time; st in
+  let zset ({ state; _ } as st) zinvec =
+    derivative state cstates ignore_der zinvec no_roots_out no_time; st in
+  let cget { state; _ } =
+    derivative state cstates ignore_der no_roots_in no_roots_out no_time; cstates in
+
+  HNode { state; fder; fzer; fout; step; reset; horizon; jump; cget; cset; zset; csize; zsize }