Skip to content

GitLab

Commit c6f19218 authored by Benoit Barbot's avatar Benoit Barbot
Browse files

improve param estimation

parent 2a28638a
Pipeline #2142 passed with stages
in 2 minutes and 37 seconds
Hide whitespace changes
Inline Side-by-side
src/Semantic.ml
View file @ c6f19218
......@@ -41,20 +41,13 @@ module type S = sig
bool ->
(t, Bt.t) ZoneGraph.t
val compute_s :
val compute_params :
?verbose:int ->
?finish:(float * (float * float * float) * float -> unit) ->
float ->
(t, Bt.t) ZoneGraph.t ->
float option
val compute_z :
?verbose:int ->
?finish:(float * float * float -> unit) ->
?seed:int ->
float ->
float option * float option ->
(t, Bt.t) ZoneGraph.t ->
float option
(float * float * float) option
* (float * (float * float * float) * float) option
end
module Make (Bt : ZoneGraphInput.BoundType) (P : Polynomial.WeightStructure) =
......@@ -162,128 +155,141 @@ struct
else if ddl - dl = 0 then v
else F.mul v F.infinity
let compute_s ?(verbose = 1) ?finish expected_duration rg =
(* Weight in initial location *)
let w = rg.statelist.(rg.init).loc_weight in
let dw = diff w s in
let ddw = diff dw s in
if verbose > 1 then (
Format.printf "\nWeight : %a@." P.print w;
Format.printf "Weight': %a@." P.print dw );
let svarexp = var ~exp:(rg.nb_poly + 1) s in
let wts = svarexp *.. w in
let dwts = svarexp *.. dw in
let ddwts = svarexp *.. ddw in
if verbose > 2 then (
Format.printf "Weight2 : %a@." P.print wts;
Format.printf "Weight2': %a@." P.print dwts );
let we = taylor_exp_s w (rg.nb_poly + 1) in
let de = taylor_exp_s dw (rg.nb_poly + 2) in
let dde = taylor_exp_s ddw (rg.nb_poly + 3) in
if verbose > 3 then (
Format.printf "Weight_taylor : %a@." P.print we;
Format.printf "Weight_taylor': %a@." P.print de );
(* computing limit for s -> -infty and s -> + infty *)
let low_w, high_w = dominating_s w in
let low_dw, high_dw = dominating_s dw in
let liminf = F.to_float @@ comp_lim_inf low_w low_dw in
let limsup = F.to_float @@ comp_lim_inf high_w high_dw in
let compv sval =
if verbose > 3 then Format.printf "Time duration for s: %g -> " sval;
let calc =
if abs_float sval < 1.0 then
(* When s very close to 0 use a taylor expansion to compute the duration*)
if abs_float sval <= 0.01 then
if sval = 0.0 then
(* when s=0 compute limit *)
let _, low_we = dominating_s we in
let _, low_dw = dominating_s de in
let lim = comp_lim_zero_p low_we low_dw in
(*Format.printf "f: %a df: %a ddf:%a@." P.print we P.print de P.print dde;
Format.printf "dominating f: %a df: %a -> %s@." print_dom low_we print_dom low_dw (F.to_string lim);*)
(P.one, P.const (F.neg lim), P.one)
else (* computing taylor exp*)
(we, de, dde)
else (wts, dwts, ddwts)
else (w, dw, ddw)
in
let v, dv = comp_exp expected_duration sval calc in
if verbose > 3 then Format.printf "(%g,%g)@." v dv;
(v, dv)
in
(* computing limit for s -> 0 *)
let exp0, _ = compv 0.0 in
let result =
if exp0 = Float.infinity then
let range = (limsup, Float.infinity, Float.infinity) in
if expected_duration < limsup then Error range
else
let x =
Math.newton_raphson_iterate ~bound:(0.0, Float.infinity)
~max_iter:1000 compv 1.0
let compute_s ?(verbose = 1) target rg =
match target with
| None -> None
| Some expected_duration -> (
(* Weight in initial location *)
let w = rg.statelist.(rg.init).loc_weight in
let dw = diff w s in
let ddw = diff dw s in
if verbose > 1 then (
Format.printf "\nWeight : %a@." P.print w;
Format.printf "Weight': %a@." P.print dw );
let svarexp = var ~exp:(rg.nb_poly + 1) s in
let wts = svarexp *.. w in
let dwts = svarexp *.. dw in
let ddwts = svarexp *.. ddw in
if verbose > 2 then (
Format.printf "Weight2 : %a@." P.print wts;
Format.printf "Weight2': %a@." P.print dwts );
let we = taylor_exp_s w (rg.nb_poly + 1) in
let de = taylor_exp_s dw (rg.nb_poly + 2) in
let dde = taylor_exp_s ddw (rg.nb_poly + 3) in
if verbose > 3 then (
Format.printf "Weight_taylor : %a@." P.print we;
Format.printf "Weight_taylor': %a@." P.print de );
(* computing limit for s -> -infty and s -> + infty *)
let low_w, high_w = dominating_s w in
let low_dw, high_dw = dominating_s dw in
let liminf = F.to_float @@ comp_lim_inf low_w low_dw in
let limsup = F.to_float @@ comp_lim_inf high_w high_dw in
let compv sval =
if verbose > 3 then Format.printf "Time duration for s: %g -> " sval;
let calc =
if abs_float sval < 1.0 then
(* When s very close to 0 use a taylor expansion to compute the duration*)
if abs_float sval <= 0.01 then
if sval = 0.0 then
(* when s=0 compute limit *)
let _, low_we = dominating_s we in
let _, low_dw = dominating_s de in
let lim = comp_lim_zero_p low_we low_dw in
(*Format.printf "f: %a df: %a ddf:%a@." P.print we P.print de P.print dde;
Format.printf "dominating f: %a df: %a -> %s@." print_dom low_we print_dom low_dw (F.to_string lim);*)
(P.one, P.const (F.neg lim), P.one)
else (* computing taylor exp*)
(we, de, dde)
else (wts, dwts, ddwts)
else (w, dw, ddw)
in
let v, _ = compv x in
Ok (x, range, v +. expected_duration)
else
let range = (limsup, exp0 +. expected_duration, liminf) in
if expected_duration >= liminf || expected_duration < limsup then
Error range
else
let x = Math.newton_raphson_iterate ~max_iter:1000 compv 1.0 in
let v, _ = compv x in
Ok (x, range, v +. expected_duration)
in
match result with
| Ok (x, y, z) ->
Option.iter (fun f -> f (x, y, z)) finish;
Some x
| Error (sup, mid, inf) ->
Format.printf "range:[%g; %g; %g]@." sup mid inf;
Format.printf "Target outside range, aborting@.";
exit 1
let compute_z ?(verbose = 1) ?finish ?seed target rgp =
let module State = Sampling.State (P) in
let module SMCMonitor = Sampling.MakeSMCMonitor (Bt) (P) (State) in
let module SMCSampler = Sampling.Make (Bt) (P) (State) (SMCMonitor) in
let make_estimate n boltz =
let rgp2 =
ZoneGraph.copy_map_weight rgp (fun p ->
P.apply_const p P.z (P.F.of_float boltz))
in
let smp = if P.is_exp_poly then compute_s 4.0 rgp2 else None in
let module Samp = SMCSampler (struct
let smp = smp
end) in
if
P.to_float ?smp @@ rgp2.statelist.(rgp2.init).loc_weight
|>>> (fun x -> x < 0.0)
|>>| true
then (0, -.target -. 1.0, 0.0)
else
Samp.sample ~out_style:target ~max_iter:100 ~boltz ?seed
~sampler:Low_disc_sampler.random rgp2
[| (Fixed, Fixed); (Fixed, Fixed) |]
n
in
let b, v =
Math.bisect_increasing ~factor:0.001 ~relax:0.1 ~up_bound:false
~low:(0.0, -.target) (0.0, 1.0) (fun b ->
let n, size, _ = make_estimate 10000 b in
if verbose > 0 then Format.printf "|%i@?" n;
if verbose > 2 then
Format.printf " samples with z:%g -> %g@." b (size +. target);
(size, 0.0))
in
Option.iter (fun f -> f (b, 0.0, v +. target)) finish;
Some b
let v, dv = comp_exp expected_duration sval calc in
if verbose > 3 then Format.printf "(%g,%g)@." v dv;
(v, dv)
in
(* computing limit for s -> 0 *)
let exp0, _ = compv 0.0 in
let result =
if exp0 = Float.infinity then
let range = (limsup, Float.infinity, Float.infinity) in
if expected_duration < limsup then Error range
else
let x =
Math.newton_raphson_iterate ~bound:(0.0, Float.infinity)
~max_iter:1000 compv 1.0
in
let v, _ = compv x in
Ok (x, range, v +. expected_duration)
else
let range = (limsup, exp0 +. expected_duration, liminf) in
if expected_duration >= liminf || expected_duration < limsup then
Error range
else
let x = Math.newton_raphson_iterate ~max_iter:1000 compv 1.0 in
let v, _ = compv x in
Ok (x, range, v +. expected_duration)
in
match result with
| Ok (x, y, z) -> Some (x, y, z)
| Error (sup, mid, inf) ->
Format.printf "range:[%g; %g; %g]@." sup mid inf;
Format.printf "Target outside range, aborting@.";
exit 1 )
let compute_params ?(verbose = 1) ?seed targets rgp =
match targets with
| Some t_length, t_duration ->
let module State = Sampling.State (P) in
let module SMCMonitor = Sampling.MakeSMCMonitor (Bt) (P) (State) in
let module SMCSampler = Sampling.Make (Bt) (P) (State) (SMCMonitor) in
let make_estimate n boltz =
let rgp2 =
ZoneGraph.copy_map_weight rgp (fun p ->
P.apply_const p P.z (P.F.of_float boltz))
in
let smp_full =
if P.is_exp_poly then compute_s ~verbose t_duration rgp2 else None
in
let smp = smp_full |>>> fun (x, _, _) -> x in
let module Samp = SMCSampler (struct
let smp = smp
end) in
if
P.to_float ?smp @@ rgp2.statelist.(rgp2.init).loc_weight
|>>> (fun x -> x < 0.0)
|>>| true
then (0, -.t_length -. 1.0, 0.0, smp_full)
else
let a, b, c =
Samp.sample ~out_style:t_length ~max_iter:100 ~boltz ?seed
~sampler:Low_disc_sampler.random rgp2
[| (Fixed, Fixed); (Fixed, Fixed) |]
n
in
(a, b, c, smp_full)
in
let b, v, smp_full =
Math.bisect_increasing ~factor:0.001 ~relax:0.1 ~up_bound:false
~low:(0.0, -.t_length, None) (0.0, 1.0) (fun b ->
let n, size, _, smp = make_estimate 10000 b in
if verbose > 0 then Format.printf "|%i@?" n;
if verbose > 2 then
Format.printf " samples with z:%g -> %g@." b (size +. t_length);
(size, 0.0, smp))
in
(Some (b, 0.0, v +. t_length), smp_full)
| None, t_duration ->
let smp_full =
if P.is_exp_poly then compute_s ~verbose t_duration rgp else None
in
(None, smp_full)
end
src/Semantic.mli
View file @ c6f19218
......@@ -18,20 +18,13 @@ module type S = sig
bool ->
(t, Bt.t) ZoneGraph.t
val compute_s :
val compute_params :
?verbose:int ->
?finish:(float * (float * float * float) * float -> unit) ->
float ->
(t, Bt.t) ZoneGraph.t ->
float option
val compute_z :
?verbose:int ->
?finish:(float * float * float -> unit) ->
?seed:int ->
float ->
float option * float option ->
(t, Bt.t) ZoneGraph.t ->
float option
(float * float * float) option
* (float * (float * float * float) * float) option
end
module Make (Bt : ZoneGraphInput.BoundType) (P : Polynomial.WeightStructure) :
......
src/compute.ml
View file @ c6f19218
......@@ -129,50 +129,66 @@ module Instantiate
(if print_rg && Param.verbose > 4 then Some input_file_name else None)
rg Param.npoly (Param.z_param <> No)
(* Compute parameters z and s *)
let compute_params ?seed rgpoly =
(* Compute z from expected size *)
let boltz =
let z_target =
match Param.z_param with
| No -> None
| Fixed f -> Some f
| Fixed _ -> None
| ToCompute target ->
if Param.verbose > 0 then
Format.printf "Computing boltzmann parameter z for E[N]=%g @?"
target;
compute_z ~verbose:Param.verbose
?seed:
(match seed with None -> None | Some s -> Some ((s * 11) + 7))
~finish:(fun (c, conv_rad, nt) ->
if Param.verbose > 0 then (
if conv_rad > 0.0 then
Format.printf "radius of convergence: %g" conv_rad;
Format.printf " -> z=%g;E[N]=%g [%.2fs]@." c nt
(Common.check_time ()) ))
target rgpoly
Some target
in
Option.iter
(fun b ->
ZoneGraph.update_weight rgpoly (fun p -> apply_const p z (F.of_float b)))
boltz;
(* Compute s from expected_duration *)
let smp =
let smp_target =
match Param.s_param with
| No -> None
| Fixed f -> Some f
| Fixed _ -> None
| ToCompute target ->
if Param.verbose > 0 then
Format.printf "Computing laplace parameter s for E[T]=%g @?" target;
compute_s ~verbose:Param.verbose
~finish:(fun (x, (low, mid, up), nt) ->
if Param.verbose > 0 then (
Format.printf "range:[%g; %g; %g]" low mid up;
Format.printf " -> s=%g;E[T]=%g [%.2fs]@." x nt
(Common.check_time ()) ))
target rgpoly
Some target
in
let boltz_full, smp_full =
compute_params ~verbose:Param.verbose
?seed:(match seed with None -> None | Some s -> Some ((s * 11) + 7))
(z_target, smp_target) rgpoly
in
let boltz =
match (Param.z_param, boltz_full) with
| No, _ -> None
| Fixed f, _ -> Some f
| ToCompute _, Some (c, conv_rad, nt) ->
if Param.verbose > 0 then (
if conv_rad > 0.0 then
Format.printf "radius of convergence: %g" conv_rad;
Format.printf " -> z=%g;E[N]=%g [%.2fs]@." c nt
(Common.check_time ()) );
Some c
| _ -> None
in
let smp =
match (Param.s_param, smp_full) with
| No, _ -> None
| Fixed f, _ -> Some f
| ToCompute _, Some (x, (low, mid, up), nt) ->
if Param.verbose > 0 then (
Format.printf "range:[%g; %g; %g]" low mid up;
Format.printf " -> s=%g;E[T]=%g [%.2fs]@." x nt
(Common.check_time ()) );
Some x
| _ -> None
in
Option.iter
(fun b ->
ZoneGraph.update_weight rgpoly (fun p -> apply_const p z (F.of_float b)))
boltz;
( match smp with
| Some sparam when not F.is_exact ->
ZoneGraph.update_weight rgpoly (fun p ->
......
src/math.ml
View file @ c6f19218
......@@ -57,18 +57,22 @@ let newton_raphson_iterate ?(factor = 1e-9) ?(max_iter = 20) ?bound f guess_p =
let rec bisect_increasing ?(factor = 1e-10) ?(relax = 0.0) ?low
?(up_bound = true) (bmin, bmax) f_to_evaluate =
(* fmin use to detect vertical asymptote *)
let xfmin, yfmin =
match low with None -> (bmin, fst @@ f_to_evaluate bmin) | Some v -> v
let xfmin, yfmin, a_val =
match low with
| None ->
let yf, _, a = f_to_evaluate bmin in
(bmin, yf, a)
| Some v -> v
in
(*Format.printf "bisect [%g; %g] low:%g, %g @?" bmin bmax xfmin yfmin;*)
if up_bound then
if bmax -. bmin < factor then (bmin, yfmin)
if bmax -. bmin < factor then (bmin, yfmin, a_val)
else
let m = 0.5 *. (bmin +. bmax) in
(*if bmax -. bmin <= factor then m*)
let fx, _ = f_to_evaluate m in
let fx, _, a = f_to_evaluate m in
(*Format.printf "eval f(%g)=%g @?" m fx;*)
(*Printf.printf "dicho min:%g; max :%g f:%g\n" bmin bmax fx;*)
(*Detect asymtote *)
......@@ -79,7 +83,7 @@ let rec bisect_increasing ?(factor = 1e-10) ?(relax = 0.0) ?low
(cf = FP_nan || cf = FP_infinite) || fx < yfmin
in
let width = bmax -. bmin in
if abs_float fx < factor || width < factor /. 10. then (*stop *) (m, fx)
if abs_float fx < factor || width < factor /. 10. then (*stop *) (m, fx, a)
else if detect_asymptote || fx > 0.0 then
(* go left *)
bisect_increasing ~relax ~factor ?low
......@@ -87,11 +91,11 @@ let rec bisect_increasing ?(factor = 1e-10) ?(relax = 0.0) ?low
f_to_evaluate
else
(* go right *)
bisect_increasing ~relax ~factor ~low:(m, fx)
bisect_increasing ~relax ~factor ~low:(m, fx, a)
(m -. (relax *. width), bmax +. (relax *. width))
f_to_evaluate
else
let fx, _ = f_to_evaluate bmax in
let fx, _, _ = f_to_evaluate bmax in
if classify_float fx = FP_nan || fx > 0.0 || fx < yfmin then
(* Go left true upperbound *)
bisect_increasing ~factor ~relax ~up_bound:true ?low (bmin, bmax)
......
src/rational_fraction.ml
View file @ c6f19218
......@@ -262,7 +262,7 @@ module Compute (Bt : ZoneGraphInput.BoundType) (P : Polynomial.S) = struct
Option.iter (fun f -> f (x, first_pole, v)) finish;
Some x
Some (x, first_pole, v)
let compute ?update:_ ?finish ?print_interm:_ rg _ _ =
let rgp = ZoneGraph.copy 1 (Bt.zero, Bt.map) zero print rg in
......@@ -270,7 +270,12 @@ module Compute (Bt : ZoneGraphInput.BoundType) (P : Polynomial.S) = struct
Option.iter (fun f -> f rgp) finish;
rgp
let compute_s ?verbose:_ ?finish:_ _ _ = None
(*let compute_s ?verbose:_ ?finish:_ _ _ = None*)
let compute_params ?(verbose = 1) ?seed (t_length, _) rg =
match t_length with
| Some t -> (compute_z ~verbose ?seed t rg, None)
| _ -> (None, None)
end
(*
......
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