glpk/tame_archive/lpproc.ml

   1 (* ========================================================================== *)
   2 (* FLYSPECK - BOOK FORMALIZATION                                              *)
   3 (*                                                                            *)
   4 (* Lemma: Easy Linear Programs     *)
   5 (* Chapter: Tame Hypermap                  *)
   6 (* Author: Thomas C. Hales                                                    *)
   7 (* Date: 2009-06-25                 *)
   8 (* ========================================================================== *)
   9
  10 (*
  11
  12 This file treats all but a few "hard" cases (hard_bb).
  13
  14 needs new mktop on platforms that do not support dynamic loading of Str.
  15
  16 ocamlmktop unix.cma nums.cma str.cma -o ocampl
  17 ./ocampl
  18
  19 glpk needs to be installed, and glpsol needs to be found in the path.
  20
  21 Execution instructions.
  22
  23 1- Download tame_archive_svn1830.txt from
  24      http://weyl.math.pitt.edu/hanoi2009/Kepler/Kepler.
  25      and save it as /tmp/tame_graph.txt
  26
  27      (* Update: 2011-5-9, This has been added to the project as tame_archive/string_archive.txt *)
  28
  29 2- make body.mod (using Parse_ineq.lpstring) if it doesn't exist already
  30      and place it in the directory glpk/tame_archive/
  31
  32   (* Update: 2012-4-8, This is part of the svn repository now.  It doesn't need to be regenerated. *)
  33
  34 3- Load this module and execute the following command.
  35
  36 let  (tame_bb,feasible_bb,hard_bb,easy_bb,remaining_easy_bb) = Lpproc.execute();;
  37
  38 *)
  39
  40 #load "str.cma";;
  41
  42 let flyspeck_dir =
  43   (try Sys.getenv "FLYSPECK_DIR" with Not_found -> Sys.getcwd());;
  44
  45 let project_root_dir = (Filename.concat (flyspeck_dir) Filename.parent_dir_name);;
  46
  47 let glpk_dir =  Filename.concat project_root_dir "glpk";;
  48
  49 let archive_dir = Filename.concat project_root_dir "graph_generator/output";;
  50
  51 let tame_dir =  Filename.concat glpk_dir "tame_archive";;
  52
  53 needs (Filename.concat glpk_dir "glpk_link.ml");;
  54
  55 module Lpproc = struct
  56
  57 open Glpk_link;;
  58 open List;;
  59
  60
  61 (* external files.
  62
  63    The archiveraw file can be downloaded from
  64    http://weyl.math.pitt.edu/hanoi2009/Kepler/Kepler (Archive of Tame Graphs).
  65
  66    body.mod is automatically generated by Parse_ineq.lpstring *)
  67
  68 (* let archiveraw = ref "/tmp/tame_graph.txt";;   (* read only *) *)
  69
  70 let archiveraw = ref (Filename.concat archive_dir "string_archive.txt");;
  71
  72 let modelbody = ref (Filename.concat tame_dir "body.mod");;
  73
  74 let model = Filename.concat tame_dir "graph_all.mod";;
  75
  76 (*
  77 let ampl_datafile = Filename.temp_file "ampl_datafile_" ".dat";;  (* only used for debugging purposes *)
  78 *)
  79
  80 let glpk_outfile = Filename.temp_file "glpk_outfile_" ".out";;
  81
  82 let make_model() =
  83   (Sys.chdir(tame_dir);
  84 Sys.command("cp head.mod "^model^"; cat "^(!modelbody)^"  >> "^
  85                      model^"; cat tail.mod >> "^model));;
  86
  87 (* conversion to list.  e.g. convert_to_ordered_list pentstring *)
  88
  89 (*
  90    The order in which the faces occur is the order of branch n bound.
  91    The order affects the effectiveness of the branch and bound.
  92    My heuristic is to branch on hexagons, then quads, then pents, then tris.
  93    Within faces of a given size, look for nodes that appear frequently.
  94 *)
  95 let order_list (h,xs) =
  96   let fl = flatten xs in
  97   let count k = length (filter ((=) k) fl) in
  98   let mc rs = maxlist0 (map count rs) in
  99   let sortfn a b = compare (mc b) (mc a) in
 100   let r k = filter (fun x -> length x = k) xs in
 101   let f k = sort sortfn (r k) in
 102    (h,f 6 @ f 4 @ f 5 @ f 3);;
 103
 104 let convert_to_ordered_list s = order_list (convert_to_list s);;
 105
 106 type hint  = Triangle_split of int list
 107               | High_low of int
 108               | Edge_split of int list;;
 109
 110 type init = No_data
 111             | File of string*Digest.t
 112             | Hash_tables of (((int,float) Hashtbl.t) * (((int*int),float) Hashtbl.t) * (((int*int),float) Hashtbl.t));;
 113
 114 (* change 4/2012 to make compatible with latest version of GLPK ( >= 4.47).
 115    New versions abort when problem is infeasible, rather than printing value 0.
 116 *)
 117
 118 type lptype = Lp_unset
 119               | Lp_infeasible
 120               | Lp_value of float;;
 121
 122 (* type for holding the parameters for a linear program.
 123     Many of the parameters are not needed for the easy cases that are treated
 124     in this module.  They will only be needed for the hard cases (hardid). *)
 125
 126 (*
 127    Fields for the hard cases: hints, diagnostic, node and edge lists
 128 *)
 129
 130 type branchnbound =
 131   {
 132     hypermap_id : string;
 133     mutable lpvalue : lptype;
 134     mutable hints : hint list;  (* hints about branching *)
 135     mutable diagnostic : init;
 136     string_rep : string;
 137     (* std_faces is the disjoint union of std_faces_not_super, std56_flat_free, std4_diag3 *)
 138     std_faces_not_super: int list list;
 139     std56_flat_free : int list list;
 140     std4_diag3 : int list list;
 141     std3_big : int list list;
 142     std3_small : int list list;
 143     (* special dart appears first in next ones *)
 144     apex_sup_flat : int list list;
 145     apex_flat : int list list;
 146     apex_A : int list list;
 147     apex5 : int list list;
 148     apex4 : int list list;
 149     (* edge lists *)
 150     d_edge_225_252 :  int list list;
 151     d_edge_200_225 :  int list list;
 152     (* node lists *)
 153     node_218_252 : int list;
 154     node_236_252 : int list;
 155     node_218_236 : int list;
 156     node_200_218 : int list;
 157   };;
 158
 159
 160
 161 let mk_bb s =
 162   let (h,face1) = convert_to_ordered_list s in
 163  {hypermap_id= h;
 164   lpvalue = Lp_unset;
 165   diagnostic = No_data;
 166   hints = [];
 167   string_rep=s;
 168   std_faces_not_super = face1;
 169   std3_big=[];
 170   std3_small=[];
 171   std4_diag3=[];
 172   std56_flat_free=[];
 173
 174   apex_flat=[];
 175   apex_sup_flat=[];
 176   apex_A=[];
 177   apex4=[];
 178   apex5=[];
 179
 180   d_edge_225_252=[];
 181   d_edge_200_225=[];
 182
 183   node_218_252=[];
 184   node_236_252=[];
 185   node_218_236=[];
 186   node_200_218=[];
 187  };;
 188
 189 let pent_bb = mk_bb pentstring;;
 190
 191 (* conversion to branchnbound.  e.g. mk_bb pentstring  *)
 192
 193 let modify_bb bb drop1std fields vfields =
 194   let add key xs t = nub ((get_values key xs) @ t)  in
 195   let jump_queue key xs vs =
 196     let ys = get_values key xs in
 197     let e = rotation ys in
 198       nub(ys @ (filter (fun t -> not(mem t e)) vs)) in
 199   let std = bb.std_faces_not_super in
 200   let jump_queue_std = jump_queue "jq" fields std in
 201 {
 202 hypermap_id = bb.hypermap_id;
 203 lpvalue = Lp_unset;
 204 diagnostic = No_data;
 205 hints = bb.hints;
 206 string_rep = bb.string_rep;
 207
 208 std_faces_not_super = if drop1std then tl std else jump_queue_std;
 209 std3_big = add "bt" fields bb.std3_big;
 210 std3_small = add "st" fields bb.std3_small;
 211 std56_flat_free = add "flat_free" fields bb.std56_flat_free;
 212 std4_diag3 = add "diag3" fields bb.std4_diag3;
 213
 214 apex_flat = add "ff" fields bb.apex_flat;
 215 apex_sup_flat = add "sf" fields bb.apex_sup_flat;
 216 apex_A = add "af" fields bb.apex_A;
 217 apex4 = add "apex4" fields bb.apex4;
 218 apex5 = add "apex5" fields bb.apex5;
 219
 220 d_edge_225_252 = add "e_225_252" fields bb.d_edge_225_252;
 221 d_edge_200_225 = add "e_200_225" fields bb.d_edge_200_225;
 222
 223 node_218_252 = add "218_252" vfields bb.node_218_252;
 224 node_236_252 = add "236_252" vfields bb.node_236_252;
 225 node_218_236 = add "218_236" vfields bb.node_218_236;
 226 node_200_218 = add "200_218" vfields bb.node_200_218;
 227 };;
 228
 229 (*
 230 Example: move [8;1;6;9] from std to std56_flat_free.
 231
 232 modify_bb pent_bb true  ["flat_free",[8;1;6;9];"ff",[9;10;11]] ["200_218",8;"218_252",3;"200_218",7];;
 233 pent_bb;;
 234 modify_bb pent_bb false ["jq",[0;3;5;4];"jq",[10;6;4]] [];;
 235 *)
 236
 237 (* functions on branch n bound *)
 238
 239 let std_faces bb = bb.std_faces_not_super @ bb.std56_flat_free @ bb.std4_diag3;;
 240 (*  @ bb.std3_big @ bb.std3_small;; *)
 241
 242 let std_tri_prebranch bb = filter
 243    (let r = rotation (bb.std3_big @ bb.std3_small) in
 244        fun t -> not(mem t r)  && (length t = 3)) bb.std_faces_not_super;;
 245
 246 (* should sort faces, so that numbering doesn't change so much when branching *)
 247
 248 let faces bb = (std_faces bb) @ bb.apex_sup_flat @ bb.apex_flat @
 249   bb.apex_A @ bb.apex5 @ bb.apex4;;
 250
 251 let triples w =
 252   let r j = nth w (j mod length w)  in
 253   let triple i =
 254       [r i; r (i+1); r(i+2)] in
 255     map triple (up (length w));;
 256
 257 let card_node bb =
 258   1+ maxlist0 (flatten (faces bb));;
 259
 260 let card_face bb = length(faces bb);;
 261
 262 let std_face_of_size bb r=
 263   let f = std_faces bb in
 264   let z = enumerate f in
 265     fst(split (filter (function _,y -> length y=r) z));;
 266
 267 let wheretriplemod xs  =
 268   let nth = List.nth in
 269   let t3 = (fun t -> [nth t 0;nth t 1;nth t 2]) in
 270   let ys = map t3 (rotation xs) in
 271   fun x ->
 272    try (
 273       whereis (t3 x) ys) mod (length xs)
 274    with Not_found -> failwith "wheretriplemod";;
 275
 276 (* generate ampl data string of branch n bound *)
 277
 278 let ampl_of_bb outs bb =
 279   let fs = faces bb in
 280   let where3 = wheremod fs in
 281   let number = map where3 in
 282   let list_of = unsplit " " string_of_int in
 283   let mk_faces xs = list_of (number xs) in
 284   let e_dart_raw  =
 285     map triples fs in
 286   let e_dart =
 287     let edata_row (i,x) = (sprintf "(*,*,*,%d) " i)^(unsplit ", " list_of x) in
 288       unsplit "\n" edata_row (enumerate e_dart_raw) in
 289   let mk_dart xs = sprintf "%d %d" (hd xs) (wheretriplemod fs xs) in
 290   let mk_darts xs = (unsplit ", " mk_dart xs) in
 291   let p = sprintf in
 292   let j = join_lines [
 293     p"param card_node := %d;" (card_node bb) ;
 294     p"param hypermap_id := %s;" bb.hypermap_id ;
 295     p"param card_face := %d;\n" (card_face bb);
 296     p"set std3 := %s;" (list_of (std_face_of_size bb 3)) ;
 297     p"set std4 := %s;" (list_of (std_face_of_size bb 4) );
 298     p"set std5 := %s;" (list_of (std_face_of_size bb 5)) ;
 299     p"set std6 := %s;\n" (list_of (std_face_of_size bb 6));
 300     p"set e_dart := \n%s;\n"  (e_dart);
 301     p"set std56_flat_free := %s;" (mk_faces bb.std56_flat_free);
 302     p"set std4_diag3 := %s;" (mk_faces bb.std4_diag3);
 303     p"set apex_sup_flat := %s;" (mk_darts bb.apex_sup_flat);
 304     p"set apex_flat := %s;" (mk_darts bb.apex_flat);
 305     p"set apex_A := %s;" (mk_darts bb.apex_A);
 306     p"set apex5 := %s;" (mk_darts bb.apex5);
 307     p"set apex4 := %s;" (mk_darts bb.apex4);
 308     p"set d_edge_225_252 := %s;" (mk_darts bb.d_edge_225_252);
 309     p"set d_edge_200_225 := %s;" (mk_darts bb.d_edge_200_225);
 310     p"set std3_big := %s;" (mk_faces bb.std3_big);
 311     p"set std3_small := %s;"  (mk_faces bb.std3_small);
 312     p"set node_218_252 := %s;" (list_of bb.node_218_252);
 313     p"set node_236_252 := %s;" (list_of bb.node_236_252);
 314     p"set node_218_236 := %s;" (list_of bb.node_218_236);
 315     p"set node_200_218 := %s;" (list_of bb.node_200_218)] in
 316     Printf.fprintf outs "%s" j;;
 317
 318 (* running of branch in glpsol *)
 319
 320 let solve_branch_verbose addhints bb =
 321   let set_some bb r = (* side effects *)
 322     if (length r = 1) then bb.lpvalue <- Lp_value (float_of_string(hd r)) else () in
 323   let (f,inp) = solve_branch_f model glpk_outfile "lnsum" ampl_of_bb bb in
 324   let _ = if (List.length f =0) then (set_some bb inp) else bb.lpvalue <- Lp_infeasible in
 325   let _ = bb.diagnostic <- File (glpk_outfile,Digest.file glpk_outfile) in
 326   let _ = addhints bb in (* hints for control flow *)
 327   let r = match bb.lpvalue with
 328     | Lp_unset -> "not set"
 329     | Lp_infeasible -> "infeasible"
 330     | Lp_value r -> sprintf "%3.3f" r in
 331   let _ = Sys.command(sprintf "echo %s: %s\n" bb.hypermap_id r) in
 332     bb;;
 333
 334 let solve_f f bb = match bb.lpvalue with
 335   | Lp_unset -> solve_branch_verbose f bb
 336   | Lp_value _ -> bb
 337   | Lp_infeasible -> bb;;
 338
 339 let solve bb = solve_f (fun t -> t) bb;;
 340
 341 (* filtering output *)
 342
 343 let is_feas bb =
 344   let feasible r = (r > 11.9999) in (* relax a bit from 12.0 *)
 345   match bb.lpvalue with
 346     | Lp_unset -> failwith "unexpected unset LP"
 347     | Lp_infeasible -> false
 348     | Lp_value r -> feasible r;;
 349
 350 let filter_feas_f f bbs =
 351   filter is_feas (map (solve_f f) bbs);;
 352
 353 let filter_feas bbs = filter_feas_f (fun t->t) bbs;;
 354
 355 (*
 356 branching on face data:
 357 switch_face does all the branching on the leading std face
 358 *)
 359
 360 let split_flatq xs i =  (* {y1,y3} is the new diagonal *)
 361   match rotateL i xs with
 362     | y1::y2::y3::ys  ->  ([y2;y3;y1],rotateR 1 (y1 :: y3 :: ys))
 363     | _ -> failwith "split_flatq match";;
 364
 365 let asplit_pent xs i = match rotateL i xs with
 366 (* y2,y4 darts of flat; y3 is the point of the A, {y1,y3}, {y3,y5} diags *)
 367   | y1::y2::y3::y4::[y5] ->  ([y2;y3;y1],[y3;y5;y1],[y4;y5;y3])
 368   | _ -> failwith "asplit_pent match error";;
 369
 370 let switch3 bb = match std_tri_prebranch bb with
 371   | [] -> failwith ("switch3 empty " ^ bb.hypermap_id)
 372   | fc::_  ->   [modify_bb bb false ["bt",fc] [];modify_bb bb false ["st",fc] []];;
 373
 374 let switch4 bb = match bb.std_faces_not_super with
 375   | [] -> failwith ("switch4 empty" ^ bb.hypermap_id)
 376   | fc::_  ->
 377   let mo s (a,b) = modify_bb bb true [s,a;s,b] [] in
 378   let f s i = mo s (split_flatq fc i) in
 379   let g s = modify_bb bb true [s,fc] [] in
 380   [f "ff" 0;f "ff" 1; f "sf" 0; f "sf" 1;g "diag3"];;
 381
 382 let switch5 bb = match  bb.std_faces_not_super with
 383   | [] -> failwith ("switch5 empty" ^ bb.hypermap_id)
 384   | fc::_ ->
 385   let mo (a,b) = modify_bb bb true ["ff",a;"apex4",b] [] in
 386   let f i = mo (split_flatq fc i) in
 387   let bbs = map f (up 5) in
 388   let mo (a,b,c) = modify_bb bb true ["ff",a;"af",b;"ff",c] [] in
 389   let f i = mo (asplit_pent fc i) in
 390   let ccs = map f (up 5) in
 391    (modify_bb bb true ["flat_free",fc] []) :: bbs @ ccs ;;
 392
 393 let switch6 bb = match bb.std_faces_not_super with
 394   | [] -> failwith ("switch6 empty" ^ bb.hypermap_id)
 395   | fc::_ ->
 396   let mo (a,b) = modify_bb bb true ["ff",a;"apex5",b] [] in
 397   let f i = mo (split_flatq fc i) in
 398    (modify_bb bb true ["flat_free",fc] []) :: (map f (up 6));;
 399
 400 let switch_face bb = match bb.std_faces_not_super with
 401   | [] -> [bb]
 402   | fc::_ ->
 403       let j = length fc in
 404       let fn = (nth [switch3;switch4;switch5;switch6] (j-3)) in
 405         fn bb;;
 406
 407 let echo bbs = Sys.command (sprintf "echo STACK %d %d" (length bbs) (nextc()));;
 408
 409 let onepass bbs =
 410   let branches = flatten (map switch_face bbs) in
 411   let _ = echo bbs in
 412     filter_feas branches;;
 413
 414 let rec allpass count bbs =
 415    let t = maxlist0 (map (fun b -> length (std_tri_prebranch b)) bbs) in
 416    if t = 0 or count <= 0 then bbs else allpass (count - 1) (onepass bbs);;
 417
 418 let hardid =
 419 [
 420 (* {3,3,3,3,3,3} *) "161847242261";
 421 (* {3,3,3,3,3,3} *) "223336279535";
 422 (* {3,3,3,3,3,3} *) "86506100695";
 423 (* one quad {3,3,4} *) "179189825656";
 424 (* two quad {3,4,4} *) "154005963125";
 425 (* {4,4,4} *) "65974205892";
 426 (* added 2010-06-24 *)
 427 "50803004532";"39599353438";"242652038506";
 428 "88089363170";"75641658977";"34970074286";
 429 (* added 2011-05-18 *)
 430 "164470574315"; "100126458338"; "215863975889"
 431 ];;
 432
 433
 434 (*
 435 tame_bb is the entire archive.
 436 feasible_bb are those that pass the first elementary linear program.
 437 hard_bb are those known to be difficult, and easy_bb is the complement.
 438 remaining_easy_bb should be empty; the easy ones should be entirely treated.
 439 *)
 440
 441 (* execute() takes about 2.5 hours to run.
 442
 443 Tested on 2010-06-24, svn 1847.
 444 with graph archive file tame_archive_svn1830.txt (June 17, 2010)
 445 It does not work on older versions of the archive because of differences
 446 in the hash code identifiers of graphs.
 447
 448 Retested 2010-07-28,
 449
 450 with graph archive file tame_archive_svn1830.txt (June 17, 2010)
 451 Digest.file "/tmp/tame_graph.txt";;
 452 val it : Digest.t = "X\221\153\0263Z\241\178\188\211S'\244\148f@"
 453
 454 project svn 1909
 455 body.mod Last Changed Rev: 1849
 456 lpproc.ml Last Changed Rev: 1850
 457
 458 Retested 2011-05-09, after deleting inequality 7676202716 from body.mod.
 459 Runs with 24K cases. All still good.
 460
 461 Retested 2011-05-15 on hex cases, using 0.6 instead of 0.7578=tameTableD[6,0].
 462 All still good.  (But a change in the tameTable would create more tame graphs,
 463 and this didn't terminate in reasonable time.)
 464
 465 *)
 466
 467 let execute() =
 468   let _ = make_model() in
 469   let _ = resetc() in
 470   let tame = strip_archive (!archiveraw) in
 471   let tame_bb = map mk_bb tame in
 472   let feasible_bb =  filter_feas (map solve tame_bb) in
 473   let (hard_bb,easy_bb) = partition (fun t -> (mem t.hypermap_id hardid)) feasible_bb in
 474   let remaining_easy_bb = allpass 20 easy_bb in
 475   (tame_bb,feasible_bb,hard_bb,easy_bb,remaining_easy_bb);;
 476
 477 end;;