hms/hsim
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

feat: re-hide inner states, properly reset sundials

Browse source
This commit is contained in:
Henri Saudubray 2025-05-12 17:27:52 +02:00
parent 80d4aef23f
commit 76dc461d44
Signed by: hms
GPG key ID: 7065F57ED8856128
14 changed files with 162 additions and 253 deletions
Download patch file Download diff file Expand all files Collapse all files

11
exm/ball.ml
View file

@ -40,17 +40,16 @@ let step ({ zin; lx; _ } as s) zfalse =
of_array [| -. 0.8 *. lx.{0}; 0.0; lx.{2}; lx.{3} |] else lx in
of_array [| s.lx.{1} |], { zin=zfalse; lx; i=false }
let bouncing_ball ()
: (state, _, _, carray, carray, carray, zarray, carray) hnode
= let yd = cmake csize in
let bouncing_ball () : (_, _, carray, carray, carray, zarray, carray) hnode =
let yd = cmake csize in
let zout = cmake zsize in
let zfalse = zmake 1 in
let fder _ _ y = fder y yd in
let fzer _ _ y = fzer y zout in
let step s _ = step s zfalse in
let init _ = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
let reset _ _ = init () in
HNode { init; fder; fzer; fout; step; reset; horizon;
let state = { zin=zfalse; lx=of_array [|y'0;y0;x'0;x0|]; i=true } in
let reset _ _ = state in
HNode { state; fder; fzer; fout; step; reset; horizon;
jump; cset; cget; zset; csize; zsize }
let errmsg = "Too many arguments for the model (needed: 0)"

3
exm/model.ml
View file

@ -3,6 +3,5 @@ open Solvers.Zls
module type Model =
sig
type state
val init : string list -> (state, 'b, 'c, carray, carray, carray, zarray, carray) hnode
val init : string list -> ('b, 'c, carray, carray, carray, zarray, carray) hnode
end

6
exm/sincos.ml
View file

@ -29,9 +29,9 @@ let sinus_cosinus theta0 omega =
let fder _ _ y = fder y yd omega in
let fzer _ _ _ = zout in
let step s _ = step s sin0 cos0 in
let init _ = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in
let reset _ _ = init () in
HNode { init; fder; fzer; fout; step; reset; horizon;
let state = { sx=of_array [| sin0; cos0; 0.0 |]; si=true } in
let reset _ _ = state in
HNode { state; fder; fzer; fout; step; reset; horizon;
jump; cset; cget; zset; csize; zsize }
let errmsg_invalid = "Invalid arguments to model (needed: 2 floats)"

2
exm/sqrt.ml
View file

@ -65,7 +65,7 @@ let sqrt () =
s_zin = zmake 1 } in
let reset _ _ = s_init in
let jump _ = true in
HNode { init = (fun _ -> s_init);
HNode { state = s_init;
fder = (fun s _ y -> fder s y yd; yd);
fzer = (fun s _ y -> fzero s y zout; zout);
fout = (fun s _ y -> fout s y);

11
exm/vdp.ml
View file

@ -28,15 +28,14 @@ let zset s _ = s
let jump _ = true
let horizon _ = max_float
let van_der_pol ()
: (state, _, _, carray, carray, carray, zarray, carray) hnode
= let yd = cmake csize in
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 fzer _ _ _ = zout in
let init _ = { lx=of_array [| x0; y0 |]; i=true } in
let reset _ _ = init () in
HNode { init; fder; fzer; fout; step; reset; horizon;
let state = { lx=of_array [| x0; y0 |]; i=true } in
let reset _ _ = state in
HNode { state; fder; fzer; fout; step; reset; horizon;
jump; cset; cget; zset; csize; zsize }
let errmsg = "Too many arguments for the model (needed: 0)"

39
src/bin/main.ml
View file

@ -37,48 +37,35 @@ let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
let () = try Arg.parse (Arg.align opts) set_model errmsg with _ -> exit 2
let m =
match !model with
try match !model with
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ball" -> (module Ball : Model.Model)
| Some "vdp" -> (module Vdp)
| Some "sincos" -> (module Sincos)
| Some "sqrt" -> (module Sqrt)
| Some "ball" -> Ball.init !modelargs
| Some "vdp" -> Vdp.init !modelargs
| Some "sincos" -> Sincos.init !modelargs
| Some "sqrt" -> Sqrt.init !modelargs
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
let z = if !inplace then (module StatefulZ.InPlace : Zsolver.ZsolverC)
else (module StatefulZ.Functional : Zsolver.ZsolverC)
let z = StatefulZ.(if !inplace then InPlace.zsolve else Functional.zsolve)
let st = if !inplace then (module State.InPlaceSimState : State.SimState)
else (module State.FunctionalSimState : State.SimState)
let () =
let open (val m) in
let m = init !modelargs in
let sim =
if !sundials then
if !greedy then
(Format.eprintf "Sundials does not support greedy simulation\n"; exit 2)
else
let open StatefulSundials in
let c = if !inplace then (module InPlace : Csolver.Csolver)
else (module Functional : Csolver.Csolver) in
let open (val c) in let open (val z) in
let s = Solver.solver csolve (d_of_dc zsolve) in
let c = if !inplace then InPlace.csolve else Functional.csolve in
let s = Solver.solver c (d_of_dc z) in
let open Sim.LazySim(val st) in run_until_n m s
else
let open StatefulRK45 in
let c = if !inplace then (module InPlace : Csolver.CsolverC)
else (module Functional : Csolver.CsolverC) in
let open (val c) in let open (val z) in
let s = Solver.solver_c csolve zsolve in
let c = if !inplace then InPlace.csolve else Functional.csolve in
let s = Solver.solver_c c z in
if !greedy then let open Sim.GreedySim(val st) in run_until_n m s
else let open Sim.LazySim(val st) in run_until_n m (d_of_dc s)
in
let open Solver in
let HNode { init; csize; zsize; fder; fzer; cget; _ } = m in
let state = init () in
let init = cget state in
let ivp = { size=csize; fder=(fun _ -> fder state ()); init; stop=1.0 } in
let zc = { size=zsize; fzer=(fun _ -> fzer state ()); init } in
sim !stop !steps ((), (ivp, zc)) (Output.print !sample)
let () = sim !stop !steps (Output.print !sample)

114
src/lib/hsim/sim.ml
View file

@ -5,32 +5,32 @@ open State
module LazySim (S : SimState) =
struct
module S = S
include S
(** "Lazy" simulation of a model with any solver. *)
let run
(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
(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 = get_init model.state solver.state in
let step s i =
let ms, ss = S.get_mstate s, S.get_sstate s in
match i, S.is_running s with
let ms, ss = get_mstate s, get_sstate s in
match i, 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, 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
let i, now, stop = get_input s, get_now s, get_stop s in
match 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
if h <= 0.0 then set_mstate ms s
else if now >= stop then set_idle s
else if model.jump ms then begin
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
@ -40,8 +40,8 @@ module LazySim (S : SimState) =
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
update ms ss (set_running ~mode ~input:i ~stop ~now s)
end else 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
@ -52,26 +52,26 @@ module LazySim (S : SimState) =
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
then set_running ~mode:Discrete ~now:h s
else set_running ~now:h s in
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
let s = set_running ~mode:Discrete ~now:h s in
update (model.zset ms z) ss s in
Some out, s in
let reset (pm, ps) s =
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
let ms = model.reset pm (get_mstate s) in
let ss = solver.reset ps (get_sstate s) in
update ms ss (set_idle s) in
DNode { init; step; reset }
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 p use =
let run_on model solver input use =
let DNode sim = run model solver in
let state = sim.step (sim.init p) (Some input) in
let state = sim.step sim.state (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
@ -79,7 +79,7 @@ module LazySim (S : SimState) =
loop state
(** Run the model on multiple inputs. *)
let run_on_n model solver inputs p use =
let run_on_n model solver inputs use =
let DNode sim = run model solver in
ignore @@ List.fold_left (fun state i ->
let state = match sim.step state (Some i) with
@ -87,7 +87,7 @@ module LazySim (S : SimState) =
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
loop state) sim.state inputs
(** Run the model autonomously until [length], or until the model stops
answering. *)
@ -106,28 +106,28 @@ module LazySim (S : SimState) =
module GreedySim (S : SimState) =
struct
module S = S
include S
(** "Greedy" simulation of a model with an appropriate solver. *)
let run
(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
(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 = get_init model.state solver.state in
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 ms, ss = get_mstate s, get_sstate s in
if not (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
step (set_running ~mode ~input:i ~now ~stop s) i
else let now, stop = get_now s, get_stop s in
match 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.set_idle s
if h <= 0.0 then step (set_mstate ms s) i
else if now >= stop then [], 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
@ -138,9 +138,9 @@ module GreedySim (S : SimState) =
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
let s = set_running ~mode ~input:i ~stop ~now s in
step (update ms ss s) i
else step (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
@ -151,41 +151,41 @@ module GreedySim (S : SimState) =
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
let s = set_running ~mode:Discrete ~now:h s in
let rest, s = step (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
let s = set_running ~now:h s in
let rest, s = step (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 = 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
let rest, s = step (update ms ss s) i in
out::rest, s in
let reset (mp, sp) s =
let ms = model.reset mp (S.get_mstate s) in
let ss = solver.reset sp (S.get_sstate s) in
S.update ms ss (S.set_idle s) in
let ms = model.reset mp (get_mstate s) in
let ss = solver.reset sp (get_sstate s) in
update ms ss (set_idle s) in
DNode { init; step; reset }
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 p use =
let run_on model solver input use =
let DNode sim = run model solver in
let o, _ = sim.step (sim.init p) input in
let o, _ = sim.step sim.state input in
List.iter use o
(** Run the model on multiple inputs. *)
let run_on_n model solver inputs p use =
let run_on_n model solver inputs 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
o::acc, state) ([], sim.state) inputs in
List.iter use (List.concat (List.rev o))
(** Run the model autonomously until [length], or until the model stops

46
src/lib/hsim/solver.ml
View file

@ -18,62 +18,60 @@ 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 ('s, 'y, 'yder) csolver =
('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode
type ('y, 'yder) csolver =
(('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 ('s, 'y, 'yder) csolver_c =
('s, ('y, 'yder) ivp, time, time * (time -> 'y)) dnode_c
type ('y, 'yder) csolver_c =
(('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 ('s, 'y, 'zin, 'zout) zsolver =
('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode
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
case greedy simulation is possible. *)
type ('s, 'y, 'zin, 'zout) zsolver_c =
('s, ('y, 'zout) zc, time * (time -> 'y), time * 'zin option) dnode_c
type ('y, 'zin, 'zout) zsolver_c =
(('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 ('s, 'y, 'yder, 'zin, 'zout) solver =
('s,
('y, 'yder) ivp * ('y, 'zout) zc,
type ('y, 'yder, 'zin, 'zout) solver =
(('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 ('s, 'y, 'yder, 'zin, 'zout) solver_c =
('s,
('y, 'yder) ivp * ('y, 'zout) zc,
type ('y, 'yder, 'zin, 'zout) solver_c =
(('y, 'yder) ivp * ('y, 'zout) zc,
time,
time * (time -> 'y) * 'zin option) dnode_c
(** Build a full solver from an ODE solver and a zero-crossing solver. *)
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 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 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 { init ; step; reset }
DNode { state; step; reset }
(** Build a full solver supporting state copies. *)
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 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 step (cstate, zstate) h =
let (h, f), cstate = csolver.step cstate h in
let (h, z), zstate = zsolver.step zstate (h, f) in
@ -82,4 +80,4 @@ let solver_c (DNodeC csolver : ('sc, '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 { init; step; reset; copy }
DNodeC { state; step; reset; copy }

26
src/lib/hsim/types.ml
View file

@ -14,21 +14,21 @@ type 'a value =
type 'a signal = 'a value option
(** A discrete node. *)
type ('s, 'p, 'a, 'b) dnode =
type ('p, 'a, 'b) dnode =
DNode :
{ init : 'p -> 's;
{ state : 's;
step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's;
} -> ('s, 'p, 'a, 'b) dnode
} -> ('p, 'a, 'b) dnode
(** A discrete node which supports a state copy. *)
type ('s, 'p, 'a, 'b) dnode_c =
type ('p, 'a, 'b) dnode_c =
DNodeC :
{ init : 'p -> 's;
{ state : 's;
step : 's -> 'a -> 'b * 's;
reset : 'p -> 's -> 's;
copy : 's -> 's;
} -> ('s, 'p, 'a, 'b) dnode_c
} -> ('p, 'a, 'b) dnode_c
(** A continuous node. *)
type ('a, 'b, 'y, 'yder) cnode =
@ -39,9 +39,9 @@ type ('a, 'b, 'y, 'yder) cnode =
} -> ('a, 'b, 'y, 'yder) cnode
(** A hybrid node. *)
type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
type ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
HNode :
{ init : 'p -> 's;
{ state: '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. *)
@ -54,18 +54,16 @@ type ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode =
zset : 's -> 'zin -> 's; (** Set zero-crossing state. *)
csize : int;
zsize : int;
} -> ('s, 'p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
} -> ('p, 'a, 'b, 'y, 'yder, 'zin, 'zout) hnode
(** The simulation of a hybrid system is a synchronous function on streams of
functions. *)
type ('s, 'p, 'a, 'b) lazy_sim =
('s, 'p, 'a signal, 'b signal) dnode
type ('p, 'a, 'b) lazy_sim = ('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 ('s, 'p, 'a, 'b) greedy_sim =
('s, 'p, 'a value, 'b value list) dnode
type ('p, 'a, 'b) greedy_sim = ('p, 'a value, 'b value list) dnode
(** Utils *)
let d_of_dc (DNodeC { init; step; reset; _ }) = DNode { init; step; reset }
let d_of_dc (DNodeC { state; step; reset; _ }) = DNode { state; step; reset }

28
src/lib/solvers/csolver.ml
View file

@ -1,28 +0,0 @@
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)

20
src/lib/solvers/statefulRK45.ml
View file

@ -3,16 +3,14 @@ open Hsim.Types
open Hsim.Solver
open Zls
module Functional : Csolver.CsolverC =
module Functional =
struct
type ('state, 'vec) state = { state: 'state; vec: 'vec }
type session = Odexx.Ode45.t
type vec = carray
let csolve : ((session, vec) state, carray, carray) csolver_c =
let csolve : (carray, carray) csolver_c =
let open Odexx.Ode45 in
let init _ =
let state =
let v = Zls.cmake 0 in
let state = initialize (fun _ _ _ -> ()) (vec v) in
set_stop_time state 1.0; { state; vec=v } in
@ -32,19 +30,17 @@ module Functional : Csolver.CsolverC =
let copy { state; vec } = { state; vec } in
DNodeC { init; step; reset; copy }
DNodeC { state; step; reset; copy }
end
module InPlace : Csolver.CsolverC =
module InPlace =
struct
type ('state, 'vec) state = { mutable state: 'state; mutable vec : 'vec }
type session = Odexx.Ode45.t
type vec = carray
let csolve : ((session, vec) state, carray, carray) csolver_c =
let csolve : (carray, carray) csolver_c =
let open Odexx.Ode45 in
let init _ =
let state =
let v = Zls.cmake 0 in
let state = initialize (fun _ _ _ -> ()) (vec v) in
set_stop_time state 1.0;
@ -65,5 +61,5 @@ module InPlace : Csolver.CsolverC =
let copy { state; vec } =
{ state = copy state; vec = Zls.copy vec } in
DNodeC { init; reset; step; copy }
DNodeC { state; reset; step; copy }
end

34
src/lib/solvers/statefulSundials.ml
View file

@ -3,29 +3,26 @@ open Hsim.Types
open Hsim.Solver
open Zls
module Functional : Csolver.Csolver =
module Functional =
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 csolve : (carray, carray) csolver =
let open Cvode in
let init { size; fder=_; _ } =
let vec = cmake size in
let state =
let vec = cmake 0 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 reset { init=i; fder; stop; _ } _ =
let f 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
Cvode.init Adams default_tolerances f 0. (Nvector_serial.wrap i) in
set_stop_time state stop;
{ state; vec } in
{ state; vec=i } in
let step ({ state; vec } as s) h =
let y = Nvector_serial.wrap vec in
@ -33,21 +30,18 @@ module Functional : Csolver.Csolver =
let f t = get_dky state y t 0; Nvector_serial.unwrap y in
(h, f), s in
DNode { init; reset; step }
DNode { state; reset; step }
end
module InPlace : Csolver.Csolver =
module InPlace =
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 csolve : (carray, carray) csolver =
let open Cvode in
let init { size; fder=_; _ } =
let vec = cmake size in
let state =
let vec = cmake 0 in
let state = init Adams default_tolerances (fun _ _ _ -> ()) 0.
(Nvector_serial.wrap vec) in
set_stop_time state 1.0;
@ -66,6 +60,6 @@ module InPlace : Csolver.Csolver =
let f t = get_dky s.state y t 0; Nvector_serial.unwrap y in
(h, f), s in
DNode { init; reset; step }
DNode { state; reset; step }
end

21
src/lib/solvers/statefulZ.ml
View file

@ -3,17 +3,14 @@ open Hsim.Types
open Hsim.Solver
open Zls
module Functional : Zsolver.ZsolverC =
module Functional =
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 zsolve : (carray, zarray, carray) zsolver_c =
let open Illinois in
let init _ =
let state =
{ state = initialize 0 (fun _ _ _ -> ()) (cmake 0);
vec = zmake 0 } in
@ -33,19 +30,17 @@ module Functional : Zsolver.ZsolverC =
let copy s = s in
DNodeC { init; step; reset; copy }
DNodeC { state; step; reset; copy }
end
module InPlace : Zsolver.ZsolverC =
module InPlace =
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 zsolve : (carray, zarray, carray) zsolver_c =
let open Illinois in
let init _ =
let state =
{ state=initialize 0 (fun _ _ _ -> ()) (cmake 0);
vec=zmake 0 } in
@ -65,5 +60,5 @@ module InPlace : Zsolver.ZsolverC =
let copy _ = raise Common.Errors.TODO in
DNodeC { init; step; reset; copy }
DNodeC { state; step; reset; copy }
end

28
src/lib/solvers/zsolver.ml
View file

@ -1,28 +0,0 @@
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)