Formal Verification of a Geometry Algorithm Yves Bertot October - PowerPoint PPT Presentation

Formal Verification of a Geometry Algorithm Yves Bertot October 2018 1 / 38

Outline ◮ Triangulation algorithm ◮ Abstract presentation and successive refinements ◮ Symmetries of the triangle and the convex hull 2 / 38

Triangulations and formal proofs ◮ Bigger context: robot motion ◮ Vorono¨ ı diagrams and Delaunay triangulation ◮ Dufourd&Bertot10 starts from an existing triangulation ◮ Dufourd&Bertot10 uses a specific data-structure ◮ darts: edge extremities ◮ two darts to make an edge, α 0 ◮ darts around a point are connected, α 1 ◮ A third permutation ( α 0 ◦ α 1 ) − 1 enumerates facets ◮ Astonishingly handy for implementation in C-like languages ◮ α 0 and α 1 are just pointers 3 / 38

Drawing with darts α 1 d 2 d 1 α 0 d 4 d 0 d 3 d 5 ◮ d 0 , d 1 , d 2 and d 3 , d 4 , d 5 form triangles ◮ unnamed darts form the outer boundary 4 / 38

A more abstract approach to triangulations ◮ Fix a finite set of points ◮ Each point has a pair of coordinates ◮ all operations on coordinates are deemed exact ◮ Only ring operations are necessary ◮ Describe a triangulation as a set of triangles ◮ Describe each triangle as a three-point set ◮ The formalized set theory we use has all basic set operations ◮ set difference, singletons, cardinal, union, intersection ◮ all sets are naturally finite and all operations are decidable 5 / 38

Crucial geometry property ◮ Given a segment [ a , b ], does c lie to the left or to the right? ◮ Knuth (1992) proposes to abstract this question as a single predicate ccw (counter-clockwise) ◮ 5 basic properties (called Knuth axioms ) 1. ccw ( a , b , c ) = ccw ( b , c , a ) 2. ccw ( a , b , c ) ⇒ ¬ ccw ( b , a , c ) 3. card { a , b , c } = 3 ⇒ ccw ( a , b , c ) ∨ ccw ( b , a , c ) 4. ccw ( a , b , d ) ∧ ccw ( b , c , d ) ∧ ccw ( c , a , d ) ⇒ ccw ( a , b , c ) 5. ccw ( a , b , c ) ∧ ccw ( a , b , d ) ∧ ccw ( a , b , e ) ∧ ccw ( a , c , d ) ∧ ccw ( a , d , e ) ⇒ ccw ( a , c , e ) ◮ Axiom 3 is a property of the data: no three points aligned 6 / 38

Illustration of Axiom 4 c d b a 7 / 38

Illustration of Axiom 5 c d b e a 8 / 38

Implementing and proving ccw ◮ To every triangle we can associate an orientedsurface ◮ ccw just says that this oriented surface is positive ◮ the oriented surface is computed by a determinant  � �  1 x a y a � � � � ◮ ccw ( a , b , c ) =  0 < 1 x b y b � �  � � 1 x c y c � � 9 / 38

Discussion on numerical issues ◮ The assumption of exact numerical computations is reasonable ◮ Operations for computing the counter-clockwise operation rely on polynomials with degree 2 ◮ Limiting the precision for the inputs and using double precision numbers is enough ◮ alternatively using 32-bit integers for inputs and 64 bit integers for results ◮ In Coq itself, all computations are exact by default ◮ We have to use a discrete number ring to be algorithmically relevant 10 / 38

Derived predicates ◮ Whether point d is inside a triangle abc can be described using ccw ◮ The convex hull of a set A of points can be described using ccw ◮ a sequence s 0 , . . . , s n such that, for every i and x different from s i , s i +1 , ccw ( s i , s i +1 , x ) and similarly for ccw ( s n , s 0 ) ◮ it can also be described by giving s 0 and the function f that maps s i to s i +1 and s n to s 0 11 / 38

Convex hull illustration s 6 s 5 f s 0 s 4 f s 1 s 3 f s 2 12 / 38

A naive triangulation algorithm ◮ Add points one by one ◮ When you have three points make the first triangle ◮ When the new point is in an existing triangle ◮ remove that one ◮ add 3 new triangles made with the edges of the old triangle and the new point ◮ Otherwise, find the boundary edges that are separating ◮ For each such separating boundary edge, make a new triangle by adding the new point 13 / 38

Naive algorithm illustration 14 / 38

Naive algorithm illustration 15 / 38

Naive algorithm illustration 16 / 38

Naive algorithm illustration 17 / 38

Mathematical formulation ◮ Blue step: find t ∈ T such that p is inside t ◮ T \{ t } ∪ { t \ q ∪ { p } | q ∈ t } ◮ Red step: ◮ T ∪ � �  { p } ∪ ( t \ q ) t ∈ T q ∈ t ,    boundary ( t \ q )  separated ( t \ q , p , q )   18 / 38

