Commit 1dd4904e authored by Benoît Barbot's avatar Benoît Barbot
Browse files

progress

parent 18afcd02
Pipeline #1662 failed with stages
in 34 seconds
......@@ -46,6 +46,7 @@ let frequency = ref None
let max_iteration = ref 20
let rational_impl = ref false
let expected_duration = ref None
let random_seed = ref None
let spec_short = [
......@@ -53,6 +54,7 @@ let spec_short = [
"--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";
"--seed", Arg.Int (fun i -> random_seed := Some i), "seed of the random sampler";
"--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";*)
......
......@@ -46,7 +46,6 @@ module type WeightStructure = sig
val apply_const : ?smp:float -> t -> var -> fl -> t
val to_float : ?smp:float -> t -> float option
val apply_all : ?smp:float -> t -> float array -> float
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : ?smp:float -> t -> var -> float -> float
end
......@@ -158,13 +157,21 @@ module Make (Bt: ZoneGraphInput.BoundType) (P: WeightStructure) = struct
!stateweight,rgp
let find_root ?smp ?factor ?max_iter ?bound p ?diffp x guess_p =
(*let p2 = match smp with
None -> p
| Some s -> apply_const ?smp p (P.var_of_int (P.nb_var-1)) (F.of_float s) in
Format.eprintf "finding root for %a @." P.print p; *)
let pp = match diffp with
None -> diff p x
| Some pp -> pp in
let f v =
let f0 = fully_apply_float ?smp p x v
and f1 = fully_apply_float ?smp pp x v in (f0,f1) in
newton_raphson_iterate ?factor ?max_iter ?bound f guess_p
and f1 = fully_apply_float ?smp pp x v in
(f0,f1) in
try newton_raphson_iterate ?factor ?max_iter ?bound f guess_p with
| Common.Not_converging (r,v) ->
Format.eprintf "Alert Not converging %a --> %g -> %g@." P.print p r v;
r
(* State instantiation for simulation
......
......@@ -24,7 +24,6 @@ module type WeightStructure =
val apply_const : ?smp:float -> t -> var -> fl -> t
val to_float : ?smp:float -> t -> float option
val apply_all : ?smp:float -> t -> float array -> float
val fully_apply : ?smp:float -> t -> var -> fl -> fl
val fully_apply_float : ?smp:float -> t -> var -> float -> float
end
type 'a transChoice = Label of 'a | Rand of float | Fixed
......
......@@ -87,8 +87,9 @@ let eval (_,clocks) = function
| Finite (a,b) -> (float a) -. clocks.(VarSet.int_of_var b-1)
end
(* Copy from Boost/math/tools/roots.hpp *)
(* Copy from Boost/math/tools/roots.hpp thus imperative style *)
exception Zero_derivative
exception Not_converging of float*float
let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?bound f guess_p =
let (min,max) = match bound with
......@@ -130,8 +131,8 @@ let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?boun
else min := !guess;
if abs_float (!result *. factor) >= abs_float !delta then raise Exit
done;
Format.eprintf "Alert Not converging --> %g@." !result;
!result
let f0,_ = f !result in
raise (Not_converging (!result,f0) )
with
Exit -> !result
......@@ -43,6 +43,7 @@ module SBound : sig
val print : VarSet.varset -> Format.formatter -> bound -> unit
val eval: 'a * float array -> bound -> float
end
exception Not_converging of float*float
val newton_raphson_iterate:
?factor:float ->
?max_iter:int ->
......
......@@ -145,6 +145,7 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
let vector_space = (zero, (+..), fun x y -> (F.of_float x) **. y )
let const x =
(* Const 0 means e-0 i.e. 1.0 *)
Poly.singleton (Const 0) (P.const x)
(* let exppoly p var =
......@@ -269,56 +270,45 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
add_mon (Const (0)) v3 acc
| Finite (const_bound,var_bound),Some s ->
assert (not F.is_exact);
let v3 = P.( **.) (F.of_float @@ exp (s*. (float (isvart *const_bound)))) v2 in
let v3 = P.( **.) (F.of_float @@ exp (s *. (float (isvart *const_bound)))) v2 in
add_mon (Finite (0,var_bound)) v3 acc)
ep Poly.empty
let apply_const ep var c =
Poly.fold (fun k v acc ->
let v2 = P.apply v var (P.const c) in
match k with
Infinite -> assert false
| Const _ -> add_mon k v2 acc
| Finite (_,var2) when var2<> var -> add_mon k v2 acc
| Finite (a,_) ->
let isvart = if var = Param.tvar then -1.0 else 1.0 in
match Param.smp with
None -> raise Polynome.Not_fully_evaluated
| Some s ->
assert (not F.is_exact);
let v3 = P.( **.) (F.of_float @@ exp (s *. isvart *. (F.to_float c))) v2 in
add_mon (Const (0)) v3 acc
)
ep Poly.empty
let elapse_time ep t =
map (fun p -> P.elapse_time p t ) ep
(* Check consistency of parameter s*)
let smpv ?smp () =
match smp,Param.smp with
None,Some x -> x
None, Some x when F.is_exact -> assert false
| None, Some x -> x
| Some x,None -> x
| Some a, Some b when a= b -> a
| Some a, Some b when a=b -> a
| Some _, Some _ -> failwith "Conflicting value for s parameter"
| None,None -> failwith "No value for s parameter"
let apply_const ?smp ep var c =
let s = match Param.smp,smp with
None,None -> raise Polynome.Not_fully_evaluated
| _,Some s -> s
| Some s,_ -> assert (not F.is_exact); s in
let s = smpv ?smp () in
if var=Param.svar then assert (F.to_float c=s);
Poly.fold (fun k v acc ->
let v2 = P.apply_const ?smp v var c in
match k with
Infinite -> assert false
| Const _ -> add_mon k v2 acc
| Finite (_,var2) when var2<> var -> add_mon k v2 acc
| Finite (a,_) ->
match k,var=Param.svar with
Infinite,_ -> assert false
| Const _,false -> add_mon k v2 acc
| Finite (_,var2),false when var2<> var -> add_mon k v2 acc
| Finite (a,_),false ->
let isvart = if var = Param.tvar then -1.0 else 1.0 in
let v3 = P.( **.) (F.of_float @@ exp (s *. isvart *. (F.to_float c))) v2 in
let v3 = P.( **.) (F.of_float @@ exp (s *. ((float a) +. isvart *. (F.to_float c)))) v2 in
add_mon (Const (0)) v3 acc
| Const a, true ->
let v3 = P.( **.) (F.of_float @@ exp ( s *. (float a))) v2 in
add_mon (Const (0)) v3 acc
| Finite (a,v),true ->
let v3 = P.( **.) (F.of_float @@ exp ( s *. (float a))) v2 in
add_mon (Finite (0,v)) v3 acc
)
ep Poly.empty
......@@ -327,7 +317,6 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
let open ZoneGraph in
let smpv2 = smpv ?smp () in
Poly.fold (fun k v acc ->
if Param.smp=None then
let v2 = P.apply_const v Param.svar (F.of_float smpv2) in
match acc,k,P.to_float v2 with
None,_,_ -> None
......@@ -336,7 +325,7 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
| _,Infinite,_ -> None
| Some f,Const e,Some vf ->
Some (f +. (exp (smpv2*.(float e)))*.vf)
else P.to_float v) p (Some 0.0)
) p (Some 0.0)
let apply_all ?smp p tab =
let _,p = Array.fold_left
......
......@@ -43,12 +43,15 @@ let phi2=1.32471795724474602596090885447809
let phin n =
let n2 = max 2 n in
Common.newton_raphson_iterate ~factor:1e-15 ~max_iter:1000 (fun v -> let nm = float (n2-1) in
(v ** (float n2) -. v -.1.0), nm*.(v** nm)-.1.0) 1.0
(v ** (float n2) -. v -.1.0), (float n2)*.(v** nm)-.1.0) 1.1
let kronecker_state = ref phi
let init n =
Random.self_init ();
let init ?seed n =
begin match seed with
None -> Random.self_init ()
| Some i -> Random.init i;
end;
kronecker_state := phin n
let kronecker_phi base index =
......
......@@ -2,13 +2,13 @@ open Format
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 =
let sample ?(outfile=Format.std_formatter) ?seed ?(verbose=1) ~out_style ~max_iter ~sampler ?(store_traj=false) rgpoly template_word nbtraj =
(* Sampling *)
(* Most of the code handle printing, Most of the work done by FunIt.sample_traj*)
(* Compute template if not provided *)
let lengthword = Array.length template_word in
Low_disc_sampler.init lengthword;
Low_disc_sampler.init ?seed lengthword;
(*printing progress and state *)
let print_bar = not (outfile == Format.std_formatter) && verbose>0 in
......
......@@ -356,6 +356,7 @@ module Make (F: Fl.FSIG)(K:sig val var_string:VarSet.varset end) = struct
*)
let laplace p smp xv =
(* smp = 1/s *)
let x = int_of_var xv in
let n = Poly.fold (fun k _ i -> max k.(x) i) p 0 in
let a = Array.make (n+1) zero in
......
......@@ -128,12 +128,12 @@ let _ =
match Weight.to_float ~smp:s w ,Weight.to_float ~smp:s dw, Weight.to_float ~smp:s ddw with
Some v,Some dv,Some ddv -> -. dv /. v -.target , -.ddv /. v +. dv*.dv /. (v*.v)
| _ -> failwith "fail to evaluate" in
let x = Common.newton_raphson_iterate compv 1.0 in
let x = Common.newton_raphson_iterate ~max_iter:100 compv 1.0 in
if !verbose>0 then printf "%g [%.2fs]@." x (check_time ());
Some x
in
let module M = MainLoop.Make(FunIt) (struct let smp = smp end) in
let module M = MainLoop.Make(FunIt) (struct let smp=smp end) in
(* Sampling *)
(* Compute template if not provided *)
......@@ -149,6 +149,7 @@ let _ =
M.sample
~outfile
?seed: !Arguments.random_seed
~verbose:!verbose
~out_style: !out_style
~max_iter: !max_iteration
......
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