Commit 467b72e9 authored by Benoît Barbot's avatar Benoît Barbot
Browse files

pb float

parent 46b83a77
Pipeline #1654 failed with stages
in 35 seconds
......@@ -42,16 +42,19 @@ let sampler = ref Low_disc_sampler.random
let sampling_meth = ref Exact
let gnuplot_driver = ref None
let export_zone_graph = ref None
let frequency = ref 0.0
let frequency = ref None
let max_iteration = ref 20
let rational_impl = ref false
let expected_duration = ref None
let spec_short = [
"--debug", Arg.Set print_rg, "Print the intermediate representation";
"--poly",Arg.Set_int npoly, "Number of polynome Iteration (default 3)";
"--traj",Arg.Set_int nbtraj, "Number of trajectories (default 10)";
"--receding", Arg.Int (fun x -> sampling_meth := Receding x), "trajectory length";
"--frequency", Arg.Set_float frequency, "frequency parameter s (default 0)";
"--frequency", Arg.Float (fun f ->frequency := Some f), "frequency parameter s (default 0)";
"--time-duration", Arg.Float (fun f -> expected_duration:= Some f), "time word duration parameter";
(*"-o", Arg.String (fun s -> outfile := Format.formatter_of_out_channel @@ open_out s), "Output file for trajectories";*)
"--template", Arg.String (fun s -> try
sampling_meth := Driven (timeword_parser s);
......
......@@ -46,7 +46,7 @@ module type WeightStructure = sig
val apply_const : t -> var -> fl -> t
val to_float : ?smp:float -> t -> float option
val apply_all : ?smp:float -> t -> float array -> float
val fully_apply : t -> var -> fl -> fl
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : t -> var -> float -> float
end
......@@ -75,16 +75,18 @@ sig
val init : (int,Bt.t) ZoneGraph.t -> state
val print_state : Format.formatter -> state -> unit
val eval_poly_state : state -> P.t -> P.t
val eval_poly_state_full : state -> P.t -> float
val eval_poly_state_full : ?smp:float -> state -> P.t -> float
val elapse_time : state -> float -> state
val fire : state -> ('b,Bt.t) ZoneGraph.transition -> state
val sample_time :
?smp:float ->
?max_iter:int ->
int ->
state -> (P.t,Bt.t) ZoneGraph.transition -> float -> float
val available_trans :
('a,Bt.t) ZoneGraph.t -> int * 'b -> ('a,Bt.t) ZoneGraph.transition list
val sample_trans :
?smp:float ->
(P.t,Bt.t) ZoneGraph.t ->
int ->
state ->
......@@ -92,6 +94,7 @@ sig
?filter:(float * float -> (P.t,Bt.t) ZoneGraph.transition -> bool) option ->
float -> (P.t,Bt.t) ZoneGraph.transition
val sample_traj :
?smp:float ->
max_iter:int ->
update_callback:(unit -> 'a) * ((state -> int -> 'b) *
(state -> int -> 'c) *
......@@ -154,7 +157,6 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
done;
!stateweight,rgp
exception Break
(* Copy from Boost/math/tools/roots.hpp *)
let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?bound p ?diffp x guess_p =
......@@ -175,7 +177,7 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
delta1 := !delta;
let f0 = fully_apply_float p x !result
and f1 = fully_apply_float pp x !result in
if f0 = 0.0 then raise Break;
if f0 = 0.0 then raise Exit;
if f1 = 0.0 then raise Zero_derivative;
delta := f0 /. f1;
if abs_float (!delta*.2.0) > abs_float !delta2 then (
......@@ -191,20 +193,20 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
if !result <= !min then (
delta := 0.5 *. (!guess -. !min);
result := !guess -. !delta;
if !result = !min || !result = !max then raise Break
if !result = !min || !result = !max then raise Exit
) else if !result >= !max then (
delta := 0.5 *. (!guess -. !max);
result := !guess -. !delta;
if !result = !min || !result = !max then raise Break
if !result = !min || !result = !max then raise Exit
);
if !delta >0.0 then max := !guess
else min := !guess;
if abs_float (!result *. factor) >= abs_float !delta then raise Break
if abs_float (!result *. factor) >= abs_float !delta then raise Exit
done;
Format.eprintf "Alert Not converging %a --> %g@." P.print p !result;
!result
with
Break -> !result
Exit -> !result
let rec find_root_old ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?bound p ?diffp x guess =
......@@ -250,8 +252,8 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
i+1,apply_const pit (var_of_int i) (F.of_float v))
(1, p) tab in
p(*apply_const p (var_of_int j) (1.0/.smp)*)
let eval_poly_state_full (_,tab) p=
P.apply_all p tab
let eval_poly_state_full ?smp (_,tab) p=
P.apply_all ?smp p tab
(* match to_float @@ eval_poly_state state p with
None -> raise Polynome.Not_fully_evaluated
| Some x -> x*)
......@@ -269,7 +271,7 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
(edge.target,tab2)
(* Sample firing time of transition trans using random number u*)
let sample_time ?max_iter:(max_iter=20) i state trans u =
let sample_time ?smp ?max_iter:(max_iter=20) i state trans u =
let low = Bt.eval_max_bound state trans.lower_bound
and up = Bt.eval_min_bound state trans.upper_bound in
assert (low <= up);
......@@ -307,7 +309,7 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
state.transition
(* sample next transition from state using random number u *)
let sample_trans rg i (id,tab) ?(tolerance=1.0e-9) ?(filter=None) u =
let sample_trans ?smp rg i (id,tab) ?(tolerance=1.0e-9) ?(filter=None) u =
let state = rg.statelist.(id) in
let totalw,accl =
......@@ -336,7 +338,7 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
1.0, trans
else
let weight,_,_ = trans.weight.(i) in
let w = eval_poly_state_full (id,tab) weight in
let w = eval_poly_state_full ?smp (id,tab) weight in
if w< -. tolerance (* should be <0 but account for numerical error *) then begin
Format.eprintf "Error: Negative weight : (%a->[%c]) %a->%g@." print_state (id,tab) trans.action P.print weight w;
assert false
......@@ -366,7 +368,7 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
| _ -> snd @@ List.nth accl (min (int_of_float (u*. (float n))) (n-1))
let sample_traj ~max_iter ~update_callback:(up0,(up1,up2,up3)) rgpoly template_word ~samplers:(samp1,samp2) =
let sample_traj ?smp ~max_iter ~update_callback:(up0,(up1,up2,up3)) rgpoly template_word ~samplers:(samp1,samp2) =
let lengthword = Array.length template_word in
(* Initialisation *)
let state = ref (init rgpoly) in
......@@ -384,15 +386,15 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
in
(* Sample the transition *)
let tr = match snd template_word.(j) with
| Rand u -> sample_trans rgpoly polyindex !state ~filter u
| Fixed | Label _ -> sample_trans rgpoly polyindex !state ~filter (samp1 j)
| Rand u -> sample_trans ?smp rgpoly polyindex !state ~filter u
| Fixed | Label _ -> sample_trans ?smp rgpoly polyindex !state ~filter (samp1 j)
in
(* Sample time *)
let time = match fst template_word.(j) with
Label t -> t
| Rand u -> sample_time ~max_iter polyindex !state tr u
| Fixed -> sample_time ~max_iter polyindex !state tr (samp2 j)
| Rand u -> sample_time ?smp ~max_iter polyindex !state tr u
| Fixed -> sample_time ?smp ~max_iter polyindex !state tr (samp2 j)
in
(* Update state *)
state := elapse_time !state time;
......
......@@ -24,7 +24,7 @@ module type WeightStructure =
val apply_const : t -> var -> fl -> t
val to_float : ?smp:float -> t -> float option
val apply_all : ?smp:float -> t -> float array -> float
val fully_apply : t -> var -> fl -> fl
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : t -> var -> float -> float
end
type 'a transChoice = Label of 'a | Rand of float | Fixed
......@@ -52,16 +52,17 @@ sig
val init : (int,Bt.t) ZoneGraph.t -> state
val print_state : Format.formatter -> state -> unit
val eval_poly_state : state -> P.t -> P.t
val eval_poly_state_full : state -> P.t -> float
val eval_poly_state_full : ?smp:float -> state -> P.t -> float
val elapse_time : state -> float -> state
val fire : state -> ('b,Bt.t) ZoneGraph.transition -> state
val sample_time :
val sample_time : ?smp:float ->
?max_iter:int ->
int ->
state -> (P.t,Bt.t) ZoneGraph.transition -> float -> float
val available_trans :
('a,Bt.t) ZoneGraph.t -> int * 'b -> ('a,Bt.t) ZoneGraph.transition list
val sample_trans :
?smp:float ->
(P.t,Bt.t) ZoneGraph.t ->
int ->
state ->
......@@ -69,6 +70,7 @@ sig
?filter:(float * float -> (P.t,Bt.t) ZoneGraph.transition -> bool) option ->
float -> (P.t,Bt.t) ZoneGraph.transition
val sample_traj :
?smp:float ->
max_iter:int ->
update_callback:(unit -> 'a) * ((state -> int -> 'b) *
(state -> int -> 'c) *
......
......@@ -40,7 +40,7 @@ let precomputation_exists ?debug filename frequency npoly exact_rational =
let fd = open_in ((Filename.remove_extension filename)^".pcmp") in
let hashmod = Digest.file filename in
let hashpcmp = (input_value fd: Digest.t) in
let freq = (input_value fd: float) in
let freq = (input_value fd: float option) in
let polyn = (input_value fd: int) in
let exact_rat = (input_value fd: bool) in
let rg = input_value fd in
......
......@@ -44,7 +44,7 @@ module type S = sig
val map_var : t -> (var -> var) -> t
val map : (P.t -> P.t) -> t -> t
val elapse_time : t -> var -> t
val fully_apply : t -> var -> fl -> fl
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : t -> var -> float -> float
val apply_const : t -> var -> fl -> t
val apply_bound: t-> var -> bound -> t
......@@ -347,10 +347,10 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
| None -> raise Polynome.Not_fully_evaluated
let fully_apply p x f =
let fully_apply ?smp p x f =
assert (not F.is_exact);
let p2 = apply_const p x f in
match to_float p2 with
match to_float ?smp p2 with
Some x -> F.of_float x
| None -> raise Polynome.Not_fully_evaluated
......
open Format
module Make (FunIt:Semantic.S) =struct
module Make (FunIt:Semantic.S) (Param:sig val smp: float option end) =struct
module OF = OutFormat.Make(FunIt)
let sample ?(outfile=Format.std_formatter) ?(verbose=1) ~out_style ~max_iter ~sampler ?(store_traj=false) rgpoly template_word nbtraj =
......@@ -26,6 +26,7 @@ module Make (FunIt:Semantic.S) =struct
(try for i =0 to nbtraj-1 do
OF.reset_traj ();
FunIt.sample_traj
?smp:Param.smp
~max_iter
~update_callback:(print_progress,print_callback)
~samplers:((fun j -> Low_disc_sampler.random j i),(fun j -> sampler j i))
......
......@@ -59,7 +59,7 @@ module type S = sig
val primitive : t -> var -> t
val integral : t -> var -> t -> t -> t
val iter_on_var : (var -> 'a -> 'a) -> 'a -> 'a
val fully_apply : t -> var -> fl -> fl
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : t -> var -> float -> float
val apply_bound: t-> var -> bound -> t
val apply_all : ?smp:float -> t -> float array -> float
......@@ -297,7 +297,7 @@ module Make (F: Fl.FSIG)(K:sig val var_string:VarSet.varset end) = struct
if snd var_string.(int_of_var c) then apply acc c (var t +.. (var c)) else acc )
p
let fully_apply p x f =
let fully_apply ?smp p x f =
match to_const @@ apply_const p x f with
Some x -> x
| None -> raise Not_fully_evaluated
......
......@@ -54,7 +54,7 @@ module type S = sig
val primitive : t -> var -> t
val integral : t -> var -> t -> t -> t
val iter_on_var : (var -> 'a -> 'a) -> 'a -> 'a
val fully_apply : t -> var -> fl -> fl
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : t -> var -> float -> float
val apply_bound: t-> var -> Common.bound -> t
val apply_all : ?smp:float -> t -> float array -> float
......
......@@ -70,9 +70,10 @@ let _ =
Some a,_ -> a.ZoneGraph.var_string
| _,Some b -> b.ZoneGraph.var_string
| None,None -> assert false end) in
let module Weight = (val if !frequency =0.0 then (module P:Semantic.WeightStructure)
let module Weight = (val if !frequency = None && !expected_duration = None then (module P:Semantic.WeightStructure)
else (module ExpPoly.Make(P)(struct
let smp = if P.F.is_exact then None else Some !frequency
let smp = if P.F.is_exact || !expected_duration <> None
then None else !frequency
let tvar = P.var_of_int 0
let svar= P.var_of_int (P.nb_var-1)
end))) in
......@@ -105,7 +106,7 @@ let _ =
close_out fdot
end;
let module M = MainLoop.Make(FunIt) in
let module M = MainLoop.Make(FunIt) (struct let smp = Some 1.0 end) in
(* Sampling *)
(* Compute template if not provided *)
......@@ -150,7 +151,7 @@ let _ =
fprintf gnuplot "plot '%s' i 0 u (binwidth*(floor(($1-binstart)/binwidth)+0.5)+binstart):(1.0) smooth freq w boxes\n" data
) else if !store_traj then (
fprintf gnuplot "set output '%s_%i_%g.%s'\n" (Filename.remove_extension data)
!npoly !frequency driver;
!npoly (!frequency |>>| 0.0) driver;
fprintf gnuplot "set xrange [-0.5 :*]\nbinwidth = 0.1\nbinstart= -0.5\nset boxwidth 0.9*binwidth\nset style fill solid 0.5\n";
fprintf gnuplot "plot '%s' i 0 u (binwidth*(floor(($2-binstart)/binwidth)+0.5)+binstart):(1.0) smooth freq w boxes\n" data
) else if !out_style = StateListFull then fprintf gnuplot "plot '%s' u 2:3:4 w points pt 7 lc variable" data
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment