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

feat: somewhat compatible with zelus output

Browse source
This commit is contained in:
Henri Saudubray 2025-06-23 15:48:58 +02:00
parent 589f89c768
commit 6d92261afd
Signed by: hms
GPG key ID: 7065F57ED8856128
19 changed files with 107 additions and 515 deletions
Download patch file Download diff file Expand all files Collapse all files

4
exm/ball.ml
View file

@ -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} |]
let fout _ _ y = of_array [| y.{1}; y.{0} |]
let jump _ = true
let horizon _ = max_float
let cget s = s.lx
@ -38,7 +38,7 @@ let zset s zin = { s with zin }
let step ({ zin; lx; _ } as s) zfalse =
let lx = if zin.{0} = 1l then
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 }
of_array [| s.lx.{1}; s.lx.{0} |], { zin=zfalse; lx; i=false }
let bouncing_ball () : (_, _, _, carray, carray, carray, zarray, carray) hrec =
let yd = cmake csize in

5
exm/dune
View file

@ -1,3 +1,8 @@
(env
(dev
(flags
(:standard -w -a))))
(library
(name examples)
(libraries hsim solvers))

131
exm/zelus/ballz/ballz.ml
View file

@ -1,131 +0,0 @@
(* The Zelus compiler, version 2.2-dev
(2025-06-16-15:24) *)
open Common
open Ztypes
open Solvers
let (+=) r v = r := !r + v
let g = 9.81
let y0 = 50.
let y'0 = 0.
type ball =
{ mutable major : bool;
mutable h : float;
mutable init : bool;
mutable z : (bool, float) zerocrossing;
mutable y' : float continuous;
mutable y : float continuous }
let ball cstate =
let ball_alloc _ =
cstate.cmax <- cstate.cmax + 2;
cstate.zmax <- cstate.zmax + 1;
{ major = false;
h = 42.;
init = false;
z = { zin = false; zout = 1. };
y' = { pos = 42.; der = 0. };
y = { pos = 42.; der = 0. }
} in
let ball_step self (_, ()) =
let cidx = cstate.cindex in let cpos = ref cidx in
let zidx = cstate.zindex in let zpos = ref zidx in
cstate.cindex <- cstate.cindex + 2;
cstate.zindex <- cstate.zindex + 1;
self.major <- cstate.major;
if cstate.major then
for i = cidx to 1 do Zls.set cstate.dvec i 0. done
else begin
self.y'.pos <- Zls.get cstate.cvec !cpos; cpos += 1;
self.y.pos <- Zls.get cstate.cvec !cpos; cpos += 1
end;
let res0, res1 =
let encore = ref false in
if self.init then self.y'.pos <- y'0;
let last_y' = self.y'.pos in
if self.z.zin then begin
encore := true; self.y'.pos <- -0.8 *. last_y'
end;
self.h <- if !encore then 0. else infinity;
if self.init then self.y.pos <- y0;
self.init <- false;
self.z.zout <- -. self.y.pos;
self.y'.der <- -. g;
self.y.der <- self.y'.pos;
self.y.pos, self.y.der
in
cstate.horizon <- min cstate.horizon self.h;
cpos := cidx;
if cstate.major then begin
Zls.set cstate.cvec !cpos self.y'.pos; cpos += 1;
Zls.set cstate.cvec !cpos self.y.pos; cpos += 1;
self.z.zin <- false
end else begin
self.z.zin <- Zls.get_zin cstate.zinvec !zpos; zpos += 1;
zpos := zidx;
Zls.set cstate.zoutvec !zpos self.z.zout; zpos += 1;
Zls.set cstate.dvec !cpos self.y'.der; cpos += 1;
Zls.set cstate.dvec !cpos self.y.der; cpos += 1
end;
Bigarray.(Array1.of_array Float64 c_layout [| res0; res1 |]) in
let ball_reset self = self.init <- true in
Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
type ('f , 'e , 'd , 'c , 'b , 'a) _main =
{ mutable i_73 : 'f ;
mutable major_62 : 'e ;
mutable h_72 : 'd ;
mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a }
let main (cstate_80:Ztypes.cstate) =
let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball
cstate_80 in
let main_alloc _ =
cstate_80.cmax <- (+) cstate_80.cmax 1;
{ major_62 = false ;
h_72 = 42. ;
i_70 = (false:bool) ;
h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. };
i_73 = i_73_alloc () (* continuous *) } in
let main_step self ((time_61:float) , ()) =
((let (cindex_81:int) = cstate_80.cindex in
let cpos_83 = ref (cindex_81:int) in
cstate_80.cindex <- (+) cstate_80.cindex 1 ;
self.major_62 <- cstate_80.major ;
(if cstate_80.major then
for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done
else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ;
cpos_83 := (+) !cpos_83 1))) ;
(let (result_85) =
let h_71 = ref (infinity:float) in
(if self.i_70 then self.h_68 <- (+.) time_61 0.) ;
(let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in
self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ;
h_71 := min !h_71 self.h_68 ;
self.h_72 <- !h_71 ;
self.i_70 <- false ;
self.t_63.der <- 1. ;
(let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in
(begin match z_69 with
| true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64
| _ -> () end) ;
Bigarray.(Array1.create Float64 c_layout 0))) in
cstate_80.horizon <- min cstate_80.horizon self.h_72 ;
cpos_83 := cindex_81 ;
(if cstate_80.major then
(((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ;
cpos_83 := (+) !cpos_83 1)))
else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ;
cpos_83 := (+) !cpos_83 1)))) ; result_85))) in
let main_reset self =
((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ):
unit) in
Node { alloc = main_alloc; step = main_step ; reset = main_reset }

