(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(*                                                                            *)
(* Chapter: Packing                                                            *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2010-02-21                                                           *)
(* ========================================================================== *)

(* definitions needed for the chapter Packing *)


flyspeck_needs "volume/vol1.hl";;

module Pack_defs = struct

  let sol = Vol1.sol;;

   let negligible_fun_p = new_definition `negligible_fun_p f S (p:real^N) = (? (C:real). (&0<= C)/\ ! (r:real). ( &1 <= r) ==> (sum (S INTER ball(p,r)) f) <= C * (r pow 2))`;;

   let negligible_fun_0 = new_definition `negligible_fun_0 f S  = negligible_fun_p f S ((vec:num->real^3) 0)`;;

   let fcc_compatible= new_definition `fcc_compatible f (S:real^N->bool) = (!v. v IN S ==> sqrt( &32) <= measure(voronoi_open S v) + f v)`;;

   let kepler_conjecture = new_definition `kepler_conjecture = (!V. packing V /\ saturated V ==> (?(c:real). !(r:real). (&1<= r ==> vol ((UNIONS {ball(v:real^3, &1) | v IN V}) INTER ball(vec 0 ,r))/ vol (ball(vec 0, r))<= pi/ sqrt( &18) + c/ r)))`;;

let HL = new_definition `hl (ul:(real^N)list) = radV (set_of_list ul)`;;

let REVERSE_TABLE  = define `(REVERSE_TABLE (f:num->A) 0 = []) /\ 
   (REVERSE_TABLE f (SUC i) = CONS (f i)  ( REVERSE_TABLE f i))`;;

let TABLE = new_definition `!(f:num->A) k. TABLE f k = REVERSE (REVERSE_TABLE f k)`;;

let left_action = new_definition `!p f x. left_action (p:A->C) (f:A->B) x = f(inverse p x)`;;

let left_action_list = new_definition `left_action_list p (ul:(A)list) = TABLE (\i. EL (inverse p i) ul) (LENGTH ul)`;;

let DROP = define `(DROP (ul:(A)list) 0  = TL ul) /\ (DROP ul (SUC i) = CONS(HD ul) (DROP (TL ul) i))`;;

let mxi = new_definition `!V ul. mxi V ul = 
  if (hl (truncate_simplex 2 ul) >= sqrt(&2)) 
  then omega_list_n V ul 2 
  else (@p. ((p IN convex hull 
    { omega_list_n V ul 2  , omega_list_n V ul 3}) /\
    (dist(p, HD ul) = sqrt(&2))))`;;

let mcell0 = 
  new_definition `!V ul. mcell0 V ul = rogers V ul DIFF ball(HD ul,sqrt(&2))`;;

(* June 17, 2011, changed to mcell1, mcell2 to
   weak ineq sqrt(&2) <= hl ul to fix an error
   reported by Ky *)

let mcell1 = new_definition `!V ul. mcell1 V ul = 
   if (sqrt(&2) <= hl ul) 
   then
     rogers V ul INTER cball(HD ul, (sqrt (&2))) DIFF 
     rcone_gt (HD ul) (HD (TL ul))
        (hl (truncate_simplex 1 ul ) / (sqrt(&2))) 
   else {}`;;

let mcell2 = new_definition `!V ul. mcell2 V ul = 
  if (hl (truncate_simplex 1 ul) < sqrt(&2)) /\ (sqrt(&2) <= hl ul)
  then
     (let a =  hl (truncate_simplex 1 ul ) / (sqrt(&2)) in
    rcone_ge (HD ul) (HD (TL ul)) a INTER
    rcone_ge  (HD (TL ul)) (HD ul) a INTER
    aff_ge { HD ul, (HD (TL ul)) } { mxi V ul, omega_list_n V ul 3 })
  else {}`;;

let mcell3 = new_definition `!V (ul:(real^3)list). mcell3 V ul = 
    if (hl (truncate_simplex 2 ul) < sqrt(&2)) /\ (sqrt(&2) <= hl ul) then
        convex hull (set_of_list (truncate_simplex 2 ul) UNION { mxi V ul })
    else {}`;;

let mcell4 = new_definition `!V ul. mcell4 (V:real^3->bool) (ul:(real^3)list) = 
 if (hl ul < sqrt(&2)) 
 then convex hull (set_of_list ul) 
 else {}`;;

let mcell = new_definition `!(V:real^3->bool) ul.  (mcell i V ul) = 
   if (i=0) then mcell0 V ul else if (i=1) then mcell1 V ul else if (i=2) then mcell2 V ul 
   else if (i=3) then mcell3 V ul else mcell4 V ul`;;

let mcell_set = new_definition `mcell_set V =
{ X |  ?i ul.  (X = mcell i V ul) /\ ul IN barV V 3 }`;;


let cell_params = new_definition `!V X. cell_params V X = @(k,ul).
    (k <= 4) /\ (ul IN barV V 3) /\ (X = mcell k V ul)`;;

let VX = new_definition
 `!V X.
         VX V X =
         (if NULLSET X
          then {}
          else let k,ul = cell_params V X in
               (if k = 0 then {} else
                set_of_list (truncate_simplex (k - 1) ul)))`;;

let edgeX = new_definition `!V X. edgeX V X =  { ({u,v}) | VX V X u /\ VX V X v /\ ~(u=v) }`;;

let total_solid = new_definition `!V X x. total_solid V X =  sum  ( VX V X  )  (\(x:real^3).   (sol x X))`;;

(* The definition of dihX is rather awkward, because a cell is essentially a quotient
    type, so it is not natural to make definitions with cells as a parameter.  
    Several lemmas will have to be proved to make these definitions useable, such as
   dihX V (mcell 4 V ul) ((EL 0 ul),(EL 1 ul)) = dihV (EL 0 ul) (EL 1 ul) (EL 2 ul) (EL 3 ul).
*)


(* deprecated 1/16/2012, tchales
let dihX2 = new_definition `!V X v0 (v1:real^3). dihX2 V X (v0,v1) = 
  (let (k,ul) = cell_params V X in
     dihV v0 v1 (mxi V ul) (omega_list_n V ul 3))`;;
p
(* The sum has a single term *)

let dihX3 = new_definition `!V X v0 (v1:real^3). dihX3 V X (v0,v1) =
  (let (k,ul) = cell_params V X in
     sum { (vl,p) | p permutes {0,1,2} /\ 
         (vl = left_action_list p ul) /\ (v0 = EL 0 vl) /\ (v1 = EL 1 vl) }
         ( \ (vl,p) . dihV v0 v1 (EL 2 vl) (mxi V ul) ))`;;

(* The sum has two terms. Take the average *)

let dihX4 = new_definition `!V X v0 (v1:real^3). dihX4 V X (v0,v1) =
  (let (k,ul) = cell_params V X in
     (&1 / &2) * sum { (vl,p) | p permutes {0,1,2,3} /\ 
         (vl = left_action_list p ul) /\ (v0 = EL 0 vl) /\ (v1 = EL 1 vl) }
         ( \ (vl,p). dihV v0 v1 (EL 2 vl) (EL 3 ul) ))`;;

let dihX = new_definition `!V X v0 (v1:real^3). dihX V X (v0,v1) = 
   if (NULLSET X) then &0 else
     (let (k,ul) = cell_params V X in
	if (k=2) then dihX2 V X (v0,v1)
	else if (k=3) then dihX3 V X (v0,v1)
	else if  (k=4) then dihX4 V X (v0,v1)
	else &0)`;;
*)

let dihu2 =  new_definition `!V ul. dihu2 V ul = 
   dihV (EL 0 ul) (EL 1 ul) (mxi V ul) (omega_list_n V ul 3)`;;

let dihu3 = new_definition `!V ul. dihu3 V ul = 
   dihV (EL 0 ul) (EL 1 ul) (EL 2 ul) (mxi V ul)`;;

let dihu4 = new_definition `!(ul:((real^3)list)). dihu4 ul = 
   dihV (EL 0 ul) (EL 1 ul) (EL 2 ul) (EL 3 ul)`;;

let cell_params_d = new_definition `!V X. cell_params_d V X vl = @(k,ul).
   (k <= 4) /\ (ul IN barV V 3) /\ (X = mcell k V ul) /\
   (initial_sublist vl ul)`;;

let dihX = new_definition `!V X v0 v1. dihX V X (v0,v1) = 
   if (NULLSET X) then &0 else
     (let (k,ul) = cell_params_d V X [v0;v1] in
	if (k=2) then dihu2 V ul
	else if (k=3) then dihu3 V ul
	else if  (k=4) then dihu4 ul
	else &0)`;;


   let sol0 = new_definition `sol0 = &3 * acs (&1 / &3)  - pi`;;
   let tau0 = new_definition `tau0 = &4 * pi - &20 * sol0`;;
   let mm1 = new_definition `mm1 = sol0 * sqrt(&8)/tau0`;;
   let mm2 = new_definition `mm2 = (&6 * sol0 - pi) * sqrt(&2) /(&6 * tau0)`;;
   let hplus = new_definition `hplus = #1.3254`;;

(*
   let marchal_quartic = new_definition `marchal_quartic h = 
      (sqrt(&2)-h)*(hplus - h)*(&17*h - &9*(h pow 2) - &3)/
      ((sqrt(&2) - &1)*(hplus - &1))`;;
*)

   let marchal =  new_definition `marchal h =
    (if (h <= sqrt(&2)) then marchal_quartic h else &0)`;;

(* changed March 16, 2012, following Ky's suggestion. *)

let gammaX = new_definition `!V X f.
    gammaX V X f = vol(X) - (&2*mm1 / pi)*total_solid V X +
     (&8*mm2/ pi) * sum (edgeX V X) (\{u,v}. (if {u,v} IN edgeX V X then (dihX V X (u ,v))* f(hl [u;v]) else (&0)))`;;

(*
   let gammaX = new_definition `!V X f.
    gammaX V X f = vol(X) - (&2*mm1 / pi)*total_solid V X +
     (&8*mm2/ pi) * sum (edgeX V X) (\{u,v}. (dihX V X (u ,v))* f(hl [u;v]))`;;
*)

(* OBSOLETE, Dec 31, 2012:
let marchal_inequality = new_definition `marchal_inequality =
    (!V X.  saturated V /\ packing V /\ mcell_set V X ==>
       gammaX V X marchal >= &0)`;;
*)

let h0 = new_definition `h0 = #1.26`;;

let lfun = new_definition `lfun h =  (h0 - h)/(h0 - &1)`;;

let lmfun = new_definition `lmfun h = if (h<=h0) then (h0 - h)/(h0 - &1) else &0`;;

let hminus = new_definition `hminus = @x. #1.2 <= x /\ x < #1.3 /\ marchal_quartic x = lmfun x`;;

let critical_edgeX =new_definition
   `critical_edgeX V X = { {u,v} | {u,v} IN edgeX V X /\ hminus <= hl [u;v] /\
					    hl[u;v] <= hplus}`;;

let subcritical_edgeX =new_definition
   `subcritical_edgeX V X = { {u,v} | {u,v} IN edgeX V X /\  hl [u;v] < hminus }`;;

let critical_weight = new_definition 
  `! V X. critical_weight V X = &1/ &(CARD (critical_edgeX V X))`;;

let bump = new_definition `!h. bump h = #0.005*(&1 - ((h- h0) pow 2)/((hplus - h0) pow 2))`;;

let critical_edge_y = new_definition `critical_edge_y y = ((&2*hminus <= y) /\ (y <= &2 *hplus))`;;

(* deprecated Dec 1, 2012 
let beta_bump_y = new_definition `beta_bump_y y1 y2 y3 y4 y5 y6 =
  (if critical_edge_y y1 then &1 else &0) *
  (if critical_edge_y y2 then &0 else &1) *
  (if critical_edge_y y3 then &0 else &1) *
  (if critical_edge_y y4 then &1 else &0) *
  (if critical_edge_y y5 then &0 else &1) *
  (if critical_edge_y y6 then &0 else &1) * 
    (bump (y1/ &2) - bump (y4 / &2))`;;
*)

let wtcount3_y = new_definition `wtcount3_y y1 y2 y3  = 
  (if critical_edge_y y1 then 1 else 0) +
  (if critical_edge_y y2 then 1 else 0) +
  (if critical_edge_y y3 then 1 else 0) `;;

let wtcount6_y = new_definition 
 `wtcount6_y y1 y2 y3 y4 y5 y6 = wtcount3_y y1 y2 y3 + wtcount3_y y4 y5 y6`;;

let cell_cluster = new_definition 
`!V e. cell_cluster V e = { X |  e IN critical_edgeX V X /\ mcell_set V X }`;;

(* Next definitions deprecated, Dec 1, 2012.
   Obsolete Jan 2, 2013. *)

(*
let beta_bump = new_definition `!V e X. beta_bump V e X = 
  (let (k,ul) = cell_params V X in  
    sum { (e',e'',p,vl) | (k=4) /\ (critical_edgeX V X = {e',e''}) /\ (e =e') /\ p permutes 0..3 
               /\ (vl = left_action_list p ul) /\ (e' = {EL 0 vl, EL 1 vl}) /\ (e'' = {EL 2 vl, EL 3 vl} ) }
    ( \ (e',e'',p,vl). (bump(hl [EL 0 vl; EL 1 vl]) - bump(hl [EL 2 vl; EL 3 vl]) )/ &4 ))`;;

let cluster_gammaX = new_definition `!V e Z. cluster_gammaX V e Z= 
   sum Z ( \ X.  gammaX V X lmfun  * critical_weight V X  + beta_bump V e X )`;;

let cell_cluster_estimate = new_definition `!V. cell_cluster_estimate V = (!e.  
  (cluster_gammaX V e (cell_cluster V e) >= &0))`;;

  (* average over all 4 = 2! x 2! possible representations *)

let beta_bumpA = new_definition `!V e X. beta_bumpA V e X = 
  (let (k,ul) = cell_params V X in  
    sum { (e',e'',p,vl) | (k=4) /\ (critical_edgeX V X = {e',e''}) /\ (e =e') /\ p permutes 0..3 /\
       (!e'''. e''' IN edgeX V X ==> ((e''' = e') \/ (e''' = e'') \/ (e''' IN subcritical_edgeX V X)))
               /\ (vl = left_action_list p ul) /\ (e' = {EL 0 vl, EL 1 vl}) /\ (e'' = {EL 2 vl, EL 3 vl} ) }
    ( \ (e',e'',p,vl). (bump(hl [EL 0 vl; EL 1 vl]) - bump(hl [EL 2 vl; EL 3 vl]) )/ &4 ))`;;

let cluster_gamma_AX = new_definition `!V e Z. cluster_gamma_AX V e Z= 
   sum Z ( \ X.  gammaX V X lmfun  * critical_weight V X  + beta_bumpA V e X )`;;

let cell_cluster_estimate_A = new_definition `!V. cell_cluster_estimate_A V = (!e.  
  (cluster_gamma_AX V e (cell_cluster V e) >= &0))`;;

*)

let beta_bump_v1 = new_definition `!V e X. beta_bump_v1 V e X = 
  let e' = VX V X DIFF e in
    (
      if (X IN mcell_set V /\ ~NULLSET X /\ e IN critical_edgeX V X /\ e' IN critical_edgeX V X /\
	    (!f. f IN edgeX V X ==> (f = e) \/ (f = e') \/ (f IN subcritical_edgeX V X)))
      then bump (radV e) - bump (radV e')
      else (&0)
    )`;;

let cluster_gamma_v1 = new_definition `!V e Z. cluster_gamma_v1 V e Z= 
   sum Z ( \ X.  gammaX V X lmfun  * critical_weight V X  + beta_bump_v1 V e X )`;;

let cell_cluster_estimate_v1 = new_definition `!V. cell_cluster_estimate_v1 V = (!e.  
  (cluster_gamma_v1 V e (cell_cluster V e) >= &0))`;;

let lmfun_inequality = new_definition 
  `!(V:real^N->bool). lmfun_inequality V = (!u.  u IN V ==>
     sum { v | v IN V /\ ~(u=v) /\ dist(u,v) <= &2*h0 } ( \v. lmfun (hl [u;v]) ) <= &12)`;;

let ball_annulus = new_definition
  `ball_annulus = cball((vec 0:real^3) , &2 * h0) DIFF ball(vec 0, &2)`;;

let lmfun_ineq_center = new_definition `!(V:real^N->bool). lmfun_ineq_center V = 
     sum V ( \v. lmfun (hl [vec 0;v]) ) <= &12`;;

let fan_of_polyhedron = new_definition `fan_of_polyhedron s = 
   { (v:real^3) | v extreme_point_of s } , 
   { {v,w}  | ~(v=w) /\ (convex hull {v,w}) face_of s }`;;


 end;;