Proving properties of this algorithm ◮ First main specification ◮ All elements of the triangulation have 3 points ◮ The union of all triangles is the input set ◮ The whole convex hull is covered by triangles ◮ Need Geometry proofs ◮ Re-use Knuth’s approach to the convex Hull 19 / 38

Examples of proofs ◮ A point inside a triangle left of a segment is also left of the segment Lemma inside_triangle_left_of_edge a b c d e f : general_position ([set a; b; c; d; e; f]) -> ccw a b c -> ccw a b d -> ccw a b e -> ccw c d f -> ccw c f e -> ccw f d e -> ccw a b f. ◮ Many cases, but also many symmetries 20 / 38

Illustration c f e d b a 21 / 38

Symmetries in this figure ◮ The statement does not say which of c , d , or e is the the leftmost point with respect to a ◮ If we assume c to be the leftmost, it does not say whether ( ccw a d e ) holds or not. ◮ There are 6 cases in the proof, each case require 3 uses of Knuth’s 5th axiom 22 / 38

Illustration c e f d b a 23 / 38

Key points concerning the convex hull ◮ In the red case, red edges are grouped together ◮ There are only two points adjacent to a blue and a red edge ◮ This requires a proof by induction (80 lines of dense Coq script) Lemma partition_red_blue a: a \in cover s’ -> ccw (f a) a p -> exists b, exists i, b \in cover s’ /\ (0 < i < size (orbit f b))%N /\ (forall j, (j < i)%N -> ccw (iter j.+1 f b) (iter j f b) p) /\ (forall j, (i <= j < size (orbit f b))%N -> ccw (iter j f b) (iter j.+1 f b) p). 24 / 38

A critique of the mathematical algorithm ◮ Complexity is terrible ◮ Generic operations on sets of sets ◮ Refinements ◮ Caching data ◮ Sparse representation for sets 25 / 38

Refinements 1. Return the convex hull with the triangulation 2. Take as input a list of points (instead of a subset) 3. Output a list of triangles (instead of a subset) 26 / 38

Convex Hull refinement ◮ Refine the output type from { set P } to { set P } * (P → P) * P ◮ The first component is the same as for the previous algorithm ◮ The second component is a function circulating the convex hull ◮ The third component is a point on the convex hull ◮ The convex hull does not change in the blue case ◮ In the red case some computation must be done ◮ Detect and suppress red edges, keep blue edges ◮ Add two new edges from and two purple points 27 / 38

Convex hull computation in the red case 28 / 38

Defining the new convex hull in the red case ◮ I chose to use the datatype of “finite functions” ◮ Just define the new convex hull as a variation over the old one ◮ Proofs about this new path ◮ values for all blue points are the same as the old path function ◮ Also true for the last purple point ◮ The new path constitutes the new convex hull Definition new_path (f : { ffun P -> P } ) (a b np : P) := finfun [fun x => f x with a |-> np, np |-> b]. 29 / 38

Proofs about the new convex hull ◮ The new path is a cycle ◮ The orbit of the new point describes the new boundary edges ◮ 400 lines of proof ◮ The new path has the convex hull property ◮ 120 lines of proof 30 / 38

Correctness statement for the first refinement Lemma naive_naive2 n (s : { set P } ) : (3 <= #|s| <= n)%N -> general_position s -> exists chf b, naive2 n s = (naive n s, chf, b) /\ convex_hull_path coords s (boundary_edge (naive n s)) chf /\ b \in cover (boundary_edge (naive n s)). 31 / 38

Input refinement ◮ Input points from a list instead of a set ◮ Move from a mathematical concept to a computer science one ◮ Simple recursion scheme for lists (none for sets) ◮ Need to replace a non-deterministic choice with a deterministic one ◮ All proofs about the initial algorithm need to be abstract wrt the pick function ◮ Lists that are equal up to permutation may lead to different triangulations ◮ Introduce a function pick seq set that chooses elements in a set according to the order given by a list, when applicable 32 / 38

Code for the second refinement Definition pick_seq_set (s : seq P) (s’ : { set P } ) : option P := if (find (mem s’) s < size s)%N then Some (nth p0 s (find (mem s’) s)) else pick_set s’. Definition naive2’ s n s’ := triangulation_algorithm.naive2 coords p0 (pick_seq_set s) pick_triangle n s’. Lemma naive2_naive3 (s : seq P) : (3 <= size s)%N -> uniq s -> naive3 s = naive2’ s (size s) [set y in s]. ◮ This proof is short: 25 lines 33 / 38

abstractions in the refinement ladder ◮ naive and naive2 are shown to return the same triangulation when using the same pick set function ◮ naive2 and naive3 are shown to return the same triangulation when naive2 is using pick seq set 34 / 38