(*

Thomas C. Hales
January 12, 2013.

Notes on the paper "Contact Graphs of Unit Sphere Packings Revisited"
by Bezdek and Reid.
preprint date October 17, 2012.


We show the following results:

There are no kissing arrangements of 12 spheres with at least 25 contact edges. 
There are no kissing arrangements of 12 spheres with at least 11 contact triples. 

****
This is not a stand-alone piece of code.
For the graph generation, it requires the installation of java.

It also requires parts of the Flyspeck project, especially
  Graph_control for graph generation, and Glpk_link to process the raw data,

We also use some functions such as (--) and chop_list from HOL-Light's lib.ml.

As the code is currently written, it assumes that the entire Flyspeck project
has been loaded, but with a bit of work, it could be made to depend on just
the indicated modules.

****

The relevant output is coord_edge and coord_triple, which exhibit inconsistent
coordinate assignments for the nodes of the graphs.

*)

flyspeck_needs ("../graph_generator/graph_control.hl");;
flyspeck_needs "../glpk/glpk_link.ml";;


module Bezdek_reid = struct

  open Graph_control;;
  open Glpk_link;;
  open List;;

  (* The file names need to be adjusted for different users *)
  let modelfile = "/Users/thomashales/Desktop/googlecode/flyspeck/projects_discrete_geom/bezdek_reid/br_model.mod";;

  (* This is the default output location for the graph generation, as set in the java code. *)
  let tmpfile = "/tmp/graph_out.txt";;

let bezdek_reid_properties_edge = 
  {
    properties_id = 
      "K. Bezdek and S. Reid: Contact Graphs of Unit Sphere Packings Revisited (2012)";
    ignore_archive=true;
    exclude_degree2=false;
    exclude_pent_qrtet=true; 
    exclude_2_in_quad=true;
    exclude_1_in_tri=true;

    (* require exactly 12 spheres *)
    vertex_count_min=12;
    vertex_count_max=12;

    (* degree at most 5 at every node *)
    node_card_max=5;
    node_card_max_at_exceptional_vertex=5;

    (* if the are more than 5.5 edges removed, then ignore *)
    squander_target=55;
    score_target= -1;

    (* a quad face = 1 edge removed, pent = 2 edges removed, etc. *)
    table_weight_d = [0;0;0;0;10;20;30;40;50];
    table_weight_a = [(0,0);(1,0);(2,0);(3,0);(4,0)];
    table_weight_b = [(0,3,30);(0,4,40);(0,5,50);
		      (1,3,30);(1,4,40);
		      (2,2,20);(2,3,30);
		      (3,2,20);
		      (4,1,10)];
  };;

let bezdek_reid_properties_triplet = 
  {
    properties_id = 
      "K. Bezdek and S. Reid: Contact Graphs of Unit Sphere Packings Revisited (2012)";
    ignore_archive=true;
    exclude_degree2=false;
    exclude_pent_qrtet=true; 
    exclude_2_in_quad=true;
    exclude_1_in_tri=true;

    (* require exactly 12 spheres *)
    vertex_count_min=12;
    vertex_count_max=12;

    (* degree at most 5 at every node *)
    node_card_max=5;
    node_card_max_at_exceptional_vertex=5;

    (* if the are more than 9.5 non-contact edges, then ignore *)
    squander_target=95;
    score_target= -1;

    (* a quad face = 1 edge removed, pent = 2 edges removed, etc. *)
    (* allow up to an 11-gon. 11-gon, 10-gon, 9-gon ruled out 
       by extra n-gons at remaining vertices.  *)
    table_weight_d = [0;0;0;0;20;30;40;50;60];
    table_weight_a = [(0,0);(1,0);(2,0);(3,0);(4,0)];
    table_weight_b = [(0,3,60);(0,4,80);(0,5,100);
		      (1,3,60);(1,4,80);
		      (2,2,40);(2,3,60);
		      (3,2,40);
		      (4,1,20)];
  };;


let max_degree hypl = 
  let fl = List.flatten (map (fun h->  ((mk_bb h).std_faces)) hypl) in
  let ll = map List.length fl in
    end_itlist max ll;;

(* initial triangle *)
let sqrt = Pervasives.sqrt;;
let p1 = (2.0,0.0,0.0);;
let p2 = (1.0,sqrt(3.0),0.0);;
let p3 = (1.0,1.0/. sqrt(3.0), sqrt(8.0) /. sqrt(3.0));;


(* We include a few vector functions from tikz.ml *)

let (+...) (x1,x2,x3) (y1,y2,y3) = (x1 +. y1, x2 +. y2, x3+. y3);;

let (-...) (x1,x2,x3) (y1,y2,y3) = (x1 -. y1, x2 -. y2, x3-. y3);;

let ( %... ) s (x1,x2,x3) = (s *. x1, s *. x2, s*. x3);; 

let ( *... ) (x1,x2,x3) (y1,y2,y3) = (x1 *. y1 +. x2 *. y2 +. x3 *. y3);;

let cross (x1,x2,x3) (y1,y2,y3) = 
  (x2 *. y3 -. x3 *. y2, x3 *. y1 -. x1 *. y3, x1 *. y2 -. x2 *. y1);;

let normalize3 x = (1.0 /. sqrt(x *... x)) %... x;;

