projects_discrete_geom/bezdek_reid/bezdek_reid.hl

   1 (*
   2
   3 Thomas C. Hales
   4 January 12, 2013.
   5
   6 Notes on the paper "Contact Graphs of Unit Sphere Packings Revisited"
   7 by Bezdek and Reid.
   8 preprint date October 17, 2012.
   9
  10
  11 We show the following results:
  12
  13 There are no kissing arrangements of 12 spheres with at least 25 contact edges.
  14 There are no kissing arrangements of 12 spheres with at least 11 contact triples.
  15
  16 ****
  17 This is not a stand-alone piece of code.
  18 For the graph generation, it requires the installation of java.
  19
  20 It also requires parts of the Flyspeck project, especially
  21   Graph_control for graph generation, and Glpk_link to process the raw data,
  22
  23 We also use some functions such as (--) and chop_list from HOL-Light's lib.ml.
  24
  25 As the code is currently written, it assumes that the entire Flyspeck project
  26 has been loaded, but with a bit of work, it could be made to depend on just
  27 the indicated modules.
  28
  29 ****
  30
  31 The relevant output is coord_edge and coord_triple, which exhibit inconsistent
  32 coordinate assignments for the nodes of the graphs.
  33
  34 *)
  35
  36 flyspeck_needs ("../graph_generator/graph_control.hl");;
  37 flyspeck_needs "../glpk/glpk_link.ml";;
  38
  39
  40 module Bezdek_reid = struct
  41
  42   open Graph_control;;
  43   open Glpk_link;;
  44   open List;;
  45
  46   (* The file names need to be adjusted for different users *)
  47   let modelfile = "/Users/thomashales/Desktop/googlecode/flyspeck/projects_discrete_geom/bezdek_reid/br_model.mod";;
  48
  49   (* This is the default output location for the graph generation, as set in the java code. *)
  50   let tmpfile = "/tmp/graph_out.txt";;
  51
  52 let bezdek_reid_properties_edge =
  53   {
  54     properties_id =
  55       "K. Bezdek and S. Reid: Contact Graphs of Unit Sphere Packings Revisited (2012)";
  56     ignore_archive=true;
  57     exclude_degree2=false;
  58     exclude_pent_qrtet=true;
  59     exclude_2_in_quad=true;
  60     exclude_1_in_tri=true;
  61
  62     (* require exactly 12 spheres *)
  63     vertex_count_min=12;
  64     vertex_count_max=12;
  65
  66     (* degree at most 5 at every node *)
  67     node_card_max=5;
  68     node_card_max_at_exceptional_vertex=5;
  69
  70     (* if the are more than 5.5 edges removed, then ignore *)
  71     squander_target=55;
  72     score_target= -1;
  73
  74     (* a quad face = 1 edge removed, pent = 2 edges removed, etc. *)
  75     table_weight_d = [0;0;0;0;10;20;30;40;50];
  76     table_weight_a = [(0,0);(1,0);(2,0);(3,0);(4,0)];
  77     table_weight_b = [(0,3,30);(0,4,40);(0,5,50);
  78                       (1,3,30);(1,4,40);
  79                       (2,2,20);(2,3,30);
  80                       (3,2,20);
  81                       (4,1,10)];
  82   };;
  83
  84 let bezdek_reid_properties_triplet =
  85   {
  86     properties_id =
  87       "K. Bezdek and S. Reid: Contact Graphs of Unit Sphere Packings Revisited (2012)";
  88     ignore_archive=true;
  89     exclude_degree2=false;
  90     exclude_pent_qrtet=true;
  91     exclude_2_in_quad=true;
  92     exclude_1_in_tri=true;
  93
  94     (* require exactly 12 spheres *)
  95     vertex_count_min=12;
  96     vertex_count_max=12;
  97
  98     (* degree at most 5 at every node *)
  99     node_card_max=5;
 100     node_card_max_at_exceptional_vertex=5;
 101
 102     (* if the are more than 9.5 non-contact edges, then ignore *)
 103     squander_target=95;
 104     score_target= -1;
 105
 106     (* a quad face = 1 edge removed, pent = 2 edges removed, etc. *)
 107     (* allow up to an 11-gon. 11-gon, 10-gon, 9-gon ruled out
 108        by extra n-gons at remaining vertices.  *)
 109     table_weight_d = [0;0;0;0;20;30;40;50;60];
 110     table_weight_a = [(0,0);(1,0);(2,0);(3,0);(4,0)];
 111     table_weight_b = [(0,3,60);(0,4,80);(0,5,100);
 112                       (1,3,60);(1,4,80);
 113                       (2,2,40);(2,3,60);
 114                       (3,2,40);
 115                       (4,1,20)];
 116   };;
 117
 118
 119 let max_degree hypl =
 120   let fl = List.flatten (map (fun h->  ((mk_bb h).std_faces)) hypl) in
 121   let ll = map List.length fl in
 122     end_itlist max ll;;
 123
 124 (* initial triangle *)
 125 let sqrt = Pervasives.sqrt;;
 126 let p1 = (2.0,0.0,0.0);;
 127 let p2 = (1.0,sqrt(3.0),0.0);;
 128 let p3 = (1.0,1.0/. sqrt(3.0), sqrt(8.0) /. sqrt(3.0));;
 129
 130
 131 (* We include a few vector functions from tikz.ml *)
 132
 133 let (+...) (x1,x2,x3) (y1,y2,y3) = (x1 +. y1, x2 +. y2, x3+. y3);;
 134
 135 let (-...) (x1,x2,x3) (y1,y2,y3) = (x1 -. y1, x2 -. y2, x3-. y3);;
 136
 137 let ( %... ) s (x1,x2,x3) = (s *. x1, s *. x2, s*. x3);;
 138
 139 let ( *... ) (x1,x2,x3) (y1,y2,y3) = (x1 *. y1 +. x2 *. y2 +. x3 *. y3);;
 140
 141 let cross (x1,x2,x3) (y1,y2,y3) =
 142   (x2 *. y3 -. x3 *. y2, x3 *. y1 -. x1 *. y3, x1 *. y2 -. x2 *. y1);;
 143
 144 let normalize3 x = (1.0 /. sqrt(x *... x)) %... x;;
 145
 146 let dist3 x y =
 147   let z = x -... y in sqrt (z *... z);;
 148
 149 let reflect a b c =
 150   let u = normalize3 (cross b c) in
 151   a -... ((2.0 *. (u *... a)) %... u);;
 152
 153 let list_of (a,b,c) = [a;b;c];;
 154
 155 let tuple_of [a;b;c] = (a,b,c);;
 156
 157 let common a b  = intersect (list_of a) (list_of b);;
 158
 159 let adj3 a b = (length (list_of a) = 3) && (length (list_of b) = 3) &&
 160   length (common a b) = 2;;
 161
 162 let rec outer x y =
 163   match x with
 164     | [] -> []
 165     | a::r -> (map (fun i -> (a,i)) y) @ (outer r y);;
 166
 167 let rec find p l =
 168   match l with
 169       [] -> failwith "find"
 170     | (h::t) -> if p(h) then h else find p t;;
 171
 172 let find_adj_pair al bl =
 173   let (s,t)  =   find (fun (s,t) -> adj3 s t) (outer al bl) in
 174   let c = common s t in
 175   let (c1,c2) = (List.nth c 0,List.nth c 1) in
 176   let (sl,tl) = (list_of s, list_of t) in
 177   let a = hd (subtract sl c) in
 178   let b = hd (subtract tl c) in
 179     (a,c1,c2,b,s,t);;
 180
 181 let update_coord (cl,br1 ,br2) =
 182   let (a,b,c,a',s,t) = find_adj_pair br1 br2 in
 183   let pa = assoc a cl in
 184   let pb = assoc b cl in
 185   let pc = assoc c cl in
 186   let pa' = reflect pa pb pc in
 187     ((a',pa')::cl , t::br1 , subtract br2 [t]);;
 188
 189 let rec mk_coords (cl,br1,br2) =
 190   if (br2 = [] or not( can (update_coord) (cl,br1,br2)))
 191   then (cl,br1,br2) else mk_coords (update_coord (cl,br1,br2));;
 192
 193 let triangles_in_archive_string ll i =
 194   let case1 = List.nth ll i in
 195   let br = (mk_bb (case1)).std_faces in
 196   let br_tri =     filter (fun t -> List.length t = 3) br in
 197     map tuple_of br_tri;;
 198
 199 let do_one_case brt =
 200   let (br1,br2) = chop_list 1 (brt) in
 201   let [(a1,a2,a3)] = br1 in
 202   let coordlist =     [(a1,p1);(a2,p2);(a3,p3)] in
 203       mk_coords (coordlist,br1,br2);;
 204
 205 let rec mk_all_coords (cll,br2) =
 206   if br2 = [] then (cll,[])
 207   else
 208     let (cl,_,br2') = do_one_case  br2 in
 209       mk_all_coords (cl::cll, br2');;
 210
 211 (* consistency checks on coordinates *)
 212
 213 let is_packing dc =
 214   let vl = setify (map fst dc) in
 215   let dis_ok = map (fun (i,j) -> if (j<=i) then true
 216                     else (1.99 < dist3 (assoc i dc) (assoc j dc)))
 217     (outer vl vl) in
 218     (setify dis_ok = [true]);;
 219
 220 let self_consis_one dc i =
 221   let eps = 0.001 in
 222   let dci = filter (fun (j,_) -> (j =i)) dc in
 223   let dr = map snd dci in
 224   let dis_ok = map (fun (a,b) -> dist3 a b < eps) (outer dr dr) in
 225     setify dis_ok = [true];;
 226
 227 let self_consis dc =
 228   let vl = setify (map fst dc) in
 229     forall (self_consis_one dc) vl;;
 230
 231 let inter_consis_one dc1 dc2 =
 232   let vl = intersect (map fst dc1) (map fst dc2) in
 233   let a1 i = assoc i dc1 in
 234   let a2 i = assoc i dc2 in
 235   let abs = Pervasives.abs_float in
 236   let eps = 0.01 in
 237     forall (fun (i,j) -> abs(dist3 (a1 i) (a1 j) -. dist3 (a2 i) (a2 j)) < eps)
 238       (outer vl vl);;
 239
 240 let inter_consis cll =
 241   forall (fun (a,b) -> inter_consis_one a b) (outer cll cll);;
 242
 243 let mk_c ll n =
 244   let br = triangles_in_archive_string ll n in
 245   let (cll,_) = mk_all_coords ([],br) in
 246     cll;;
 247
 248 let test_coords ll n =
 249   let cll = mk_c ll n in
 250     (forall (is_packing) cll) &&
 251       (forall (self_consis) cll) &&
 252       (inter_consis cll);;
 253
 254 (* we see by inspection from the explicit coordinates that they are
 255    overdetermined. test_edge and test_triplet should both evaluate to [false] *)
 256
 257 (* edge case processing *)
 258
 259 let hypermap_edge =
 260   let archiveraw_edge = tmpfile  in
 261   let _ = Graph_control.run bezdek_reid_properties_edge in
 262     Glpk_link.strip_archive archiveraw_edge;;
 263
 264 let count_edge = List.length hypermap_edge;; (* 18 cases *)
 265 let max_degree_edge = max_degree hypermap_edge;;
 266 let range_edge =  (0--(List.length hypermap_edge - 1));;
 267 let test_edge = setify (map (test_coords hypermap_edge) range_edge);;
 268
 269
 270 (* triplet case *)
 271
 272 let hypermap_triplet =
 273   let archiveraw_triplet = tmpfile in
 274   let _ = Graph_control.run bezdek_reid_properties_triplet in
 275     Glpk_link.strip_archive archiveraw_triplet;;
 276 let count_triplet = List.length hypermap_triplet;; (* 1335 cases *)
 277 let max_degree_edge = max_degree hypermap_triplet;;  (* 8 *)
 278 let range_triplet = (0--(List.length hypermap_triplet - 1));;
 279 let test_triplet = setify (map (test_coords hypermap_triplet) range_triplet);;
 280
 281
 282
 283 end;;