Commit 94535ca6 authored by Benoit Barbot's avatar Benoit Barbot
Browse files

progress

parent 718af3d3
Pipeline #1673 failed with stages
in 40 seconds
......@@ -68,7 +68,7 @@ let is_time = function
let build_bound a b =
if b=0 then Const ( a) else Finite ( a,b)
module SBound = struct
module SBound = struct
let map clockmap = function
Infinite -> Infinite
| Const a -> Const a
......@@ -84,16 +84,15 @@ let print vs f= function
let eval (_,clocks) = function
Infinite -> infinity
| Const a -> float a
| Finite (a,b) -> (float a) -. clocks.(VarSet.int_of_var b-1)
| Finite (a,b) -> (float a) -. clocks.(VarSet.int_of_var b-1)
end
(* Copy from Boost/math/tools/roots.hpp thus imperative style *)
exception Zero_derivative
exception Not_converging of float*float
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
None -> ref min_float, ref max_float
None -> ref (-.max_float), ref max_float
| Some (x,y) -> ref x, ref y in
let guess = ref guess_p in
......@@ -105,10 +104,13 @@ let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?boun
delta2 := !delta1;
delta1 := !delta;
let f0,f1 = f !result in
(*Format.printf "nr x:%g f:%g df:%g minmax: %g;%g@." !result f0 f1 !min !max;*)
if f0 = 0.0 then raise Exit;
if f1 = 0.0 then raise Zero_derivative;
delta := f0 /. f1;
(*Format.printf "delta: %g delta2:%g@." !delta !delta2;*)
if abs_float (!delta*.2.0) > abs_float !delta2 then (
(*Format.printf "test@.";*)
let shift = if !delta>0.0 then (!result -. !min)/. 2.0 else (!result -. !max)/. 2.0 in
if !result <>0.0 && (abs_float shift > abs_float !result) then
delta := copysign (0.9*.(abs_float !result)) !delta
......@@ -116,6 +118,7 @@ let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?boun
delta1:= 3.0*. !delta;
delta2 := 3.0 *. !delta
);
(*Format.printf "test2@.";*)
guess := !result;
result := !result -. !delta;
if !result <= !min then (
......@@ -135,4 +138,3 @@ let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?boun
raise (Not_converging (!result,f0) )
with
Exit -> !result
......@@ -235,15 +235,16 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
let taylor_exp_s ep n =
Poly.fold (fun k v acc -> match k with
| Infinite -> assert false
| Const c ->
let p = P.( **.) (F.of_int c) (P.var (Param.svar)) in
P.( +..) acc (P.( *..) (snd @@ taylor_exp_exp p n) v)
| Finite _ ->
let p= P.( *..) (P.of_bound k) (P.var (Param.svar)) in
P.( +..) acc (P.( *..) (snd @@ taylor_exp_exp p n) v)
) ep P.zero
let cons = Poly.fold (fun k v acc -> match k with
| Infinite -> assert false
| Const c ->
let p = P.( **.) (F.of_int c) (P.var (Param.svar)) in
P.( +..) acc (P.( *..) (snd @@ taylor_exp_exp p n) v)
| Finite _ ->
let p= P.( *..) (P.of_bound k) (P.var (Param.svar)) in
P.( +..) acc (P.( *..) (snd @@ taylor_exp_exp p n) v)
) ep P.zero in
Poly.singleton (Const 0) cons
let apply_bound ep var bound =
......@@ -302,17 +303,17 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
let isvart = if var = Param.tvar then -1.0 else 1.0 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 ->
| 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
let to_float ?smp p =
let open ZoneGraph in
let smpv2 = smpv ?smp () in
......@@ -326,17 +327,17 @@ module Make (P:Polynome.S) (Param:sig val smp:float option val tvar:P.var val sv
| Some f,Const e,Some vf ->
Some (f +. (exp (smpv2*.(float e)))*.vf)
) p (Some 0.0)
let apply_all ?smp p tab =
let _,p = Array.fold_left
(fun (i,pit) v ->
i+1,apply_const ?smp pit (var_of_int i) (F.of_float v))
(1, p) tab in
match to_float ?smp p with
match to_float ?smp p with
Some x -> x
| None -> raise Polynome.Not_fully_evaluated
let fully_apply ?smp p x f =
(*assert (not F.is_exact);*)
let p2 = apply_const ?smp p x f in
......
......@@ -22,7 +22,7 @@ with this program; if not, write to the Free Software Foundation, Inc.,
open Common
module ZGI = ZoneGraphInput.Make (ZoneGraphInput.LinearBound)
module ZGI = ZoneGraphInput.Make (ZoneGraphInput.LinearBound)
(*Building the zone graph *)
let rg = ZGI.input_zone_graph ~verbose:0 Sys.argv.(1)
......@@ -52,16 +52,16 @@ let compute_weight rg =
let compv (w,d) s =
if s=0.0 then
let svar = P.var_of_int (P.nb_var-1) in
let we = Weight.taylor_exp_s w (poly) in
(*Format.printf "vte: %a@." P.print we;*)
let de = Weight.taylor_exp_s d (poly+1) in
Format.printf "vte: %a@." P.print de;
let w0 = P.apply we svar (P.const P.F.zero) in
let d0 = P.apply de svar (P.const P.F.zero) in
match P.to_float w0,P.to_float d0 with
Some v1,Some v2 -> v1,v2
| _ -> failwith "fail to evaluate"
let svar = Weight.var_of_int (P.nb_var-1) in
let we = Weight.taylor_exp_s w (poly) in
(*Format.printf "vte: %a@." P.print we;*)
let de = Weight.taylor_exp_s d (poly+1) in
Format.printf "vte: %a@." Weight.print de;
let w0 = Weight.apply_const ~smp:0.0 we svar (Weight.F.zero) in
let d0 = Weight.apply_const ~smp:0.0 de svar (Weight.F.zero) in
match Weight.to_float w0,Weight.to_float d0 with
Some v1,Some v2 -> v1,v2
| _ -> failwith "fail to evaluate"
else match Weight.to_float ~smp:s w ,Weight.to_float ~smp:s d with
Some v1,Some v2 -> v1,v2
| _ -> failwith "fail to evaluate"
......
......@@ -112,12 +112,18 @@ let _ =
close_out fdot
end;
let comp_exp target s (f,df,ddf) =
match Weight.to_float ~smp:s f ,Weight.to_float ~smp:s df, Weight.to_float ~smp:s ddf with
Some v,Some dv,Some ddv -> -. dv /. v -.target , -.ddv /. v +. dv*.dv /. (v*.v)
| _ -> failwith "fail to evaluate" in
(* Compute s from expected_duration *)
let smp = match !expected_duration with
None -> !frequency
| Some target ->
if !verbose>0 then (
printf "Computing laplace parameters s : %g -> " target; flush stdout;
printf "Computing laplace parameters s : target %g@." target; flush stdout;
print_flush ();
);
let w = weights.(rgpoly.init) in
......@@ -125,22 +131,26 @@ let _ =
let dw = Weight.diff w svar in
let ddw = Weight.diff dw svar in
let compv s =
if s=0.0 then
let svar = P.var_of_int (P.nb_var-1) in
let we = Weight.taylor_exp_s w (!npoly) in
(*Format.printf "vte: %a@." P.print we;*)
let de = Weight.taylor_exp_s dw (!npoly+1) in
Format.printf "vte: %a@." P.print de;
let w0 = P.apply we svar (P.const P.F.zero) in
let d0 = P.apply de svar (P.const P.F.zero) in
match P.to_float w0,P.to_float d0 with
Some v1,Some v2 -> v1,v2
| _ -> failwith "fail to evaluate"
else
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 ~max_iter:100 compv 1.0 in
(*Format.printf "test s: %g@." s;*)
let calc =
if abs_float s <= 0.1 then
let svar = Weight.var_of_int (P.nb_var-1) in
let we = Weight.taylor_exp_s w (!npoly) in
(*Format.printf "vte: %a@." P.print we;*)
let de = Weight.taylor_exp_s dw (!npoly+1) in
let dde = Weight.taylor_exp_s ddw (!npoly+2) in
let w0 = Weight.apply_const ~smp:s we svar (Weight.F.of_float s) in
let d0 = Weight.apply_const ~smp:s de svar (Weight.F.of_float s) in
let dd0 = Weight.apply_const ~smp:s dde svar (Weight.F.of_float s) in
w0,d0,dd0
else w,dw,ddw in
comp_exp target s calc in
let exp0,_ = compv 0.0 in
let expinf,_ = compv 1e30 in
let mexpinf,_ = compv (-100.0) in
if !verbose>0 then printf "Expected time range: [%g; %g; %g]@." (expinf+.target) (exp0+.target) (mexpinf+.target);
let x = Common.newton_raphson_iterate ~max_iter:100 compv (1.0) in
let v,_ = compv x in
if !verbose>0 then printf "%g achieved %g [%.2fs]@." x (v+.target) (check_time ());
Some x
......
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