feat: a lot of stuff

src/lib/common/errors.ml

@@ -0,0 +1,2 @@
+exception TODO
+exception Internal of string

src/lib/hsim/dune

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

src/lib/hsim/sim.ml

@@ -5,13 +5,14 @@ open State
 
 module LazySim (S : SimState) =
   struct
+    module S = S
 
     (** "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)
-      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
-    = let state = S.get_init model.state solver.state in
+      (HNode model : ('ms, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNode solver : ('ss, 'y, 'yder, 'zin, 'zout) solver)
+      : (('a, 'ms, 'ss) S.state, 'p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) lazy_sim
+    = let init (p, s) = S.get_init (model.init p) (solver.init s) in
 
       let step s i =
         let ms, ss = S.get_mstate s, S.get_sstate s in
@@ -30,10 +31,10 @@ module LazySim (S : SimState) =
                   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 init   = model.cget ms and stop = stop -. now 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 ivp    = { fder; stop; init; size = model.csize } 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;
@@ -64,32 +65,29 @@ module LazySim (S : SimState) =
         let ss = solver.reset ps (S.get_sstate s) in
         S.update ms ss (S.set_idle s) in
 
-      DNode { state; step; reset }
+      DNode { init; 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 run_on model solver input p use =
       let DNode sim = run model solver in
-      let state = match sim.step sim.state (Some input) 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   -> ()
-        | Some o -> use o; loop (DNode { s with state }) in
-      loop (DNode { sim with state })
+      let state = sim.step (sim.init p) (Some input) in
+      let state = match state with None, s -> s | _ -> assert false in
+      let rec loop state =
+        let o, state = sim.step state None in
+        match o with None -> () | Some o -> use o; loop state in
+      loop state
 
     (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      ignore @@ List.fold_left (fun (DNode sim) i ->
-        let state = match sim.step sim.state (Some i) with
+    let run_on_n model solver inputs p use =
+      let DNode sim = run model solver in
+      ignore @@ List.fold_left (fun state i ->
+        let state = match sim.step 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
+        let rec loop state =
+          let o, state = sim.step state None in
+          match o with None -> state | Some o -> use o; loop state in
+        loop state) (sim.init p) inputs
 
     (** Run the model autonomously until [length], or until the model stops 
         answering. *)
@@ -108,12 +106,14 @@ module LazySim (S : SimState) =
 
 module GreedySim (S : SimState) =
   struct
+    module S = S
+
     (** "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)
-      : ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) greedy_sim
-    = let state = S.get_init model.state solver.state in
+      (HNode model : ('ms, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode)
+      (DNodeC solver : ('ss, 'y, 'yder, 'zin, 'zout) solver_c)
+      : (('a, 'ms, 'ss) S.state, 'p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) greedy_sim
+    = let init (m, s) = S.get_init (model.init m) (solver.init s) in
 
       let rec step s i =
         let ms, ss = S.get_mstate s, S.get_sstate s in
@@ -132,7 +132,7 @@ module GreedySim (S : SimState) =
                 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 ivp    = { fder; stop = stop -. now; init; size = model.csize } 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; 
@@ -171,20 +171,21 @@ module GreedySim (S : SimState) =
         let ss = solver.reset sp (S.get_sstate s) in
         S.update ms ss (S.set_idle s) in
 
-      DNode { state; step; reset }
+      DNode { init; 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 run_on model solver input p use =
       let DNode sim = run model solver in
-      let o, _ = sim.step sim.state input in
+      let o, _ = sim.step (sim.init p) input in
       List.iter use o
 
     (** Run the model on multiple inputs. *)
-    let run_on_n model solver inputs use =
-      let o, _ = List.fold_left (fun (acc, DNode sim) i ->
-        let o, state = sim.step sim.state i in
-        o::acc, DNode { sim with state }) ([], run model solver) inputs in
+    let run_on_n model solver inputs p use =
+      let DNode sim = run model solver in
+      let o, _ = List.fold_left (fun (acc, state) i ->
+        let o, state = sim.step state i in
+        o::acc, state) ([], sim.init p) inputs in
       List.iter use (List.concat (List.rev o))
 
     (** Run the model autonomously until [length], or until the model stops 

src/lib/hsim/solver.ml

@@ -5,7 +5,8 @@ open Types
 type ('y, 'yder) ivp =
   { init : 'y;                         (** [y₀]: initial value of y.          *)
     fder : time -> 'y -> 'yder;        (** [dy/dt]: derivative of y.          *)
-    stop : time }                      (** Stop time.                         *)
+    stop : time;                       (** Stop time.                         *)
+    size : int }
 
 (** A zero-crossing expression. *)
 type ('y, 'zout) zc =
@@ -17,72 +18,62 @@ type ('y, 'zout) zc =
   - an initial value problem as parameter;
   - an horizon to reach as input;
   - an actual time reached and dense solution as output *)
-type ('y, 'yder) csolver =
-  (('y, 'yder) ivp, time, time * (time -> 'y)) dnode
+type ('s, 'y, 'yder) csolver =
+  ('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode
 
 (** An ODE solver can optionally provide a state copy method, in which case
     greedy simulation is possible. *)
-type ('y, 'yder) csolver_c =
-  (('y, 'yder) ivp, time, time * (time -> 'y)) dnode_c
+type ('s, 'y, 'yder) csolver_c =
+  ('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode_c
 
 (** A zero-crossing solver is a synchronous function with:
   - a zero-crossing expression as parameter;
   - a time and dense solution as input;
   - an actual time reached and optional zero-crossing as output *)
-type ('y, 'zin, 'zout) zsolver =
-  (('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
+type ('s, 'y, 'zin, 'zout) zsolver =
+  ('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
 
 (** 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
+type ('s, 'y, 'zin, 'zout) zsolver_c =
+  ('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode_c
 
 (** A solver is a synchronous function with:
   - an initial value problem and zero-crossing expression as parameter;
   - an horizon to reach as input;
   - an actual time, dense solution and optional zero-crossing as output *)
-type ('y, 'yder, 'zin, 'zout) solver =
-  (('y, 'yder) ivp * ('y, 'zout) zc,
+type ('s, 'y, 'yder, 'zin, 'zout) solver =
+  ('s,
+   ('y, 'yder) ivp * ('y, 'zout) zc,
    time,
    time * (time -> 'y) * 'zin option) dnode
 
 (** 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,
+type ('s, 'y, 'yder, 'zin, 'zout) solver_c =
+  ('s,
+   ('y, 'yder) ivp * ('y, 'zout) zc,
    time,
    time * (time -> 'y) * 'zin option) dnode_c
 
-let csolver_from_c (DNodeC csolver : ('y, 'yder) csolver_c)
-  : ('y, 'yder) csolver
-= DNode { state = csolver.state; step = csolver.step; reset = csolver.reset }
-
-let zsolver_from_c (DNodeC zsolver : ('y, 'zin, 'zout) zsolver_c)
-  : ('y, 'zin, 'zout) zsolver
-= DNode { state = zsolver.state; step = zsolver.step; reset = zsolver.reset }
-
-let solver_from_c (DNodeC solver : ('y, 'yder, 'zin, 'zout) solver_c)
-  : ('y, 'yder, 'zin, 'zout) solver
-= DNode { state = solver.state; step = solver.step; reset = solver.reset }
-
 (** Build a full solver from an ODE solver and a zero-crossing solver. *)
-let solver (DNode csolver : ('y, 'yder) csolver)
-           (DNode zsolver : ('y, 'zin, 'zout) zsolver)
-    : ('y, 'yder, 'zin, 'zout) solver =
-  let state = csolver.state, zsolver.state in
+let solver (DNode csolver : ('sc, 'y, 'yder) csolver)
+           (DNode zsolver : ('sz, 'y, 'zin, 'zout) zsolver)
+    : ('sc * 'sz, 'y, 'yder, 'zin, 'zout) solver =
+  let init (ivp, zc) = csolver.init ivp, zsolver.init zc in
   let step (cstate, zstate) h =
     let (h, f), cstate = csolver.step cstate h 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
-  DNode { state; step; reset }
+  DNode { init ; step; reset }
 
 (** Build a full solver supporting state copies. *)
-let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
-             (DNodeC zsolver : ('y, 'zin, 'zout) zsolver_c)
-    : ('y, 'yder, 'zin, 'zout) solver_c =
-  let state = csolver.state, zsolver.state in
+let solver_c (DNodeC csolver : ('sc, 'y, 'yder) csolver_c)
+             (DNodeC zsolver : ('sz, 'y, 'zin, 'zout) zsolver_c)
+    : ('sc * 'sz, 'y, 'yder, 'zin, 'zout) solver_c =
+  let init (ivp, zc) = csolver.init ivp, zsolver.init zc in
   let step (cstate, zstate) h =
     let (h, f), cstate = csolver.step cstate h in
     let (h, z), zstate = zsolver.step zstate (h, f) in
@@ -91,5 +82,4 @@ let solver_c (DNodeC csolver : ('y, 'yder) csolver_c)
     csolver.reset ivp cstate, zsolver.reset zc zstate in
   let copy (cstate, zstate) =
     csolver.copy cstate, zsolver.copy zstate in
-  DNodeC { state; step; reset; copy }
-
+  DNodeC { init; step; reset; copy }

src/lib/hsim/statefulZ.ml

@@ -1,34 +0,0 @@
-
-open Types
-open Solvers
-open Solver
-
-module Functional =
-  struct
-
-    type ('state, 'vec) state = { state: 'state; vec: 'vec }
-
-    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
-      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;
-          vec = if Zls.length vec = size then vec else Zls.zmake size } in
-
-      let step ({ state; vec } as s) (h, fder) =
-        let y1 = fder h in
-        let fder h _ = let y = fder h in Zls.blit y y1 in
-        Illinois.step state h y1;
-        let v = Illinois.has_roots state in
-        if v then
-          let h = Illinois.find state (fder, y1) vec in
-          (h, Some vec), s
-        else (h, None), s in
-
-      let copy s = s in
-
-      DNodeC { state; step; reset; copy }
-
-  end

src/lib/hsim/types.ml

@@ -14,21 +14,21 @@ type 'a value =
 type 'a signal = 'a value option
 
 (** A discrete node. *)
-type ('p, 'a, 'b) dnode =
+type ('s, 'p, 'a, 'b) dnode =
   DNode :
-    { state : 'ds;
-      step  : 'ds -> 'a -> 'b * 'ds;
-      reset : 'p -> 'ds -> 'ds;
-    } -> ('p, 'a, 'b) dnode
+    { init  : 'p -> 's;
+      step  : 's -> 'a -> 'b * 's;
+      reset : 'p -> 's -> 's;
+    } -> ('s, 'p, 'a, 'b) dnode
 
 (** A discrete node which supports a state copy. *)
-type ('p, 'a, 'b) dnode_c =
+type ('s, 'p, 'a, 'b) dnode_c =
   DNodeC :
-    { state : 'ds;
-      step  : 'ds -> 'a -> 'b * 'ds;
-      reset : 'p -> 'ds -> 'ds;
-      copy  : 'ds -> 'ds;
-    } -> ('p, 'a, 'b) dnode_c
+    { init  : 'p -> 's;
+      step  : 's -> 'a -> 'b * 's;
+      reset : 'p -> 's -> 's;
+      copy  : 's -> 's;
+    } -> ('s, 'p, 'a, 'b) dnode_c
 
 (** A continuous node. *)
 type ('a, 'b, 'y, 'yder) cnode =
@@ -39,33 +39,33 @@ type ('a, 'b, 'y, 'yder) cnode =
     } -> ('a, 'b, 'y, 'yder) cnode
 
 (** A hybrid node. *)
-type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
+type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
   HNode :
-    { state : 'hs;
-      step : 'hs -> 'a -> 'b * 'hs;    (** Discrete step function.            *)
-      fder : 'hs -> 'a -> 'y -> 'yder; (** Continuous derivative function.    *)
-      fout : 'hs -> 'a -> 'y -> 'b;    (** Continuous output function.        *)
-      fzer : 'hs -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *)
-      reset : 'p -> 'hs -> 'hs;        (** Reset function.                    *)
-      horizon : 'hs -> time;           (** Next integration horizon.          *)
-      jump : 'hs -> bool;              (** Discontinuity flag.                *)
-      cget : 'hs -> 'y;                (** Get continuous state.              *)
-      cset : 'hs -> 'y -> 'hs;         (** Set continuous state.              *)
-      zset : 'hs -> 'zin -> 'hs;       (** Set zero-crossing state.           *)
+    { init : 'p -> 's;
+      step : 's -> 'a -> 'b * 's;    (** Discrete step function.            *)
+      fder : 's -> 'a -> 'y -> 'yder; (** Continuous derivative function.    *)
+      fout : 's -> 'a -> 'y -> 'b;    (** Continuous output function.        *)
+      fzer : 's -> 'a -> 'y -> 'zout; (** Continuous zero-crossing function. *)
+      reset : 'p -> 's -> 's;        (** Reset function.                    *)
+      horizon : 's -> time;           (** Next integration horizon.          *)
+      jump : 's -> bool;              (** Discontinuity flag.                *)
+      cget : 's -> 'y;                (** Get continuous state.              *)
+      cset : 's -> 'y -> 's;         (** Set continuous state.              *)
+      zset : 's -> 'zin -> 's;       (** Set zero-crossing state.           *)
+      csize : int;
       zsize : int;
-    } -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
+    } -> ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
 
 (** The simulation of a hybrid system is a synchronous function on streams of
     functions. *)
-type ('p, 'a, 'b) lazy_sim =
-  ('p, 'a signal, 'b signal) dnode
+type ('s, 'p, 'a, 'b) lazy_sim =
+  ('s, 'p, 'a signal, 'b signal) dnode
 
 (** Greedy simulation takes in an input and computes as many solver and
     subsystem steps as needed to reach the input's horizon. *)
-type ('p, 'a, 'b) greedy_sim =
-  ('p, 'a value, 'b value list) dnode
+type ('s, 'p, 'a, 'b) greedy_sim =
+  ('s, 'p, 'a value, 'b value list) dnode
 
 (** Utils *)
 
-let d_of_dc (DNodeC { state; step; reset; _ }) =
-  DNode { state; step; reset }
+let d_of_dc (DNodeC { init; step; reset; _ }) = DNode { init; step; reset }

src/lib/solvers/csolver.ml

@@ -0,0 +1,28 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module type Csolver =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val csolve : ((session, vec) state, carray, carray) csolver
+  end
+
+module type CsolverC =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val csolve : ((session, vec) state, carray, carray) csolver_c
+  end
+
+module CsolverOfC =
+  functor (S : CsolverC) -> (struct
+    type ('a, 'b) state = ('a, 'b) S.state
+    type session = S.session
+    type vec = S.vec
+    let csolve = d_of_dc S.csolve
+  end : Csolver)

src/lib/solvers/dune

@@ -1,4 +1,8 @@
-(env (dev (flags (:standard -w -9-27-32))))
+(env (dev (flags (:standard -w -32))))
 
 (library
-  (name solvers))
+  (name solvers)
+  (libraries
+    hsim
+    ;sundialsml
+    ))

src/lib/solvers/illinois.ml

@@ -27,7 +27,7 @@ let get_check_root rdir =
   let check_down   x0 x1 = if x0 > 0.0 && x1 <= 0.0 then -1l else 0l in
   let check_either x0 x1 = if x0 < 0.0 && x1 >= 0.0 then  1l else
                            if x0 > 0.0 && x1 <= 0.0 then -1l else 0l in
-  let no_check x0 x1 = 0l in
+  let no_check _x0 _x1 = 0l in
 
   match rdir with
   | Up     -> check_up
@@ -118,7 +118,7 @@ type t = {
   }
 
 (* Called from find when bothf_valid = false to initialise f1. *)
-let reinitialize ({ g; f1 = f1; t1 = t1 } as s) t c =
+let reinitialize ({ g; f1 = f1; t1 = t1; _ } as s) t c =
   s.t1 <- t;
   g t1 c f1;   (* fill f1, because it is immediately copied into f0 by next_mesh *)
   if !debug then (printf "z|---------- init(%.24e, ... ----------@." t;
@@ -148,10 +148,10 @@ let initialize nroots g c =
   s
 
 
-let num_roots { f0 } = Zls.length f0
+let num_roots { f0; _ } = Zls.length f0
 
 (* f0/t0 take the previous values of f1/t1, f1/t1 are refreshed by g *)
-let step ({ g; f0 = f0; f1 = f1; t1 = t1 } as s) t c =
+let step ({ g; f0 = f0; f1 = f1; t1 = t1; _ } as s) t c =
   (* swap f0 and f1; f0 takes the previous value of f1 *)
   s.f0 <- f1;
   s.t0 <- t1;
@@ -184,7 +184,7 @@ let resolve_intervals r1 r2 =
                    (possible) zero-crossing in (f_mid, f_right]
  *)
 let check_interval calc_zc f_left f_mid =
-  let check i r x0 x1 =
+  let check _i r x0 x1 =
     let rv = calc_zc x0 x1 in
     let r' = if rv = 0l then SearchRight
              else if x1 = 0.0 then FoundMid
@@ -340,17 +340,17 @@ 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 has_roots { bothf_valid; f0; f1; calc_zc; _ } =
+  bothf_valid && (check_interval calc_zc f0 f1 <> SearchRight)
 
-let takeoff { bothf_valid = bothf_valid; f0; f1 } =
+let takeoff { bothf_valid; f0; f1; _ } =
   bothf_valid && (takeoff f0 f1)
 
 (* returns true if a signal has moved from zero to a stritly positive value *)
 (* Added by MP. Ask Tim if this code is necessary, that is, what happens *)
 (* 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 find_takeoff ({ f0; f1; _ } as s) roots =
   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

src/lib/solvers/odexx.ml

@@ -156,7 +156,7 @@ struct (* {{{1 *)
   (* NB: y must be the initial state vector (y_0)
    *     k(0) must be the initial deriviatives vector (dy_0) *)
   let initial_stepsize { initial_step_size; abs_tol; rel_tol; max_step;
-                         time; y; hmax; k } =
+                         time; y; hmax; k; _ } =
     let hmin = 16.0 *. epsilon_float *. abs_float time in
     match initial_step_size with
     | Some h -> minmax hmin max_step h
@@ -168,7 +168,8 @@ struct (* {{{1 *)
       in
       max hmin (if hmax *. rh > 1.0 then 1.0 /. rh else hmax)
 
-  let reinitialize ?rhsfn ({ stop_time; min_step; max_step; sysf } as s) t ny =
+  let reinitialize
+      ?rhsfn ({ stop_time; min_step; max_step; sysf; _ } as s) t ny =
     Bigarray.Array1.blit ny s.y;
     s.time      <- t;
     s.last_time <- t;
@@ -250,9 +251,9 @@ struct (* {{{1 *)
 
   (* TODO: add stats: nfevals, nfailed, nsteps *)
   let step s t_limit user_y =
-    let { stop_time; min_step; abs_tol; rel_tol;
+    let { stop_time; abs_tol; rel_tol;
           sysf = f; time = t; h = h; hmax = hmax;
-          k = k; y = y; yold = ynew; } = s in
+          k = k; y = y; yold = ynew; _ } = s in
 
     (* First Same As Last (FSAL) swap; doing it after the previous
        step invalidates the interpolation routine. *)
@@ -323,7 +324,7 @@ struct (* {{{1 *)
     s.h <- nexth;
     s.time
 
-  let get_dky { last_time = t; time = t'; h = h; yold = y; k = k } yi ti kd =
+  let get_dky { last_time = t; time = t'; yold = y; k; _ } yi ti kd =
 
     if kd > 0 then
       failwith
@@ -355,11 +356,11 @@ struct (* {{{1 *)
     done
 
   (* copy functions *)
-  let copy ({ last_time; time; h; yold; k } as s) =
+  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 }
-      ({ last_time; time; h; yold; k } as s2) =
+  let blit { last_time = l1; time = t1; yold = yhold1; k = k1; _ }
+      ({ yold; k; _ } as s2) =
     s2.last_time <- l1; s2.time <- t1;
     Zls.blit yhold1 yold; Zls.blit_matrix k1 k
 

src/lib/solvers/statefulRK45.ml

@@ -1,25 +1,24 @@
 
-open Types
-open Solvers
-open Solver
+open Hsim.Types
+open Hsim.Solver
+open Zls
 
-module Functional =
+module Functional : Csolver.CsolverC =
   struct
     type ('state, 'vec) state = { state: 'state; vec: 'vec }
+    type session = Odexx.Ode45.t
+    type vec = carray
 
-    let csolve : (Zls.carray, Zls.carray) csolver_c =
+    let csolve : ((session, vec) state, carray, carray) csolver_c =
       let open Odexx.Ode45 in
 
-      let state =
+      let init _ =
         let v = Zls.cmake 0 in
         let state = initialize (fun _ _ _ -> ()) (vec v) in
         set_stop_time state 1.0; { state; vec=v } in
 
-      let reset
-        ({ fder; init; stop }: (Zls.carray, Zls.carray) ivp)
-        (_: (t, Zls.carray) state)
-        : (t, Zls.carray) state
-      = let fder t cvec dvec = Zls.blit (fder t cvec) dvec in
+      let reset { fder; init; stop; _ } _ = 
+        let fder t cvec dvec = Zls.blit (fder t cvec) dvec in
         let state = initialize fder (vec init) in
         set_stop_time state stop;
         { state; vec = init } in
@@ -33,25 +32,25 @@ module Functional =
 
       let copy { state; vec } = { state; vec } in
 
-      DNodeC { state; step; reset; copy }
+      DNodeC { init; step; reset; copy }
   end
 
-module InPlace =
+module InPlace : Csolver.CsolverC =
   struct
+    type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
+    type session = Odexx.Ode45.t
+    type vec = carray
 
-    type ('state, 'vec) state =
-      { mutable state: 'state; mutable vec : 'vec }
-
-    let csolve : (Zls.carray, Zls.carray) csolver_c =
+    let csolve : ((session, vec) state, carray, carray) csolver_c =
       let open Odexx.Ode45 in
 
-      let state =
+      let init _ =
         let v = Zls.cmake 0 in
         let state = initialize (fun _ _ _ -> ()) (vec v) in
         set_stop_time state 1.0;
         { state; vec=v } in
 
-      let reset { fder: time -> Zls.carray -> Zls.carray; init; stop } s =
+      let reset { fder; init; stop; _ } s =
         let fder t cvec dvec =
           let dvec' = fder t cvec in Zls.blit dvec' dvec in
         let state = initialize fder (vec init) in
@@ -66,5 +65,5 @@ module InPlace =
       let copy { state; vec } =
         { state = copy state; vec = Zls.copy vec } in
 
-      DNodeC { state; reset; step; copy }
+      DNodeC { init; reset; step; copy }
   end

src/lib/solvers/statefulSundials.ml

@@ -0,0 +1,71 @@
+(*
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module Functional : Csolver.Csolver =
+  struct
+    type ('state, 'vec) state = { state : 'state; vec : 'vec }
+    type session = (Sundials_RealArray.t, Nvector_serial.kind) Cvode.session
+    type vec = carray
+
+    let csolve : ((session, vec) state, carray, carray) csolver =
+      let open Cvode in
+
+      let init { size; fder=_; _ } =
+        let vec = cmake size in
+        let state = init Adams default_tolerances (fun _ _ _ -> ()) 0.
+                         (Nvector_serial.wrap vec) in
+        set_stop_time state 1.0;
+        { state; vec } in
+
+      let reset { init=i; fder; stop; _ } { vec; _ } =
+        let fder t cvec dvec =
+          let dvec' = fder t cvec in blit dvec' dvec in
+        let state =
+          Cvode.init Adams default_tolerances fder 0. (Nvector_serial.wrap i) in
+        set_stop_time state stop;
+        { state; vec } in
+
+      let step ({ state; vec } as s) h =
+        let y = Nvector_serial.wrap vec in
+        let h, _ = solve_one_step state h y in
+        let f t = get_dky state y t 0; Nvector_serial.unwrap y in
+        (h, f), s in
+
+      DNode { init; reset; step }
+  end
+
+module InPlace : Csolver.Csolver =
+  struct
+    type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
+
+    type session = (Sundials_RealArray.t, Nvector_serial.kind) Cvode.session
+    type vec = carray
+
+    let csolve : ((session, vec) state, carray, carray) csolver =
+      let open Cvode in
+
+      let init { size; fder=_; _ } =
+        let vec = cmake size in
+        let state = init Adams default_tolerances (fun _ _ _ -> ()) 0.
+                         (Nvector_serial.wrap vec) in
+        set_stop_time state 1.0;
+        { state; vec } in
+
+      let reset { init=i; fder; _ } s =
+        let fder t cvec dvec =
+          let dvec' = fder t cvec in blit dvec' dvec in
+        let state =
+          Cvode.init Adams default_tolerances fder 0. (Nvector_serial.wrap i) in
+        set_stop_time state 1.0; s.state <- state; s.vec <- i; s in
+
+      let step s h =
+        let y = Nvector_serial.wrap s.vec in
+        let h, _ = solve_one_step s.state h y in
+        let f t = get_dky s.state y t 0; Nvector_serial.unwrap y in
+        (h, f), s in 
+
+      DNode { init; reset; step }
+  end
+*)

src/lib/solvers/statefulZ.ml

@@ -0,0 +1,69 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module Functional : Zsolver.ZsolverC =
+  struct
+
+    type ('state, 'vec) state = { state: 'state; vec: 'vec }
+    type session = Illinois.t
+    type vec = zarray
+
+    let zsolve : ((session, vec) state, carray, vec, carray) zsolver_c =
+      let open Illinois in
+
+      let init _ =
+        { state = initialize 0 (fun _ _ _ -> ()) (cmake 0);
+          vec = zmake 0 } in
+
+      let reset { fzer; init; size } { vec; _ } =
+        let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in
+        { state = initialize size fzer init;
+          vec = if length vec = size then vec else zmake size } in
+
+      let step ({ state; vec } as s) (h, fder) =
+        let y1 = fder h in
+        let fder h _ = let y = fder h in blit y y1 in
+        step state h y1;
+        if has_roots state then
+          let h = find state (fder, y1) vec in
+          (h, Some vec), s
+        else (h, None), s in
+
+      let copy s = s in
+
+      DNodeC { init; step; reset; copy }
+  end
+
+module InPlace : Zsolver.ZsolverC =
+  struct
+    type ('state, 'vec) state = { mutable state : 'state; mutable vec : 'vec }
+    type session = Illinois.t
+    type vec = zarray
+
+    let zsolve : ((session, vec) state, carray, vec, carray) zsolver_c =
+      let open Illinois in
+
+      let init _ =
+        { state=initialize 0 (fun _ _ _ -> ()) (cmake 0);
+          vec=zmake 0 } in
+
+      let reset { size; init; fzer } s =
+        let fzer t cvec zout = let zout' = fzer t cvec in blit zout' zout in
+        s.state <- initialize size fzer init;
+        if length s.vec <> size then s.vec <- zmake size; s in
+
+      let step ({ state; vec } as s) (h, fder) =
+        let y = fder h in
+        let fder h _ = let y' = fder h in blit y' y in
+        step state h y;
+        if has_roots state then
+          let h = find state (fder, y) vec in
+          (h, Some vec), s
+        else (h, None), s in
+
+      let copy _ = raise Common.Errors.TODO in
+      
+      DNodeC { init; step; reset; copy }
+  end

src/lib/solvers/zsolver.ml

@@ -0,0 +1,28 @@
+
+open Hsim.Types
+open Hsim.Solver
+open Zls
+
+module type Zsolver =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val zsolve : ((session, vec) state, carray, zarray, carray) zsolver
+  end
+
+module type ZsolverC =
+  sig
+    type ('a, 'b) state
+    type session
+    type vec
+    val zsolve : ((session, vec) state, carray, zarray, carray) zsolver_c
+  end
+
+module ZsolverOfC =
+  functor (S : ZsolverC) -> (struct
+    type ('a, 'b) state = ('a, 'b) S.state
+    type session = S.session
+    type vec = S.vec
+    let zsolve = d_of_dc S.zsolve
+  end : Zsolver)