137
exm/zelus/ballz/ballz.ml.bak
View file

@ -1,137 +0,0 @@
(* The Zelus compiler, version 2.2-dev
(2025-06-16-15:24) *)
open Common
open Ztypes
open Solvers
let g = 9.81
let y0 = 50.
let y'0 = 0.
type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball =
{ mutable major_50 : 'g ;
mutable h_60 : 'f ;
mutable h_58 : 'e ;
mutable i_56 : 'd ;
mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a }
let ball (cstate_74:Ztypes.cstate) =
let ball_alloc _ =
cstate_74.cmax <- (+) cstate_74.cmax 2 ;
cstate_74.zmax <- (+) cstate_74.zmax 1;
{ major_50 = false ;
h_60 = 42. ;
h_58 = (42.:float) ;
i_56 = (false:bool) ;
x_55 = { zin = false; zout = 1. } ;
y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in
let ball_step self ((_time_49:float) , ((y0_48:float) , (y'0_47:float))) =
Printf.printf "STEP (%d)\n" cstate_74.cindex;
((let (cindex_75:int) = cstate_74.cindex in
let cpos_77 = ref (cindex_75:int) in
let (zindex_76:int) = cstate_74.zindex in
let zpos_78 = ref (zindex_76:int) in
cstate_74.cindex <- (+) cstate_74.cindex 2 ;
cstate_74.zindex <- (+) cstate_74.zindex 1 ;
self.major_50 <- cstate_74.major ;
(if cstate_74.major then
for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done
else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ;
cpos_77 := (+) !cpos_77 1) ;
(self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ;
cpos_77 := (+) !cpos_77 1))) ;
(let (result_79:float) =
let h_59 = ref (infinity:float) in
let encore_57 = ref (false:bool) in
(if self.i_56 then self.y'_52.pos <- y'0_47) ;
(let (l_54:float) = self.y'_52.pos in
(begin match self.x_55.zin with
| true ->
encore_57 := true ;
self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end)
;
self.h_58 <- (if !encore_57 then 0. else infinity) ;
h_59 := min !h_59 self.h_58 ;
self.h_60 <- !h_59 ;
(if self.i_56 then self.y_51.pos <- y0_48) ;
self.i_56 <- false ;
self.x_55.zout <- (~-.) self.y_51.pos ;
self.y'_52.der <- (~-.) g ;
self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in
cstate_74.horizon <- min cstate_74.horizon self.h_60 ;
cpos_77 := cindex_75 ;
(if cstate_74.major then
(((Printf.printf "idx: %d\n" !cpos_77;
Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ;
cpos_77 := (+) !cpos_77 1) ;
(Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ;
cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false)))
else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ;
zpos_78 := (+) !zpos_78 1)) ;
zpos_78 := zindex_76 ;
((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ;
zpos_78 := (+) !zpos_78 1)) ;
((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ;
cpos_77 := (+) !cpos_77 1) ;
(Zls.set cstate_74.dvec !cpos_77 self.y_51.der ;
cpos_77 := (+) !cpos_77 1)))) ; result_79)):float) in
let ball_reset self =
(self.i_56 <- true:unit) in
Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
type ('f , 'e , 'd , 'c , 'b , 'a) _main =
{ mutable main_ball : 'f ;
mutable main_major : 'e ;
mutable h_72 : 'd; mutable i_70 : 'c; mutable h : 'b;
mutable t_63 : 'a }
let main cs =
let Node
{ alloc = ball_alloc;
step = ball_step;
reset = ball_reset } = ball cs in
let main_alloc _ =
cs.cmax <- cs.cmax + 1;
{ main_major = false;
h_72 = 42.0; i_70 = false; h = 42.0;
t_63 = { pos = 42.; der = 0. };
main_ball = ball_alloc () } in
let main_step self (time, ()) =
let cindex = cs.cindex in
let cpos = ref cindex in
Printf.printf "main:cindex: %d\n" cs.cindex;
cs.cindex <- cs.cindex + 1;
self.main_major <- cs.major;
if cs.major then for i = cindex to 0 do Zls.set cs.dvec i 0. done
else begin self.t_63.pos <- Zls.get cs.cvec !cpos; cpos := !cpos + 1 end;
let result =
if self.i_70 then self.h <- time;
let z = self.main_major && (time >= self.h) in
if z then self.h <- self.h +. 0.01;
self.h_72 <- min infinity self.h;
self.i_70 <- false;
self.t_63.der <- 1.;
let y_64 = ball_step self.main_ball (time, (y0, y'0)) in
if z then begin
print_float self.t_63.pos;
print_string "\t";
print_float y_64;
print_newline ()
end;
Bigarray.(Array1.create Float64 c_layout 0) in
cs.horizon <- min cs.horizon self.h_72;
cpos := cindex;
if cs.major then Zls.set cs.cvec !cpos self.t_63.pos
else Zls.set cs.dvec !cpos self.t_63.der;
result in
let main_reset self =
self.i_70 <- true;
self.t_63.pos <- 0.;
ball_reset self.main_ball in
Node { alloc = main_alloc; step = main_step ; reset = main_reset }

135
exm/zelus/ballz/ballz.ml.bak2
View file

@ -1,135 +0,0 @@
(* The Zelus compiler, version 2.2-dev
(2025-06-16-15:24) *)
open Common
open Ztypes
open Solvers
let g = 9.81
let y0 = 50.
let y'0 = 0.
type ('g , 'f , 'e , 'd , 'c , 'b , 'a) _ball =
{ mutable major_50 : 'g ;
mutable h_60 : 'f ;
mutable h_58 : 'e ;
mutable i_56 : 'd ;
mutable x_55 : 'c ; mutable y'_52 : 'b ; mutable y_51 : 'a }
let ball (cstate_74:Ztypes.cstate) =
let ball_alloc _ =
cstate_74.cmax <- (+) cstate_74.cmax 2 ;
cstate_74.zmax <- (+) cstate_74.zmax 1;
{ major_50 = false ;
h_60 = 42. ;
h_58 = (42.:float) ;
i_56 = (false:bool) ;
x_55 = { zin = false; zout = 1. } ;
y'_52 = { pos = 42.; der = 0. } ; y_51 = { pos = 42.; der = 0. } } in
let ball_step self ((_time_49:float) , ()) =
((let (cindex_75:int) = cstate_74.cindex in
let cpos_77 = ref (cindex_75:int) in
let (zindex_76:int) = cstate_74.zindex in
let zpos_78 = ref (zindex_76:int) in
cstate_74.cindex <- (+) cstate_74.cindex 2 ;
cstate_74.zindex <- (+) cstate_74.zindex 1 ;
self.major_50 <- cstate_74.major ;
(if cstate_74.major then
for i_1 = cindex_75 to 1 do Zls.set cstate_74.dvec i_1 0. done
else ((self.y'_52.pos <- Zls.get cstate_74.cvec !cpos_77 ;
cpos_77 := (+) !cpos_77 1) ;
(self.y_51.pos <- Zls.get cstate_74.cvec !cpos_77 ;
cpos_77 := (+) !cpos_77 1))) ;
(let (result_79:float) =
let h_59 = ref (infinity:float) in
let encore_57 = ref (false:bool) in
(if self.i_56 then self.y'_52.pos <- y'0) ;
(let (l_54:float) = self.y'_52.pos in
begin match self.x_55.zin with
| true ->
encore_57 := true ;
self.y'_52.pos <- ( *. ) (-0.8) l_54 | _ -> () end;
self.h_58 <- (if !encore_57 then 0. else infinity) ;
h_59 := min !h_59 self.h_58 ;
self.h_60 <- !h_59 ;
(if self.i_56 then self.y_51.pos <- y0) ;
self.i_56 <- false ;
self.x_55.zout <- (~-.) self.y_51.pos ;
self.y'_52.der <- (~-.) g ;
self.y_51.der <- self.y'_52.pos ; self.y_51.pos) in
cstate_74.horizon <- min cstate_74.horizon self.h_60 ;
cpos_77 := cindex_75 ;
(if cstate_74.major then
(((Zls.set cstate_74.cvec !cpos_77 self.y'_52.pos ;
cpos_77 := (+) !cpos_77 1) ;
(Zls.set cstate_74.cvec !cpos_77 self.y_51.pos ;
cpos_77 := (+) !cpos_77 1)) ; ((self.x_55.zin <- false)))
else (((self.x_55.zin <- Zls.get_zin cstate_74.zinvec !zpos_78 ;
zpos_78 := (+) !zpos_78 1)) ;
zpos_78 := zindex_76 ;
((Zls.set cstate_74.zoutvec !zpos_78 self.x_55.zout ;
zpos_78 := (+) !zpos_78 1)) ;
((Zls.set cstate_74.dvec !cpos_77 self.y'_52.der ;
cpos_77 := (+) !cpos_77 1) ;
(Zls.set cstate_74.dvec !cpos_77 self.y_51.der ;
cpos_77 := (+) !cpos_77 1)))) ;
Bigarray.(Array1.of_array Float64 c_layout [| result_79 |])))) in
let ball_reset self =
(self.i_56 <- true:unit) in
Node { alloc = ball_alloc; step = ball_step ; reset = ball_reset }
type ('f , 'e , 'd , 'c , 'b , 'a) _main =
{ mutable i_73 : 'f ;
mutable major_62 : 'e ;
mutable h_72 : 'd ;
mutable i_70 : 'c ; mutable h_68 : 'b ; mutable t_63 : 'a }
let main (cstate_80:Ztypes.cstate) =
let Node { alloc = i_73_alloc; step = i_73_step ; reset = i_73_reset } = ball
cstate_80 in
let main_alloc _ =
cstate_80.cmax <- (+) cstate_80.cmax 1;
{ major_62 = false ;
h_72 = 42. ;
i_70 = (false:bool) ;
h_68 = (42.:float) ; t_63 = { pos = 42.; der = 0. };
i_73 = i_73_alloc () (* continuous *) } in
let main_step self ((time_61:float) , ()) =
((let (cindex_81:int) = cstate_80.cindex in
let cpos_83 = ref (cindex_81:int) in
cstate_80.cindex <- (+) cstate_80.cindex 1 ;
self.major_62 <- cstate_80.major ;
(if cstate_80.major then
for i_1 = cindex_81 to 0 do Zls.set cstate_80.dvec i_1 0. done
else ((self.t_63.pos <- Zls.get cstate_80.cvec !cpos_83 ;
cpos_83 := (+) !cpos_83 1))) ;
(let (result_85) =
let h_71 = ref (infinity:float) in
(if self.i_70 then self.h_68 <- (+.) time_61 0.) ;
(let (z_69:bool) = (&&) self.major_62 ((>=) time_61 self.h_68) in
self.h_68 <- (if z_69 then (+.) self.h_68 0.01 else self.h_68) ;
h_71 := min !h_71 self.h_68 ;
self.h_72 <- !h_71 ;
self.i_70 <- false ;
self.t_63.der <- 1. ;
(let (y_64:float) = (i_73_step self.i_73 (time_61 , ())).{0} in
(begin match z_69 with
| true -> Printf.printf "%.10e\t%.10e\n" self.t_63.pos y_64
| _ -> () end) ;
Bigarray.(Array1.create Float64 c_layout 0))) in
cstate_80.horizon <- min cstate_80.horizon self.h_72 ;
cpos_83 := cindex_81 ;
(if cstate_80.major then
(((Zls.set cstate_80.cvec !cpos_83 self.t_63.pos ;
cpos_83 := (+) !cpos_83 1)))
else (((Zls.set cstate_80.dvec !cpos_83 self.t_63.der ;
cpos_83 := (+) !cpos_83 1)))) ; result_85))) in
let main_reset self =
((self.i_70 <- true ; self.t_63.pos <- 0. ; i_73_reset self.i_73 ):
unit) in
Node { alloc = main_alloc; step = main_step ; reset = main_reset }

BIN
exm/zelus/ballz/ballz.zci
View file

Binary file not shown.

11
exm/zelus/ballz/ballz.zls
View file

@ -2,18 +2,21 @@ let g = 9.81
let y0 = 50.0
let y'0 = 0.0
let hybrid ball (y0, y'0) = y where
let hybrid ball (y0, y'0) = (y, y', z) where
rec der y = y' init y0
and der y' = -. g init y'0 reset z -> -0.8 *. (last y')
and z = up(-. y)
let hybrid main () =
let der t = 1.0 init 0.0 in
let y = ball (y0, y'0) in
let z = period(0.01) in
present z -> (
let rec der p = 1.0 init -0.01 reset s -> -0.01
and s = up(p) in
let (y, y', z) = ball (y0, y'0) in
present z | s -> (
print_float t;
print_string "\t";
print_float y;
print_string "\t";
print_float y';
print_newline ()
); ()

6
exm/zelus/ballz/dune Normal file
View file

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

6
exm/zelus/sincos/dune Normal file
View file

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

45
exm/zelus/sincos/sincosz.ml
View file

@ -1,45 +0,0 @@
(* The Zelus compiler, version 2.2-dev
(2025-06-16-15:24) *)
open Common
open Ztypes
open Solvers
type ('c , 'b , 'a) _f =
{ mutable major_11 : 'c ; mutable sin_13 : 'b ; mutable cos_12 : 'a }
let f (cstate_14:Ztypes.cstate) =
let f_alloc _ =
cstate_14.cmax <- (+) cstate_14.cmax 2;
{ major_11 = false ;
sin_13 = { pos = 42.; der = 0. } ; cos_12 = { pos = 42.; der = 0. } } in
let f_step self ((_time_10:float) , ()) =
((let (cindex_15:int) = cstate_14.cindex in
let cpos_17 = ref (cindex_15:int) in
cstate_14.cindex <- (+) cstate_14.cindex 2 ;
self.major_11 <- cstate_14.major ;
(if cstate_14.major then
for i_1 = cindex_15 to 1 do Zls.set cstate_14.dvec i_1 0. done
else ((self.sin_13.pos <- Zls.get cstate_14.cvec !cpos_17 ;
cpos_17 := (+) !cpos_17 1) ;
(self.cos_12.pos <- Zls.get cstate_14.cvec !cpos_17 ;
cpos_17 := (+) !cpos_17 1))) ;
(let (result_19) =
self.cos_12.der <- (~-.) self.sin_13.pos ;
self.sin_13.der <- self.cos_12.pos ;
Bigarray.(Array1.of_array Float64 c_layout
[| self.sin_13.pos; self.cos_12.pos |]) in
cpos_17 := cindex_15 ;
(if cstate_14.major then
(((Zls.set cstate_14.cvec !cpos_17 self.sin_13.pos ;
cpos_17 := (+) !cpos_17 1) ;
(Zls.set cstate_14.cvec !cpos_17 self.cos_12.pos ;
cpos_17 := (+) !cpos_17 1)))
else (((Zls.set cstate_14.dvec !cpos_17 self.sin_13.der ;
cpos_17 := (+) !cpos_17 1) ;
(Zls.set cstate_14.dvec !cpos_17 self.cos_12.der ;
cpos_17 := (+) !cpos_17 1)))) ; result_19))) in
let f_reset self =
((self.sin_13.pos <- 0. ; self.cos_12.pos <- 1.):unit) in
Node { alloc = f_alloc; step = f_step ; reset = f_reset }

BIN
exm/zelus/sincos/sincosz.zci
View file

Binary file not shown.

4
exm/ztypes.ml Normal file
View file

@ -0,0 +1,4 @@
include Common
include Ztypes
include Solvers

1
src/bin/lift.ml
View file

@ -38,7 +38,6 @@ let lift f =
(* the function that compute the zero-crossings *)
let fzer { state; _ } input y =
Common.Debug.print "LIFT :: Calling [fzer]";
cstate.major <- false;
cstate.zinvec <- no_roots_in;
cstate.dvec <- ignore_der;

36
src/bin/main.ml
View file

@ -19,6 +19,7 @@ let maxstep = ref None
let mintol = ref None
let maxtol = ref None
let no_print = ref false
let zelus = ref false
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
@ -44,6 +45,7 @@ let opts = [
"-mintol", Arg.String (opt mintol), "\tSet minimum solver tolerance";
"-maxtol", Arg.String (opt maxtol), "\tSet maximum solver tolerance";
"-no-print", Arg.Set no_print, "\tDo not print output values";
"-zelus", Arg.Set zelus, "\tUse the output of the Zélus compiler";
]
let errmsg = "Usage: " ^ Sys.executable_name ^ " [OPTIONS] MODEL\nOptions are:"
@ -53,18 +55,30 @@ 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 step s () = let _, s = m.step s () in ret, s in
HNode { m with fout; step }
let m =
try match !model with
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ball" -> Ball.init args
| Some "vdp" -> Vdp.init args
| Some "sincos" -> Sincos.init args
| Some "sqrt" -> Sqrt.init args
| Some "sin1x" -> Sin1x.init args
| Some "sin1xd" -> Sin1x_der.init args
| Some "ballz" -> lift Ballz.ball
| Some "sincosz" -> lift Sincosz.f
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
try
if !zelus then
match !model with
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ballz" -> wrap_zelus (lift Ballz.main)
| Some "sincosz" -> wrap_zelus (lift Sincosz.f)
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
else
match !model with
| None -> Format.eprintf "Missing model\n"; exit 2
| Some "ball" -> Ball.init args
| Some "vdp" -> Vdp.init args
| Some "sincos" -> Sincos.init args
| Some "sqrt" -> Sqrt.init args
| Some "sin1x" -> Sin1x.init args
| Some "sin1xd" -> Sin1x_der.init args
| Some s -> Format.eprintf "Unknown model: %s\n" s; exit 2
with Invalid_argument s -> Format.eprintf "%s\n" s; exit 2
let st = if !inplace then (module State.InPlaceSimState : State.SimState)

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

@ -7,12 +7,13 @@ module Sim (S : SimState) =
struct
include S
let step_discrete s step hor fder fzer cget csize zsize jump reset =
Common.Debug.print "SIMU :: DISCRETE :: start";
let ms, ss = get_mstate s, get_sstate s in
let step_discrete s step hor fder fzer cget zset csize zsize jump reset =
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 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
@ -27,30 +28,25 @@ module Sim (S : SimState) =
let mode, stop, now = Continuous, i.h, 0.0 in
update ms ss (set_running ~mode ~input:i ~stop ~now s)
end else set_running ~mode:Continuous s in
Common.Debug.print "SIMU :: DISCRETE :: end";
Utils.dot o, s
let step_continuous s step cset fout zset =
Common.Debug.print "SIMU :: CONTINUOUS :: start";
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 (h, f, z), ss = step ss stop in
let stop = min stop (hor ms) in
let (h, f, z), ss = step ss (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 s, c = match z with
| None ->
let s, c = if h >= stop
then set_running ~mode:Discrete ~now:h s, Discontinuous
else set_running ~now: h s, Continuous in
update ms ss s, c
| Some z ->
let s = set_running ~mode:Discrete ~now:h s in
update (zset ms z) ss s, Discontinuous in
if h >= stop
then set_running ~mode:Discrete ~now:h s, Discontinuous
else set_running ~now: h s, Continuous
| Some _ -> set_running ~mode:Discrete ~now:h s, Discontinuous in
let h = h -. now in
Common.Debug.print "SIMU :: CONTINUOUS :: end";
{ h; u=fout; c }, s, { h; c; u=fms }
{ h; u=fout; c }, update ms ss (set_zin z s), { h; c; u=fms }
(** Simulation of a model with any solver. *)
let run
@ -59,11 +55,11 @@ module Sim (S : SimState) =
: ('p * (('y, 'yder) ivp * ('y, 'zout) zc), 'a, 'b) sim
= let state = get_init m.state s.state in
let step_discrete st =
let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget
let o, s = step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset
m.csize m.zsize m.jump s.reset in
Some o, s in
let step_continuous st =
let o, s, _ = step_continuous st s.step m.cset m.fout m.zset in
let o, s, _ = step_continuous st s.step m.cset m.fout m.horizon in
Some o, s in
let step st = function
@ -95,8 +91,8 @@ module Sim (S : SimState) =
let step_discrete (st, al) =
let m=m.body in
let o, st =
step_discrete st m.step m.horizon m.fder m.fzer m.cget m.csize m.zsize
m.jump s.reset in
step_discrete st m.step m.horizon m.fder m.fzer m.cget m.zset m.csize
m.zsize m.jump s.reset in
let al = List.map (fun (DNode a) ->
let _, state = a.step a.state @@ Some (Utils.dot @@ get_mstate st) in
DNode { a with state }) al in
@ -104,7 +100,7 @@ module Sim (S : SimState) =
let step_continuous (st, al) =
let ({ h; _ } as o), st, u =
step_continuous st s.step m.body.cset m.body.fout m.body.zset in
step_continuous st s.step m.body.cset m.body.fout m.body.horizon in
let al = List.map (fun (DNode a) ->
(* Step assertions repeatedly until they reach the horizon. *)
let rec step s =
@ -140,10 +136,10 @@ module Sim (S : SimState) =
= let state = get_init m.state s.state in
let step_discrete st =
let o, st = step_discrete st m.step m.horizon m.fder m.fzer m.cget
m.csize m.zsize m.jump s.reset in
m.zset m.csize m.zsize m.jump s.reset in
Some o, st in
let step_continuous st =
let o, st, _ = step_continuous st s.step m.cset m.fout m.zset in
let o, st, _ = step_continuous st s.step m.cset m.fout m.horizon in
o, st in
let rec step st = function

54
src/lib/hsim/state.ml
View file

@ -15,59 +15,65 @@ module type SimState =
- Idle: waiting for input;
- Running: currently integrating; in this case, we have access to the
step mode, current input, timestamp and stop time. *)
type ('a, 'ms, 'ss) state
type ('a, 'ms, 'ss, 'zin) state
(** Get the model state. *)
val get_mstate : ('a, 'ms, 'ss) state -> 'ms
val get_mstate : ('a, 'ms, 'ss, 'zin) state -> 'ms
(** Get the solver state. *)
val get_sstate : ('a, 'ms, 'ss) state -> 'ss
val get_sstate : ('a, 'ms, 'ss, 'zin) state -> 'ss
(** Get the last zero-crossing value. *)
val get_zin : ('a, 'ms, 'ss, 'zin) state -> 'zin option
(** Get the current step mode.
⚠ Should only be called when running (see [is_running]). *)
val get_mode : ('a, 'ms, 'ss) state -> mode
val get_mode : ('a, 'ms, 'ss, 'zin) state -> mode
(** Get the current input.
⚠ Should only be called when running (see [is_running]). *)
val get_input : ('a, 'ms, 'ss) state -> 'a value
val get_input : ('a, 'ms, 'ss, 'zin) state -> 'a value
(** Get the current timestamp.
⚠ Should only be called when running (see [is_running]). *)
val get_now : ('a, 'ms, 'ss) state -> time
val get_now : ('a, 'ms, 'ss, 'zin) state -> time
(** Get the current stop time.
⚠ Should only be called when running (see [is_running]). *)
val get_stop : ('a, 'ms, 'ss) state -> time
val get_stop : ('a, 'ms, 'ss, 'zin) state -> time
(** Build an initial state. *)
val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss) state
val get_init : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state
(** Is the simulation running or idle ? *)
val is_running : ('a, 'ms, 'ss) state -> bool
val is_running : ('a, 'ms, 'ss, 'zin) state -> bool
(** Update the model state. *)
val set_mstate : 'ms -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
val set_mstate : 'ms -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the solver state. *)
val set_sstate : 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
val set_sstate : 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the zero-crossing value. *)
val set_zin : 'zin option -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update both the solver and model states. *)
val update : 'ms -> 'ss -> ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
val update : 'ms -> 'ss -> ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the status to running. *)
val set_running :
?mode:mode -> ?input:'a value -> ?now:time -> ?stop:time ->
('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
(** Update the status to idle. *)
val set_idle : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
val set_idle : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
end
module type SimStateCopy =
sig
include SimState
val copy : ('a, 'ms, 'ss) state -> ('a, 'ms, 'ss) state
val copy : ('a, 'ms, 'ss, 'zin) state -> ('a, 'ms, 'ss, 'zin) state
end
module FunctionalSimState : SimState =
@ -88,15 +94,17 @@ module FunctionalSimState : SimState =
(** Internal state of the simulation node: model state, solver state and
current simulation status. *)
type ('a, 'ms, 'ss) state =
type ('a, 'ms, 'ss, 'zin) state =
{ status : 'a status; (** Current simulation status. *)
mstate : 'ms; (** Model state. *)
sstate : 'ss } (** Solver state. *)
sstate : 'ss; (** Solver state. *)
zin : 'zin option; } (** Last zero-crossing vector *)
exception Not_running
let get_mstate state = state.mstate
let get_sstate state = state.sstate
let get_zin state = state.zin
let is_running state =
match state.status with Running _ -> true | Idle -> false
@ -120,6 +128,7 @@ module FunctionalSimState : SimState =
let set_mstate mstate state = { state with mstate }
let set_sstate sstate state = { state with sstate }
let set_zin zin state = { state with zin }
let update mstate sstate state = { state with mstate; sstate }
@ -132,7 +141,7 @@ module FunctionalSimState : SimState =
let get_stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running
let get_init mstate sstate = { status = Idle; mstate; sstate }
let get_init mstate sstate = { status = Idle; mstate; sstate; zin = None }
end
module InPlaceSimState : SimState =
@ -146,15 +155,17 @@ module InPlaceSimState : SimState =
mutable stop : time;
} -> 'a status
type ('a, 'ms, 'ss) state =
type ('a, 'ms, 'ss, 'zin) state =
{ mutable status : 'a status;
mutable mstate : 'ms;
mutable sstate : 'ss }
mutable sstate : 'ss;
mutable zin : 'zin option }
exception Not_running
let get_mstate state = state.mstate
let get_sstate state = state.sstate
let get_zin state = state.zin
let is_running state =
match state.status with Running _ -> true | Idle -> false
@ -179,6 +190,7 @@ module InPlaceSimState : SimState =
let set_mstate mstate state = state.mstate <- mstate; state
let set_sstate sstate state = state.sstate <- sstate; state
let set_zin zin state = state.zin <- zin; state
let update mstate sstate state =
state.mstate <- mstate; state.sstate <- sstate; state
@ -192,6 +204,6 @@ module InPlaceSimState : SimState =
let get_stop s =
match s.status with Running r -> r.stop | Idle -> raise Not_running
let get_init mstate sstate = { status = Idle; mstate; sstate }
let get_init mstate sstate = { status = Idle; mstate; sstate; zin=None }
end

1
src/lib/solvers/illinois.ml
View file

@ -154,7 +154,6 @@ 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 =
Common.Debug.print "ZSOL :: Calling [step]";
(* swap f0 and f1; f0 takes the previous value of f1 *)
s.f0 <- f1;
s.t0 <- t1;

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

@ -22,8 +22,6 @@ module Functional =
{ state; vec = init } in
let step ({ state ; vec=v } as s) h =
Common.Debug.print "SOLVER STEP";
Common.Debug.print_entry v;
let y_nv = vec v in
let h = step state h y_nv in
let state = copy state in

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

@ -17,8 +17,6 @@ module Functional =
let reset { fzer; init; size } { vec; _ } =
let fzer t cvec zout =
let zout' = fzer t cvec in blit zout' zout in
Common.Debug.print "ZSolver Reset";
Common.Debug.print_entry init;
{ state = initialize size fzer init;
vec = if length vec = size then vec else zmake size } in