common.ml 4.29 KB
Newer Older
Benoit Barbot's avatar
Benoit Barbot committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
(*

Copyright (C) 2019 LACL Universite Paris Est Creteil Val de Marne
Authors: Benoît Barbot
Licence: GPL-3.0-or-later

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

*)

Benoit Barbot's avatar
Benoit Barbot committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
let (|>>) x f = match x with
    Some y -> f y
  | None -> None

let (|>>>) x f = match x with
    Some y -> begin try Some (f y) with _ -> None end
  | None -> None

let (|>>|) x v = match x with
    Some y -> y
  | None -> v

let (|<) x f = let () = f x in x
let (|<>|) f (x,y) = f x y


module StringSet = Set.Make(String)

let check_time =
  let init_timer = ref (Sys.time ()) in
  fun () ->
  let t2 = Sys.time () in
  let c = t2 -. !init_timer in
  init_timer := t2;
  c

module VarSet = struct
  type var = int
  type varset = (string*bool) array
  let var_of_int x =x
  let int_of_var x =x
  let zero_var = -1
end

type bound =
Benoit Barbot's avatar
Benoit Barbot committed
58
59
  Finite of int*VarSet.var
| Const of int
Benoit Barbot's avatar
Benoit Barbot committed
60
61
62
63
64
65
| Infinite

let is_time = function
      Finite (_,0) -> true
    | _ -> false

Benoit Barbot's avatar
Benoit Barbot committed
66

Benoit Barbot's avatar
Benoit Barbot committed
67
let build_bound a b =
Benoit Barbot's avatar
Benoit Barbot committed
68
  if b=0 then Const ( a) else Finite ( a,b)
Benoit Barbot's avatar
Benoit Barbot committed
69

Benoit Barbot's avatar
Benoit Barbot committed
70
module SBound = struct
Benoît Barbot's avatar
Benoît Barbot committed
71
let map clockmap = function
Benoit Barbot's avatar
Benoit Barbot committed
72
73
74
75
76
    Infinite -> Infinite
  | Const a -> Const a
  | Finite (a,b) -> let c = clockmap.(b) in
                    if c <> 0 then Finite (a,c) else Const a

Benoît Barbot's avatar
Benoît Barbot committed
77
let print vs f= function
Benoit Barbot's avatar
Benoit Barbot committed
78
    Infinite -> Format.fprintf f "inf"
Benoit Barbot's avatar
Benoit Barbot committed
79
80
81
  | Const (a) -> Format.fprintf f "%i" a
  | Finite (0,b) -> Format.fprintf f "-%s" (fst vs.(b))
  | Finite (a,b) -> Format.fprintf f "%i-%s" a (fst vs.(b))
Benoît Barbot's avatar
Benoît Barbot committed
82
83
84
85

let eval (_,clocks) = function
      Infinite -> infinity
    | Const a -> float a
Benoit Barbot's avatar
Benoit Barbot committed
86
    | Finite (a,b) -> (float a) -. clocks.(VarSet.int_of_var b-1)
Benoît Barbot's avatar
Benoît Barbot committed
87
end
Benoit Barbot's avatar
Benoit Barbot committed
88

Benoît Barbot's avatar
Benoît Barbot committed
89
(* Copy from Boost/math/tools/roots.hpp thus imperative style *)
90
exception Zero_derivative
Benoit Barbot's avatar
Benoit Barbot committed
91
exception Not_converging of float*float
92
93
let newton_raphson_iterate ?factor:(factor= 1e-10) ?max_iter:(max_iter=20) ?bound f guess_p =
    let (min,max) = match bound with
Benoit Barbot's avatar
Benoit Barbot committed
94
        None -> ref (-.max_float), ref max_float
95
96
97
98
99
100
101
102
103
104
105
      | Some (x,y) -> ref x, ref y in

    let guess = ref guess_p in
    let result = ref (if guess_p= infinity then 1.0 else guess_p) in
    let delta = ref max_float
    and delta1 = ref max_float
    and delta2 = ref max_float in
    try for i =0 to max_iter do
          delta2 := !delta1;
          delta1 := !delta;
          let f0,f1 = f !result in
Benoit Barbot's avatar
Benoit Barbot committed
106
          (*Format.printf "nr x:%g f:%g df:%g minmax: %g;%g@." !result f0 f1 !min !max;*)
107
108
109
          if f0 = 0.0 then raise Exit;
          if f1 = 0.0 then raise Zero_derivative;
          delta := f0 /. f1;
Benoit Barbot's avatar
Benoit Barbot committed
110
          (*Format.printf "delta: %g delta2:%g@." !delta !delta2;*)
111
          if abs_float (!delta*.2.0) > abs_float !delta2 then (
Benoit Barbot's avatar
Benoit Barbot committed
112
            (*Format.printf "test@.";*)
113
114
115
116
117
118
119
            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
            else delta := shift;
            delta1:= 3.0*. !delta;
            delta2 := 3.0 *. !delta
          );
Benoit Barbot's avatar
Benoit Barbot committed
120
          (*Format.printf "test2@.";*)
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
          guess := !result;
          result := !result -. !delta;
          if !result <= !min then (
            delta := 0.5 *. (!guess -. !min);
            result := !guess -. !delta;
            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 Exit
          );
          if !delta >0.0 then max := !guess
          else min := !guess;
          if abs_float (!result *. factor) >= abs_float !delta then raise Exit
        done;
Benoît Barbot's avatar
Benoît Barbot committed
136
137
        let f0,_ = f !result in
        raise (Not_converging (!result,f0) )
138
139
    with
      Exit -> !result