examples_in_progress/convex_hull.mlw

   1
   2
   3 theory ConvexSet
   4
   5 use export int.Int
   6 use export real.RealInfix
   7
   8 type pt = { x : real ; y : real }
   9
  10 function crossproduct (p1 p2 p3:pt) : real =
  11   (p2.x -. p1.x) *. (p3.y -. p1.y) -. (p3.x -. p1.x) *. (p2.y -. p1.y)
  12
  13 predicate colinear (p1 p2 p3:pt) = crossproduct p1 p2 p3 = 0.0
  14
  15 predicate turns_left (p1 p2 p3:pt) = crossproduct p1 p2 p3 >. 0.0
  16
  17 lemma knuth1:
  18   forall p1 p2 p3:pt. turns_left p1 p2 p3 -> turns_left p2 p3 p1
  19
  20 lemma knuth2:
  21   forall p1 p2 p3:pt. turns_left p1 p2 p3 -> not (turns_left p2 p1 p3)
  22
  23 lemma knuth4:
  24   forall p1 p2 p3 p4:pt.
  25     turns_left p1 p2 p4 ->
  26     turns_left p2 p3 p4 ->
  27     turns_left p3 p1 p4 ->
  28     turns_left p1 p2 p3
  29
  30 lemma knuth5:
  31   forall p1 p2 p3 p4 p5:pt.
  32     turns_left p1 p2 p3 ->
  33     turns_left p1 p2 p4 ->
  34     turns_left p1 p2 p5 ->
  35     turns_left p3 p2 p4 ->
  36     turns_left p4 p2 p5 ->
  37     turns_left p3 p2 p5
  38
  39 lemma knuth3:
  40   forall p1 p2 p3:pt.
  41    not (colinear p1 p2 p3) -> turns_left p1 p2 p3 \/ turns_left p1 p3 p2
  42
  43
  44 (* test: p1 = (0,0), p2 = (1,0) p3 = (1,1) *)
  45
  46 goal test1 :
  47    turns_left {x=0.0;y=0.0} {x=1.0;y=0.0} {x=1.0;y=1.0}
  48
  49
  50 use list.List
  51 use list.Length
  52 use list.NthNoOpt
  53
  54 type path = list pt
  55
  56 predicate is_inside_convex_hull (p:pt) (pa:path) =
  57   let l = length pa in
  58   forall i:int. 0 <= i < l ->
  59     let i' = if i = l-1 then 0 else i+1 in
  60     let p0 = nth i pa in
  61     let p1 = nth i' pa in
  62     turns_left p0 p1 p
  63
  64 predicate is_ccw_convex (pa:path) =
  65   let l = length pa in
  66   forall i:int. 0 <= i < l ->
  67     let i' = if i = l-1 then 0 else i+1 in
  68     let p0 = nth i pa in
  69     let p1 = nth i' pa in
  70     forall j:int. 0 <= j < l /\ j <> i /\ j <> i' ->
  71        turns_left p0 p1 (nth j pa)
  72
  73 end
  74
  75
  76 module Incremental
  77
  78
  79
  80 end
  81
  82
  83 module Jarvis
  84
  85 (*
  86
  87 on calcule le point p d'ordonnee minimale, et d'abscisse minimale
  88 parmi ceux-ci
  89
  90 on recherche le point q minimum pour la relation R x y = turns_left p x y
  91
  92 p0 = p
  93 p1 = q
  94 si q = p0: fini
  95   sinon
  96     recommencer avec p := q
  97
  98 *)
  99
 100 use ConvexSet
 101 clone set.SetApp with type elt = pt
 102 use ref.Ref
 103
 104 val predicate lower (p q:pt)
 105   ensures { result <-> p.y <. q.y \/ (p.y = q.y /\ p.x <=. q.x) }
 106
 107 let lowest_pt (s:set) : (pt, set)
 108   requires { not (is_empty s) }
 109   ensures { let (p,r) = result in
 110     s = add p r /\ forall q:pt. mem q r -> lower p q }
 111 =
 112   let p = ref (choose s) in
 113   let r = ref (remove !p s) in
 114   while not (is_empty !r) do
 115      variant { cardinal !r }
 116      invariant { mem !p s }
 117      invariant { subset !r s }
 118      invariant { forall q:pt. mem q s /\ not (mem q !r) -> lower !p q }
 119      let q = choose !r in
 120      if lower q !p then p := q;
 121      r := remove q !r
 122   done;
 123   (!p,remove !p s)
 124
 125 (*
 126 let rightest_pt (p:pt) (s:set pt) : (pt, set pt)
 127   requires { not (is_empty s) }
 128   ensures { let (p,r) = result in
 129     s = add p r /\ forall q:pt. mem q r -> lower p q }
 130 =
 131   let p = ref (choose s) in
 132   let r = ref (remove !p s) in
 133   while not (is_empty !r) do
 134      invariant { mem !p s }
 135      invariant { subset !r s }
 136      invariant { forall q:pt. mem q s /\ not (mem q !r) -> lower !p q }
 137      let q = choose !r in
 138      if lower q !p then p := q;
 139      r := remove q !r
 140   done;
 141   (!p,remove !p s)
 142
 143
 144 let jarvis (s:set pt) : list pt =
 145    let p0 = lowest_pt s in
 146 *)
 147
 148
 149 end
 150
 151
 152
 153 module Graham
 154
 155
 156 use ConvexSet
 157
 158 (*
 159 let convex_hull (l:path) : path =
 160     (* all pts of the result are points of the input *)
 161   ensures { forall p:pt. List.mem p result -> List.mem p l }
 162     (* the output forms a ccw convex circuit *)
 163   ensures { is_ccw_convex result }
 164     (* all pt of the input are inside the convex hull of the output *)
 165   ensures { forall p:pt. List.mem p l /\ not List.mem p result ->
 166     is_inside_convex_hull p result
 167   }
 168   let min,rem = find_minimum_pt l in
 169   let sorted = sorted_increasing_angle min rem in
 170   match sorted with
 171     | Nil -> l
 172     | Cons p r ->
 173        let done = ref (Cons p min) in
 174        let todo = ref r in
 175        try while true do
 176          match !todo with
 177          | Nil ->
 178        with Exit -> !done
 179
 180 *)
 181
 182 end