formal_lp/glpk/build_certificates.hl

   1 needs "../formal_lp/hypermap/main/lp_certificate.hl";;
   2 needs "../formal_lp/glpk/lp_binary_certificate.hl";;
   3 needs "../formal_lp/hypermap/computations/informal_computations.hl";;
   4 needs "../glpk/tame_archive/lpproc.ml";;
   5 needs "../glpk/tame_archive/hard_lp.ml";;
   6
   7
   8 module Build_certificates = struct
   9
  10
  11 (* Temporary directory where certificates are computed *)
  12 let tmp_dir = ref (flyspeck_dir ^ "/../formal_lp/glpk/tmp");;
  13 (* Output directory for certificates *)
  14 let output_dir = ref (flyspeck_dir ^ "/../formal_lp/glpk/binary");;
  15 (* Path to LP_HL.exe *)
  16 let lp_hl_dir = ref (flyspeck_dir ^ "/../formal_lp/LP-HL/LP-HL/bin/Release");;
  17 (* Main model file *)
  18 let model2_path = ref (tame_dir ^ "/model2.mod");;
  19
  20
  21
  22 open Glpk_link;;
  23 open Lp_informal_computations;;
  24 open Lpproc;;
  25 open Hard_lp;;
  26 open List;;
  27 open Lp_certificate;;
  28 open Lp_binary_certificate;;
  29
  30
  31 let mem_stat () =
  32   let stat = Gc.stat() in
  33   let word = float_of_int (Sys.word_size / 8) in
  34   let free = float_of_int stat.Gc.free_words *. word /. 1024.0 in
  35   let total = float_of_int stat.Gc.heap_words *. word /. 1024.0 in
  36   let allocated = total -. free in
  37   let str = sprintf "allocated = %f (free = %f; total_size = %f; %f)\n"
  38     allocated free total (free /. total) in
  39     print_string str;;
  40
  41
  42 let run_command dir com =
  43   let cur_dir = Sys.getcwd() in
  44   let _ = Sys.chdir dir in
  45   let result = Sys.command com in
  46   let _ = Sys.chdir cur_dir in
  47     result;;
  48
  49
  50 let make_models include_main =
  51   let _ = make_model() in
  52   let sed1 = "s/maximize objective:.*/maximize objective: sum{i in node} ln[i];/" and
  53       sed2 = "s/lnsum_def:.*//" and
  54       sed3 = "s/main:.*//" in
  55   let _ = if include_main then () else
  56     let cmd = sprintf "sed -e '%s' %s > %s" sed3 model !model2_path in
  57     let cmd2 = sprintf "cp %s %s" !model2_path model in
  58     let _ = run_command tame_dir cmd in
  59     let _ = run_command tame_dir cmd2 in
  60       () in
  61   let cmds = [sed1; sed2] in
  62   let str1 = itlist (fun cmd str -> sprintf "-e '%s' %s" cmd str) cmds "" in
  63   let cmd = sprintf "sed %s %s > %s" str1 model !model2_path in
  64   let _ = run_command tame_dir cmd in
  65     ();;
  66
  67
  68 let ampl_of_bb' fname bb =
  69   let out = open_out fname in
  70   let _ = ampl_of_bb out bb in
  71     close_out out;;
  72
  73
  74 let save_string fname str =
  75   let out = open_out fname in
  76   let _ = Printf.fprintf out "%s" str in
  77     close_out out;;
  78
  79
  80 let find_index p list = fst (find (fun _,x -> p x) (Glpk_link.enumerate list));;
  81
  82 let find_face_index fs f =
  83   let rots = Glpk_link.rotation [f] in
  84   let rec find fs =
  85     match fs with
  86       | [] -> failwith "find_face_index"
  87       | h :: t ->
  88           let eq = fold_right (or) (map ((=) h) rots) false in
  89             if eq then 0 else find t + 1 in
  90     find fs;;
  91
  92
  93 let build_permutation fs0 fs = map (fun f -> find_face_index fs f) fs0;;
  94
  95
  96 let save_info out_dir name infeasible hyp_list bb =
  97   let f2 = faces bb in
  98   let perm = build_permutation hyp_list f2 in
  99   let _ = ampl_of_bb' (sprintf "%s/%s_pars.txt" out_dir name) bb in
 100   let perm_str = unsplit ", " string_of_int perm in
 101   let hyp_str = unsplit ";" (unsplit "," string_of_int) hyp_list in
 102   let p = sprintf in
 103   let lines = [
 104     p "name: %s" name;
 105     p "infeasible: %b" infeasible;
 106     p "hypermap: %s" hyp_str;
 107     p "faces: %s" perm_str;
 108   ] in
 109     save_string (out_dir ^ "/flyspeck-" ^ name ^ ".txt") (join_lines lines);;
 110
 111
 112 let clean_out_dir =
 113   let permanent_files = [
 114     "000.txt";
 115     "string_archive.txt";
 116   ] in
 117     fun out_dir ->
 118       let files = Array.to_list (Sys.readdir out_dir) in
 119       let files' = filter (fun name -> not (mem name permanent_files)) files in
 120       let _ = map (fun name -> Sys.remove (sprintf "%s/%s" out_dir name)) files' in
 121         ();;
 122
 123 (************************************)
 124
 125 (* Builds a terminal case from the given hypermap *)
 126 let build_terminal_case =
 127   let create_infeasible_solution out_name =
 128     let sed1_cmd = sprintf "sed -f %s/../sed_add_slack.sed %s.lp > %s_slack.lp" !tmp_dir out_name out_name and
 129         sed2_cmd = sprintf "sed -f %s/../sed_build_objective.sed %s.lp > /dev/null" !tmp_dir out_name and
 130         sed3_cmd = sprintf "sed -f %s/../sed_replace_objective.sed %s_slack.lp > %s_slack2.lp" !tmp_dir out_name out_name in
 131     let _ = run_command !tmp_dir sed1_cmd and
 132         _ = run_command !tmp_dir sed2_cmd and
 133         _ = run_command !tmp_dir sed3_cmd in
 134     let solve_com2 = sprintf "glpsol --cpxlp %s_slack2.lp -w %s.txt" out_name out_name in
 135       run_command !tmp_dir solve_com2 in
 136
 137     fun hyp_list infeasible bb ->
 138       let error_glpsol = sprintf "build_terminal_case: glpsol failed for %s" bb.hypermap_id in
 139       let _ = clean_out_dir !tmp_dir in
 140       let out_name = "out" in
 141         (* Create info and parameter files *)
 142       let _ = save_info !tmp_dir out_name infeasible hyp_list bb in
 143       let data_file = sprintf "%s_pars.txt" out_name in
 144         (* Create .lp file and solve a feasible problem *)
 145       let solve_com = sprintf "glpsol -m %s -d %s -w %s.txt --wcpxlp %s.lp > /dev/null"  !model2_path data_file out_name out_name in
 146       let result = run_command !tmp_dir solve_com in
 147       let _ = result = 0 or failwith error_glpsol in
 148         (* Solve an infeasible problem *)
 149       let result = if infeasible then create_infeasible_solution out_name else 0 in
 150       let _ = result = 0 or failwith error_glpsol in
 151         (* Run LP-HL.exe *)
 152       let com = sprintf "mono %s/LP-HL.exe %s" !lp_hl_dir (sprintf "flyspeck-%s.txt" out_name) in
 153       let result = run_command !tmp_dir com in
 154       let _ = result = 0 or failwith (sprintf "LP-HL.exe failed for %s" bb.hypermap_id) in
 155         (* Create a certificate *)
 156       let terminal = read_binary_terminal (sprintf "%s/%s_out.bin" !tmp_dir out_name) in
 157         terminal;;
 158
 159
 160 (*************************************)
 161
 162 (* If true, then terminal cases are not constructed *)
 163 let test_mode = ref false;;
 164 let print_progress = ref false;;
 165
 166
 167 let reset_build_counters, next_build_case, next_terminal_case, report_build_progress =
 168   let cases = ref 0 and
 169       terminals = ref 0 in
 170     (fun () -> cases := 0; terminals := 0),
 171   (fun () -> cases := !cases + 1),
 172   (fun () -> terminals := !terminals + 1),
 173   (fun () -> if !print_progress then report (sprintf "cases = %d; terminals = %d" !cases !terminals) else ());;
 174
 175
 176 let set_face_numerics_info bb =
 177   let opp xs = nub (xs @ map (C opposite_edge bb) xs) in
 178   let edge_of_small = opp (rotation bb.std3_small) in
 179   let short_edge = opp bb.d_edge_200_225 in
 180   let long_edge = opp bb.d_edge_225_252 in
 181   let _ =  (intersect edge_of_small long_edge = []) or failwith "set_face_numerics" in
 182   let shortadds =  subtract (edge_of_small @ short_edge) bb.d_edge_200_225 in
 183   let shortfields = (map (fun t-> ("e_200_225",t)) shortadds) in
 184   let longadds =  subtract long_edge bb.d_edge_225_252 in
 185   let longfields = (map (fun t-> ("e_225_252",t)) longadds) in
 186   let r = filter (fun t -> mem t (std_faces bb) & (length t = 3) )
 187           (nub (map (C face_of_dart bb) long_edge)) in
 188   let _ = (intersect (rotation bb.std3_small) r =[]) or failwith "set_face_numerics" in
 189   let bigfields = map (fun t -> ("bt",t)) (subtract r bb.std3_big) in
 190   let fields = shortfields @ longfields @ bigfields in
 191     if fields=[] then bb, [] else
 192       let new_bb = modify_bb bb false fields [] in
 193         (* Don't need to add symmetric inequalities and short edges for small triangles: *)
 194         (* it is done at other steps. Add new big triangles only. *)
 195       let new_big_faces = subtract r bb.std3_big in
 196       let long_edge_faces = zip (map (C face_of_dart bb) long_edge) long_edge in
 197       let darts = map (C assoc long_edge_faces) new_big_faces in
 198       let info_list = map (fun d -> {split_type = "add_big"; split_face = d}) darts in
 199         new_bb, info_list;;
 200
 201
 202 let set_node_numerics_info bb =
 203   if not(card_node bb = 13) then bb, [] else
 204   let n_high = length bb.node_236_252 in
 205   let n_mid = length bb.node_218_236 in
 206   let n_highish = length (highish bb) in
 207   if (n_high =0 )  & (n_mid +n_highish < 2) then bb, [] else
 208   let _ = (n_mid * 18 + n_highish * 18 + n_high *36 <= 52) or failwith "set_node_numerics" in
 209   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
 210   let vfields_low = map (fun t -> ("200_218",t)) node_new_low in
 211   let vfields_mid = map(fun t->("218_236",t)) (highish bb) in
 212   let vfields = vfields_low @ vfields_mid in
 213     if vfields = [] then bb, [] else
 214       let new_bb = modify_bb bb false [] vfields in
 215       let info =
 216         if n_high > 0 then
 217           {split_type = "high"; split_face = bb.node_236_252}
 218         else
 219           {split_type = "mid"; split_face = bb.node_218_236 @ highish bb} in
 220         new_bb, [info];;
 221
 222
 223 let rec build hard_flag (hyp_list, bb0) =
 224   let _ = next_build_case() in
 225   let _ = report_build_progress() in
 226
 227   let bb, info_list = if not hard_flag then bb0, []
 228   else
 229     let bb1, list1 = set_face_numerics_info bb0 in
 230     let bb2, list2 = set_node_numerics_info bb1 in
 231       bb2, list1 @ list2 in
 232
 233   let f_hint = if hard_flag then add_hints_include_flat else (fun t -> ()) in
 234   let result = solve_f f_hint bb in
 235   let certificate =
 236     if not (is_feas result) then
 237       let _ = next_terminal_case() in
 238       let infeasible =
 239         match result.lpvalue with
 240           | Lp_infeasible -> true
 241           | _ -> false in
 242         (* terminal case *)
 243       let terminal =
 244         if !test_mode then
 245           empty_terminal
 246         else
 247           build_terminal_case hyp_list infeasible bb in
 248         Lp_terminal terminal
 249     else
 250       (* split case *)
 251       let n = length (hd bb.std_faces_not_super) in
 252         if n = 3 then
 253           if hard_flag then
 254             split3_hard hyp_list bb
 255           else
 256             split3 hard_flag hyp_list bb
 257         else if n = 4 then
 258           split4 hard_flag hyp_list bb
 259         else if n = 5 then
 260           split5 hard_flag hyp_list bb
 261         else if n = 6 then
 262           split6 hard_flag hyp_list bb
 263         else
 264           failwith (sprintf "build: incorrect face size - %d" n) in
 265     itlist (fun info c -> Lp_split (info, [c])) info_list certificate
 266
 267 (* split3 - hard cases *)
 268 and split3_hard hyp_list bb =
 269   let _ = assert (bb.hints <> []) in
 270   let bbs, split_type, list =
 271     match hd (bb.hints) with
 272       | Triangle_split d -> switch_std3 d bb, "tri", face_of_dart d bb
 273       | Edge_split d -> switch_edge d bb, "edge", d
 274       | High_low i -> switch_node bb i, (if mem i (highish bb) then "236" else "218"), [i] in
 275   let _ = map clear_hint bbs in
 276   let case_args = zip [hyp_list; hyp_list] bbs in
 277   let cases = map (build true) case_args in
 278   let info = {split_type = split_type; split_face = list} in
 279     Lp_split (info, cases)
 280
 281 (* split3 *)
 282 and split3 hard_flag hyp_list bb =
 283   let split_face = hd (std_tri_prebranch bb) in
 284   let _ = assert (length split_face = 3) in
 285   let bbs = switch3 bb in
 286   let case_args = zip [hyp_list; hyp_list] bbs in
 287   let cases = map (build hard_flag) case_args in
 288   let info = {split_type = "tri"; split_face = split_face} in
 289     Lp_split (info, cases)
 290
 291 (* split4 *)
 292 and split4 hard_flag hyp_list bb =
 293   let split_face = hd bb.std_faces_not_super in
 294   let _ = assert (length split_face = 4) in
 295   let bbs = switch4 bb in
 296   let dart13 = nth split_face 1, nth split_face 2 and
 297       dart24 = nth split_face 0, nth split_face 1 in
 298   let split13 = split_list hyp_list dart13 and
 299       split24 = split_list hyp_list dart24 in
 300   let case_args = zip [split13; split24; split13; split24; hyp_list] bbs in
 301   let cases = map (build hard_flag) case_args in
 302   let info = {split_type = "quad"; split_face = split_face} in
 303     Lp_split (info, cases)
 304
 305 (* split5 *)
 306 and split5 hard_flag hyp_list bb =
 307   let split_face = hd bb.std_faces_not_super in
 308   let _ = assert (length split_face = 5) in
 309   let bbs = switch5 bb in
 310   let darts = rotateL 1 (list_pairs split_face) in
 311   let splits_one = map (split_list hyp_list) darts in
 312   let splits_two = map2 split_list splits_one (rotateL 2 darts) in
 313   let case_args = zip (hyp_list :: (splits_one @ splits_two)) bbs in
 314   let cases = map (build hard_flag) case_args in
 315   let info = {split_type = "pent"; split_face = split_face} in
 316     Lp_split (info, cases)
 317
 318 (* split6 *)
 319 and split6 hard_flag hyp_list bb =
 320   let split_face = hd bb.std_faces_not_super in
 321   let _ = assert (length split_face = 6) in
 322   let bbs = switch6 bb in
 323   let darts = Glpk_link.rotateL 1 (list_pairs split_face) in
 324   let splits = map (split_list hyp_list) darts in
 325   let case_args = zip (hyp_list :: splits) bbs in
 326   let cases = map (build hard_flag) case_args in
 327   let info = {split_type = "hex"; split_face = split_face} in
 328     Lp_split (info, cases);;
 329
 330
 331 (* Moves all hex faces to std56_flat_free *)
 332 (* This operation can be done on all hypermaps since branching on hexes is not required: *)
 333 (* there are no inequalities for std56_flat_free INTER std6 *)
 334 let modify_hex_cases bb =
 335   let faces6 = filter (fun f -> length f = 6) bb.std_faces_not_super in
 336     itlist (fun fc bb -> modify_bb bb true ["flat_free", fc] []) faces6 bb;;
 337
 338
 339 (* Builds an lp certificate *)
 340 let build_certificate modify_hex bb =
 341   let _ = reset_build_counters() in
 342   let bb = if modify_hex then modify_hex_cases bb else bb in
 343   let hyp_list = (snd o convert_to_list) bb.string_rep in
 344   let hard_flag = bb.hints <> [] in
 345   let root = build hard_flag (hyp_list, bb) in
 346   let _ = report_build_progress() in
 347   let certificate = {hypermap_string = bb.string_rep; root_case = root} in
 348     certificate;;
 349
 350
 351
 352
 353 (************************)
 354
 355 (* Builds and saves certificates for all given cases *)
 356 let build_and_save_all =
 357   let counter = ref 0 and
 358       total = ref 0 in
 359   let buf = ref [] and
 360       file_counter = ref 0 and
 361       buf_terminals = ref 0 in
 362
 363   let save_buf name_prefix =
 364     if length !buf > 0 then
 365       let _ = file_counter := !file_counter + 1 in
 366       let fname = sprintf "%s/%s_%d.dat" !output_dir name_prefix !file_counter in
 367       let _ = write_lp_certificates fname !buf in
 368         buf := [];
 369         buf_terminals := 0
 370     else
 371       () in
 372
 373   let process name_prefix max modify_hex bb =
 374     let _ = counter := !counter + 1 in
 375     let _ = report (sprintf "%d/%d" !counter !total) in
 376     let certificate = build_certificate modify_hex bb in
 377     let _ = buf := certificate :: !buf in
 378     let _ = buf_terminals := !buf_terminals + count_terminals certificate.root_case in
 379       if !buf_terminals >= max then
 380         save_buf name_prefix
 381       else
 382         () in
 383
 384     fun name_prefix max modify_hex bbs ->
 385       let _ = total := length bbs and
 386           _ = counter := 0 and
 387           _ = file_counter := 0 and
 388           _ = buf := [] and
 389           _ = buf_terminals := 0 in
 390       let _ = map (process name_prefix max modify_hex) bbs in
 391         save_buf name_prefix;;
 392
 393
 394 end;;