let dist3 x y = 
  let z = x -... y in sqrt (z *... z);;

let reflect a b c = 
  let u = normalize3 (cross b c) in
  a -... ((2.0 *. (u *... a)) %... u);;

let list_of (a,b,c) = [a;b;c];;

let tuple_of [a;b;c] = (a,b,c);;

let common a b  = intersect (list_of a) (list_of b);;

let adj3 a b = (length (list_of a) = 3) && (length (list_of b) = 3) && 
  length (common a b) = 2;;

let rec outer x y = 
  match x with
    | [] -> []
    | a::r -> (map (fun i -> (a,i)) y) @ (outer r y);;

let rec find p l =
  match l with
      [] -> failwith "find"
    | (h::t) -> if p(h) then h else find p t;;

let find_adj_pair al bl = 
  let (s,t)  =   find (fun (s,t) -> adj3 s t) (outer al bl) in
  let c = common s t in
  let (c1,c2) = (List.nth c 0,List.nth c 1) in
  let (sl,tl) = (list_of s, list_of t) in
  let a = hd (subtract sl c) in
  let b = hd (subtract tl c) in
    (a,c1,c2,b,s,t);;

let update_coord (cl,br1 ,br2) = 
  let (a,b,c,a',s,t) = find_adj_pair br1 br2 in
  let pa = assoc a cl in
  let pb = assoc b cl in
  let pc = assoc c cl in
  let pa' = reflect pa pb pc in
    ((a',pa')::cl , t::br1 , subtract br2 [t]);;

let rec mk_coords (cl,br1,br2) = 
  if (br2 = [] or not( can (update_coord) (cl,br1,br2)))
  then (cl,br1,br2) else mk_coords (update_coord (cl,br1,br2));;

let triangles_in_archive_string ll i = 
  let case1 = List.nth ll i in
  let br = (mk_bb (case1)).std_faces in
  let br_tri =     filter (fun t -> List.length t = 3) br in
    map tuple_of br_tri;;

let do_one_case brt = 
  let (br1,br2) = chop_list 1 (brt) in
  let [(a1,a2,a3)] = br1 in
  let coordlist =     [(a1,p1);(a2,p2);(a3,p3)] in
      mk_coords (coordlist,br1,br2);;

let rec mk_all_coords (cll,br2) = 
  if br2 = [] then (cll,[])
  else
    let (cl,_,br2') = do_one_case  br2 in
      mk_all_coords (cl::cll, br2');;

(* consistency checks on coordinates *)

let is_packing dc = 
  let vl = setify (map fst dc) in
  let dis_ok = map (fun (i,j) -> if (j<=i) then true 
		    else (1.99 < dist3 (assoc i dc) (assoc j dc)))
    (outer vl vl) in
    (setify dis_ok = [true]);;

let self_consis_one dc i = 
  let eps = 0.001 in
  let dci = filter (fun (j,_) -> (j =i)) dc in
  let dr = map snd dci in
  let dis_ok = map (fun (a,b) -> dist3 a b < eps) (outer dr dr) in
    setify dis_ok = [true];;

let self_consis dc = 
  let vl = setify (map fst dc) in
    forall (self_consis_one dc) vl;;

let inter_consis_one dc1 dc2 = 
  let vl = intersect (map fst dc1) (map fst dc2) in
  let a1 i = assoc i dc1 in
  let a2 i = assoc i dc2 in
  let abs = Pervasives.abs_float in
  let eps = 0.01 in
    forall (fun (i,j) -> abs(dist3 (a1 i) (a1 j) -. dist3 (a2 i) (a2 j)) < eps)
      (outer vl vl);;
  
let inter_consis cll = 
  forall (fun (a,b) -> inter_consis_one a b) (outer cll cll);;

let mk_c ll n = 
  let br = triangles_in_archive_string ll n in
  let (cll,_) = mk_all_coords ([],br) in
    cll;;

let test_coords ll n = 
  let cll = mk_c ll n in 
    (forall (is_packing) cll) &&
      (forall (self_consis) cll) &&
      (inter_consis cll);;

(* we see by inspection from the explicit coordinates that they are
   overdetermined. test_edge and test_triplet should both evaluate to [false] *)

(* edge case processing *)

let hypermap_edge = 
  let archiveraw_edge = tmpfile  in
  let _ = Graph_control.run bezdek_reid_properties_edge in
    Glpk_link.strip_archive archiveraw_edge;;

let count_edge = List.length hypermap_edge;; (* 18 cases *)
let max_degree_edge = max_degree hypermap_edge;;
let range_edge =  (0--(List.length hypermap_edge - 1));;
let test_edge = setify (map (test_coords hypermap_edge) range_edge);;


(* triplet case *)

let hypermap_triplet = 
  let archiveraw_triplet = tmpfile in
  let _ = Graph_control.run bezdek_reid_properties_triplet in
    Glpk_link.strip_archive archiveraw_triplet;;
let count_triplet = List.length hypermap_triplet;; (* 1335 cases *)
let max_degree_edge = max_degree hypermap_triplet;;  (* 8 *)
let range_triplet = (0--(List.length hypermap_triplet - 1));;
let test_triplet = setify (map (test_coords hypermap_triplet) range_triplet);;



end;;