glpk/tame_archive/hard_lp.ml

   1 (* ========================================================================== *)
   2 (* FLYSPECK - CODE FORMALIZATION                                              *)
   3 (*                                                                            *)
   4 (* Program: Linear Programs                                                   *)
   5 (* Author: Thomas C. Hales                                                    *)
   6 (* Date: 2010-08-01                                                           *)
   7 (* ========================================================================== *)
   8
   9
  10 (*
  11 The purpose of this file is to treat the 12 hard tame hypermaps that
  12 are not treated in lpproc.ml.
  13 *)
  14
  15 module Hard_lp = struct
  16
  17 (* code for the hard cases... *)
  18
  19 let hardidref = ref Lpproc.hardid;;
  20
  21 open Str;;
  22 open List;;
  23 open Glpk_link;;
  24 open Lpproc;;
  25
  26 let sqrt = Pervasives.sqrt;;
  27 let dih_y = Sphere_math.dih_y;;
  28
  29 (*
  30 The purpose of this first section is to build up tables of what the
  31 LP solution gives as the dihedral angles
  32 versus the floating point calculation of dihedral angles.
  33
  34 The idea is that branching should occur where the discrepancy is the biggest.
  35 *)
  36
  37 (* build up hashtables of all the variables assignments from the glpk_outfile *)
  38
  39 let new_tables () = (Hashtbl.create 13,Hashtbl.create 70,Hashtbl.create 70);;
  40
  41 let tables bb = match bb.diagnostic with
  42     Hash_tables (a,b,c) -> (a,b,c);;
  43
  44 let ynhash bb= (fun (a,_,_) -> a) (tables bb);;
  45 let yehash bb = (fun (_,a,_) -> a) (tables bb);;
  46 let azimhash bb = (fun (_,_,a) -> a) (tables bb);;
  47
  48 let yn bb i = Hashtbl.find (ynhash bb) i;;
  49 let ye bb (i,j) = Hashtbl.find (yehash bb) (i,j);;
  50 let azim bb (i,j) = Hashtbl.find (azimhash bb) (i,j);;
  51
  52 (*
  53 There is state involved here.  This code is rather fragile.
  54 The glpk_outfile gets overwritten the next time solve is called.
  55 So we want to generate the hash tables right after solve is called.
  56 *)
  57
  58 (* this renews the glpk_outfile *)
  59 let resolve bb = solve_branch_verbose (fun t->t) bb;;
  60
  61 let init_hash_verified_glpk_outfile glpk_outfile bb =
  62   let _ = (bb.diagnostic <- Hash_tables (new_tables())) in
  63   let _ = Hashtbl.clear (ynhash bb) in
  64   let _ = Hashtbl.clear (yehash bb) in
  65   let _ = Hashtbl.clear (azimhash bb) in
  66   let com = sprintf "cat %s | grep -v ':'  | grep '=' | tr '[\[\]=,]' ' ' | sed 's/.val//' | sed 's/\( [0-9]*\)$/\1.0/g'" glpk_outfile in
  67   let is = int_of_string in
  68   let fs = float_of_string in
  69   let (ic,oc) = Unix.open_process(com) in
  70   let _  = close_out oc in
  71   let inp = load_and_close_channel false ic in
  72   let _ = Unix.close_process(ic,oc) in
  73   let split_sp=  Str.split (regexp " +") in
  74   let inpa = map split_sp inp in
  75   let ynproc [a;b;c] = Hashtbl.add (ynhash bb) (is b) (fs c) in
  76   let yeproc [a;b;c;d] = Hashtbl.add (yehash bb) ( (is b), (is c)) (fs d) in
  77   let azimproc [a;b;c;d] = Hashtbl.add (azimhash bb) ( (is b), (is c)) (fs d) in
  78   let proc1 xs =
  79     begin
  80       let a = hd xs in
  81       if (a = "yn") then ynproc xs
  82       else if (a = "ye") then yeproc xs
  83       else if (a = "azim") then azimproc xs
  84      end in
  85   let _ = try (map proc1 inpa) with Failure _ -> failwith "init_has_v" in ();;
  86
  87 let init_hash bb =
  88   match bb.diagnostic with
  89     | File (glpk_outfile,dig) ->  if not(dig = Digest.file glpk_outfile)
  90       then (bb.diagnostic <- No_data; failwith "init_hash:stale")
  91       else init_hash_verified_glpk_outfile glpk_outfile bb
  92     | _ -> ();;
  93
  94 (* look at edge lengths and azimuth angles of a triangle *)
  95 let int_of_face xs bb = wheremod (faces bb) xs;;
  96
  97 (* get_azim_table dih_y [2;4;3] bb;; *)
  98
  99 let get_azim_dart_diff f xs bb =
 100    let [y1;y2;y3] = map (yn bb) xs in
 101    let [y6;y4;y5] = map (fun i -> (ye bb (i,int_of_face xs bb))) xs in
 102    let a1 = azim bb (hd xs, int_of_face xs bb) in
 103    let b1 = f y1 y2 y3 y4 y5 y6 in
 104    abs_float (a1 -. b1);;
 105 (* get_azim_dart_diff dih_y [2;4;3] bb;;  *)
 106
 107 (* first entry is the triangle with the biggest azim error. *)
 108 let sorted_azim_diff p bb =
 109   let ys = p bb in
 110   let u = map (fun t->  (get_azim_dart_diff dih_y t bb  ,t))  ys in
 111   let v = List.sort (fun a b -> - compare (fst a) (fst b)) u in
 112    v;;
 113 (* sorted_azim_diff darts_of_std_tri bb;;   *)
 114
 115 (* if we always do branching on the worst triangle, we quickly run out of
 116     possibilities because we can get stuck on a single triangle.  We need a
 117     heuristic that moves emphasis away from triangles where extensive branching
 118     has already been done.  Here is my heuristic to do so.   Many different
 119     strategies are possible.  This works sufficiently well in practice.  *)
 120
 121 let heuristic_weight d bb =
 122   if not(mem d (rotation (bb.std3_small @ bb.std3_big))) then 1.0 else
 123   if not(subset d (bb.node_200_218 @ bb.node_218_252)) then 0.7 else
 124   if not(subset d (bb.node_200_218 @ bb.node_218_236 @ bb.node_236_252)) then 0.49 else
 125   if not(subset (rotation [d]) (bb.d_edge_200_225 @ bb.d_edge_225_252 @ map (rotateL 1) bb.apex_flat)) then 0.3
 126   else 0.0;;
 127
 128 (* first entry is the triangle with the biggest azim weighted error. *)
 129
 130 let sorted_azim_weighted_diff p bb =
 131   let ys = p bb in
 132   let u = map (fun t->  ((heuristic_weight t bb *. get_azim_dart_diff dih_y t bb)  ,t))  ys in
 133   let v = List.sort (fun a b -> - compare (fst a) (fst b)) u in
 134    v;;
 135
 136
 137 (* ------------------------ *)
 138 (* HINTS *)
 139
 140 (* add_hints is called right after the LP is solved and the lpvalue set.
 141     The glpk_outfile has been initialized. *)
 142
 143 let darts_of_std_tri bb =
 144   rotation(filter (fun t -> length t = 3) bb.std_faces_not_super);;
 145
 146 let darts_of_std_tri_and_flats bb =
 147   rotation(filter (fun t -> length t = 3) (bb.std_faces_not_super @ bb.apex_flat));;
 148
 149 let highish bb = subtract bb.node_218_252 (bb.node_218_236 @ bb.node_236_252);;
 150
 151 let face_of_dart fc bb =
 152   let xss = map triples (faces bb) in
 153   nth (find (fun t -> mem fc t) xss) 0;;
 154
 155 let add_hints_force bb =
 156   try(
 157     let _ = init_hash bb in
 158   let dart = snd(hd(sorted_azim_weighted_diff darts_of_std_tri bb)) in
 159   let f = face_of_dart dart bb in
 160   if not(mem f (bb.std3_big @ bb.std3_small)) then bb.hints <-  [Triangle_split f] else
 161   let f1 = subtract f  (bb.node_200_218 @ bb.node_218_252) in
 162   if not(f1 = []) then bb.hints <- [High_low (hd f1)] else
 163     let f2 = intersect (highish bb) f in
 164   if not(f2 = []) then  bb.hints <- [High_low (hd f2)] else
 165     let d1 = subtract (rotation [dart]) (bb.d_edge_200_225 @ bb.d_edge_225_252) in
 166   if not (d1 = []) then bb.hints <- [Edge_split (hd d1)] else ()
 167   ) with | Failure s -> failwith (s^"in add_hints")
 168     | Not_found -> failwith "Not_found add_hints_force";;
 169
 170
 171 let add_hints_force_include_flat bb =
 172   try(
 173     let _ = init_hash bb in
 174     let dart = snd(hd(sorted_azim_weighted_diff darts_of_std_tri_and_flats bb)) in
 175     let f = face_of_dart dart bb in
 176       if (mem f (rotation bb.std_faces_not_super)) then add_hints_force bb else
 177         let f1 = subtract f  (bb.node_200_218 @ bb.node_218_252) in
 178           if not(f1 = []) then bb.hints <- [High_low (hd f1)] else
 179             let f2 = intersect (highish bb) f in
 180               if not(f2 = []) then  bb.hints <- [High_low (hd f2)] else
 181                 add_hints_force bb
 182   ) with | Failure s -> failwith ( s^" in add_hints_force_include_flat")
 183     | Not_found -> failwith ( " Not_found in add_hints_force_include_flat");;
 184
 185
 186 let add_hints_include_flat bb =
 187   if not(is_feas bb) then () else
 188   if not(bb.hints = []) then () else
 189     add_hints_force_include_flat bb;;
 190
 191 (* ------------------------ *)
 192
 193 let is_none bb = match bb.lpvalue with (* for debugging *)
 194   | Lp_value _ -> false
 195   | _ -> true;;
 196
 197
 198 let calc_max bbs = fold_right
 199   (fun bb x -> match bb.lpvalue with
 200     |Lp_value y -> max x y
 201     |_  -> x)
 202   bbs 0.0 ;;
 203
 204 let find_max bbs =
 205   let r = Lp_value (calc_max bbs) in
 206     find (fun bb -> r = bb.lpvalue) bbs;;
 207
 208 let findid s  = find (fun t -> s = t.hypermap_id);;
 209
 210 let findid_list s = filter (fun t -> s = t.hypermap_id);;
 211
 212
 213 (*
 214 We are proving the bound L(V) <= 12.  When there are 13 nodes,
 215 this is automatic when the nodes excess heights (|v|-2) add up to at least 0.52.
 216 Fail if the bb already satisfies L(V) <= 12.
 217 *)
 218
 219 let node_list bb = 0 -- (card_node bb - 1);;
 220
 221 let set_node_numerics bb =
 222   if not(card_node bb = 13) then bb else
 223   let n_high = length bb.node_236_252 in
 224   let n_mid = length bb.node_218_236 in
 225   let n_highish = length (highish bb) in
 226   if (n_high =0 )  & (n_mid +n_highish < 2) then bb else
 227   let _ = (n_mid * 18 + n_highish * 18 + n_high *36 <= 52) or failwith "set_node_numerics" in
 228   let node_new_low = subtract (node_list bb) (unions [bb.node_200_218 ;bb.node_218_236; bb.node_236_252;bb.node_218_252]) in
 229   let vfields_low = map (fun t -> ("200_218",t)) node_new_low in
 230   let vfields_mid = map(fun t->("218_236",t)) (highish bb) in
 231   let vfields = vfields_low @ vfields_mid in
 232     if vfields = [] then bb else     modify_bb bb false [] vfields;;
 233
 234 (*
 235 let t1 = modify_bb  (pent_bb) false [] ["236_252",5];;
 236 set_node_numerics t1;;
 237 let t1 = modify_bb  (pent_bb) false [] [("236_252",5);("218_236",3)];;
 238 can set_node_numerics t1;;
 239 let t1 = modify_bb  (pent_bb) false [] [("218_236",5);("218_236",3)];;
 240 set_node_numerics t1;;
 241 let t1 = modify_bb  (pent_bb) false [] [("218_236",5);("218_252",3)];;
 242 set_node_numerics t1;;
 243 *)
 244
 245 let opposite_edge [i;j;k] bb =
 246   let f = find (fun t -> (nth t 0 = j) & (nth t 1 = i)) (rotation (faces bb)) in
 247    [j;i;nth f 2];;
 248
 249 (*
 250 if a face is small the three edges lie in the range [2,2.25].
 251 if an edge is in the range (2.25,2.52], then the face is big.
 252 A directed edge is in the same category as the oppositely directed edge.
 253 *)
 254
 255 let set_face_numerics bb =
 256   let opp xs = nub (xs @ map (C opposite_edge bb) xs) in
 257   let edge_of_small = opp (rotation bb.std3_small) in
 258   let short_edge = opp bb.d_edge_200_225 in
 259   let long_edge = opp bb.d_edge_225_252 in
 260   let _ =  (intersect edge_of_small long_edge = []) or failwith "set_face_numerics" in
 261   let shortadds =  subtract (edge_of_small @ short_edge) bb.d_edge_200_225 in
 262   let shortfields = (map (fun t-> ("e_200_225",t)) shortadds) in
 263   let longadds =  subtract long_edge bb.d_edge_225_252 in
 264   let longfields = (map (fun t-> ("e_225_252",t)) longadds) in
 265   let r = filter (fun t -> mem t (std_faces bb) & (length t = 3) )
 266           (nub (map (C face_of_dart bb) long_edge)) in
 267   let _ = (intersect (rotation bb.std3_small) r =[]) or failwith "set_face_numerics" in
 268   let bigfields = map (fun t -> ("bt",t)) (subtract r bb.std3_big) in
 269   let fields = shortfields @ longfields @ bigfields in
 270     if fields=[] then bb else     modify_bb bb false fields [];;
 271
 272
 273 (*
 274 let t1 = modify_bb pent_bb false [("st",[9;6;10]);("e_225_252",[7;8;12])] [];;
 275 set_face_numerics t1;;
 276 let t1 = modify_bb pent_bb false ["st",[8;2;1]] [];;
 277 set_face_numerics t1;;
 278 let t1 = modify_bb pent_bb false ["e_225_252",[7;12;11]] [];;
 279 set_face_numerics t1;;
 280 *)
 281
 282 (*
 283 Hints are given as a list.  However, I never ended up using more
 284 than a single hint at a time.  It could be restructured as
 285 Some _ | None.
 286 *)
 287
 288 (* ------------------------ *)
 289 (* switch and selection functions *)
 290
 291 let clear_hint bb =
 292   if (bb.hints = []) then () else (bb.hints <- tl bb.hints);;
 293
 294 let switch_std3 d bb =
 295   let c = face_of_dart d bb in
 296   [modify_bb bb false ["bt",c] [];modify_bb bb false ["st",c] []];;
 297
 298 let switch_edge d bb =
 299       [modify_bb bb false ["e_225_252",d] [];modify_bb bb false ["e_200_225",d] []];;
 300
 301 let switch_node bb i =
 302   if (mem i (highish bb)) then
 303   [modify_bb bb false [] ["218_236",i];modify_bb bb false [] ["236_252",i]] else
 304     let settable = subtract (node_list bb ) (bb.node_200_218 @ bb.node_218_236 @ bb.node_236_252) in
 305   if not(mem i settable) then failwith "switch_node" else
 306   [modify_bb bb false [] ["218_252",i];modify_bb bb false [] ["200_218",i]];;
 307
 308 let follow_hint bb =
 309   let _ =assert(not(bb.hints = [])) in
 310   let hint = hd (bb.hints) in
 311 (*  let _ = clear_hint bb in *)
 312   let sbb = match hint with
 313   | Triangle_split d -> switch_std3 d bb
 314   | High_low i -> switch_node bb i
 315   | Edge_split d -> switch_edge d bb in
 316   let _ = map clear_hint sbb in sbb;;
 317
 318 let filter_feas_hint_include_flat bbs = filter_feas_f add_hints_include_flat bbs;;
 319
 320 let switch_hint bb =
 321    if (length bb.std_faces_not_super > 0) &
 322       (length (hd bb.std_faces_not_super) > 3) then switch_face bb
 323    else if not(bb.hints = []) then follow_hint bb else [bb];;
 324
 325 (* For debugging purposes, when we interrupt a calculation *)
 326
 327 let onepass_backup = ref [];;
 328
 329 let sortbb bbs =
 330   let eval bb = (match bb.lpvalue with Lp_value r -> r | _ -> 0.0) in
 331   let v = List.sort (fun a b -> - compare (eval a) (eval b)) bbs in
 332    v;;
 333
 334
 335
 336 (* One iteration of the main loop *)
 337
 338 let onepass_hint_include_flat = function
 339   [] -> []
 340   | bb::bbss as bbs ->
 341   let _ = onepass_backup := bbs in
 342   let brs =  switch_hint bb in
 343   let brs1 = map set_face_numerics brs in
 344   let brs2 = map set_node_numerics brs1 in
 345   let _ = echo bbs in
 346     sortbb ((filter_feas_hint_include_flat brs2) @ bbss);;
 347
 348
 349 (* The main loop *)
 350
 351 let rec allpass_hint_include_flat count bbs =
 352    if  count <= 0 then bbs else allpass_hint_include_flat (count - 1)
 353    (onepass_hint_include_flat bbs);;
 354
 355 let hard_string_rep() =
 356    List.find_all (fun s -> mem (fst (Glpk_link.convert_to_list s)) !hardidref)
 357    (Glpk_link.strip_archive (!Lpproc.archiveraw));;
 358
 359 let resolve_with_hints_include_flat t =
 360   let u = resolve t in
 361   let _ = add_hints_force_include_flat u in
 362     u;;
 363
 364 let hard_bb() =
 365   let r = map mk_bb (hard_string_rep()) in
 366   map resolve_with_hints_include_flat r;;
 367
 368 let hard i = List.nth (hard_bb()) i;;
 369
 370 (* if successful, all lists will be empty.  This takes several hours
 371     to run on my laptop.  *)
 372
 373 let execute() =
 374   let _ = Lpproc.make_model() in
 375   let _ = resetc() in
 376   map (fun t-> allpass_hint_include_flat 20000 [t]) (hard_bb());;
 377
 378 end;;