needs "../formal_lp/hypermap/main/lp_certificate.hl";;
needs "../formal_lp/glpk/lp_binary_certificate.hl";;
needs "../formal_lp/hypermap/computations/informal_computations.hl";;
needs "../glpk/tame_archive/lpproc.ml";;
needs "../glpk/tame_archive/hard_lp.ml";;
module Build_certificates = struct
(* Temporary directory where certificates are computed *)
let tmp_dir = ref (flyspeck_dir ^ "/../formal_lp/glpk/tmp");;
(* Output directory for certificates *)
let output_dir = ref (flyspeck_dir ^ "/../formal_lp/glpk/binary");;
(* Path to LP_HL.exe *)
let lp_hl_dir = ref (flyspeck_dir ^ "/../formal_lp/LP-HL/LP-HL/bin/Release");;
(* Main model file *)
let model2_path = ref (tame_dir ^ "/model2.mod");;
open Glpk_link;;
open Lp_informal_computations;;
open Lpproc;;
open Hard_lp;;
open List;;
open Lp_certificate;;
open Lp_binary_certificate;;
let mem_stat () =
let stat = Gc.stat() in
let word = float_of_int (Sys.word_size / 8) in
let free = float_of_int stat.Gc.free_words *. word /. 1024.0 in
let total = float_of_int stat.Gc.heap_words *. word /. 1024.0 in
let allocated = total -. free in
let str = sprintf "allocated = %f (free = %f; total_size = %f; %f)\n"
allocated free total (free /. total) in
print_string str;;
let run_command dir com =
let cur_dir = Sys.getcwd() in
let _ = Sys.chdir dir in
let result = Sys.command com in
let _ = Sys.chdir cur_dir in
result;;
let make_models include_main =
let _ = make_model() in
let sed1 = "s/maximize objective:.*/maximize objective: sum{i in node} ln[i];/" and
sed2 = "s/lnsum_def:.*//" and
sed3 = "s/main:.*//" in
let _ = if include_main then () else
let cmd = sprintf "sed -e '%s' %s > %s" sed3 model !model2_path in
let cmd2 = sprintf "cp %s %s" !model2_path model in
let _ = run_command tame_dir cmd in
let _ = run_command tame_dir cmd2 in
() in
let cmds = [sed1; sed2] in
let str1 = itlist (fun cmd str -> sprintf "-e '%s' %s" cmd str) cmds "" in
let cmd = sprintf "sed %s %s > %s" str1 model !model2_path in
let _ = run_command tame_dir cmd in
();;
let ampl_of_bb' fname bb =
let out = open_out fname in
let _ = ampl_of_bb out bb in
close_out out;;
let save_string fname str =
let out = open_out fname in
let _ = Printf.fprintf out "%s" str in
close_out out;;
let find_index p list = fst (find (fun _,x -> p x) (Glpk_link.enumerate list));;
let find_face_index fs f =
let rots = Glpk_link.rotation [f] in
let rec find fs =
match fs with
| [] -> failwith "find_face_index"
| h :: t ->
let eq = fold_right (or) (map ((=) h) rots) false in
if eq then 0 else find t + 1 in
find fs;;
let build_permutation fs0 fs = map (fun f -> find_face_index fs f) fs0;;
let save_info out_dir name infeasible hyp_list bb =
let f2 = faces bb in
let perm = build_permutation hyp_list f2 in
let _ = ampl_of_bb' (sprintf "%s/%s_pars.txt" out_dir name) bb in
let perm_str = unsplit ", " string_of_int perm in
let hyp_str = unsplit ";" (unsplit "," string_of_int) hyp_list in
let p = sprintf in
let lines = [
p "name: %s" name;
p "infeasible: %b" infeasible;
p "hypermap: %s" hyp_str;
p "faces: %s" perm_str;
] in
save_string (out_dir ^ "/flyspeck-" ^ name ^ ".txt") (join_lines lines);;
let clean_out_dir =
let permanent_files = [
"000.txt";
"string_archive.txt";
] in
fun out_dir ->
let files = Array.to_list (Sys.readdir out_dir) in
let files' = filter (fun name -> not (mem name permanent_files)) files in
let _ = map (fun name -> Sys.remove (sprintf "%s/%s" out_dir name)) files' in
();;
(************************************)
(* Builds a terminal case from the given hypermap *)
let build_terminal_case =
let create_infeasible_solution out_name =
let sed1_cmd = sprintf "sed -f %s/../sed_add_slack.sed %s.lp > %s_slack.lp" !tmp_dir out_name out_name and
sed2_cmd = sprintf "sed -f %s/../sed_build_objective.sed %s.lp > /dev/null" !tmp_dir out_name and
sed3_cmd = sprintf "sed -f %s/../sed_replace_objective.sed %s_slack.lp > %s_slack2.lp" !tmp_dir out_name out_name in
let _ = run_command !tmp_dir sed1_cmd and
_ = run_command !tmp_dir sed2_cmd and
_ = run_command !tmp_dir sed3_cmd in
let solve_com2 = sprintf "glpsol --cpxlp %s_slack2.lp -w %s.txt" out_name out_name in
run_command !tmp_dir solve_com2 in
fun hyp_list infeasible bb ->
let error_glpsol = sprintf "build_terminal_case: glpsol failed for %s" bb.hypermap_id in
let _ = clean_out_dir !tmp_dir in
let out_name = "out" in
(* Create info and parameter files *)
let _ = save_info !tmp_dir out_name infeasible hyp_list bb in
let data_file = sprintf "%s_pars.txt" out_name in
(* Create .lp file and solve a feasible problem *)
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
let result = run_command !tmp_dir solve_com in
let _ = result = 0 or failwith error_glpsol in
(* Solve an infeasible problem *)
let result = if infeasible then create_infeasible_solution out_name else 0 in
let _ = result = 0 or failwith error_glpsol in
(* Run LP-HL.exe *)
let com = sprintf "mono %s/LP-HL.exe %s" !lp_hl_dir (sprintf "flyspeck-%s.txt" out_name) in
let result = run_command !tmp_dir com in
let _ = result = 0 or failwith (sprintf "LP-HL.exe failed for %s" bb.hypermap_id) in
(* Create a certificate *)
let terminal = read_binary_terminal (sprintf "%s/%s_out.bin" !tmp_dir out_name) in
terminal;;
(*************************************)
(* If true, then terminal cases are not constructed *)
let test_mode = ref false;;
let print_progress = ref false;;
let reset_build_counters, next_build_case, next_terminal_case, report_build_progress =
let cases = ref 0 and
terminals = ref 0 in
(fun () -> cases := 0; terminals := 0),
(fun () -> cases := !cases + 1),
(fun () -> terminals := !terminals + 1),
(fun () -> if !print_progress then report (sprintf "cases = %d; terminals = %d" !cases !terminals) else ());;
let set_face_numerics_info bb =
let opp xs = nub (xs @ map (C opposite_edge bb) xs) in
let edge_of_small = opp (rotation bb.std3_small) in
let short_edge = opp bb.d_edge_200_225 in
let long_edge = opp bb.d_edge_225_252 in
let _ = (intersect edge_of_small long_edge = []) or failwith "set_face_numerics" in
let shortadds = subtract (edge_of_small @ short_edge) bb.d_edge_200_225 in
let shortfields = (map (fun t-> ("e_200_225",t)) shortadds) in
let longadds = subtract long_edge bb.d_edge_225_252 in
let longfields = (map (fun t-> ("e_225_252",t)) longadds) in
let r = filter (fun t -> mem t (std_faces bb) & (length t = 3) )
(nub (map (C face_of_dart bb) long_edge)) in
let _ = (intersect (rotation bb.std3_small) r =[]) or failwith "set_face_numerics" in
let bigfields = map (fun t -> ("bt",t)) (subtract r bb.std3_big) in
let fields = shortfields @ longfields @ bigfields in
if fields=[] then bb, [] else
let new_bb = modify_bb bb false fields [] in
(* Don't need to add symmetric inequalities and short edges for small triangles: *)
(* it is done at other steps. Add new big triangles only. *)
let new_big_faces = subtract r bb.std3_big in
let long_edge_faces = zip (map (C face_of_dart bb) long_edge) long_edge in
let darts = map (C assoc long_edge_faces) new_big_faces in
let info_list = map (fun d -> {split_type = "add_big"; split_face = d}) darts in
new_bb, info_list;;
let set_node_numerics_info bb =
if not(card_node bb = 13) then bb, [] else
let n_high = length bb.node_236_252 in
let n_mid = length bb.node_218_236 in
let n_highish = length (highish bb) in
if (n_high =0 ) & (n_mid +n_highish < 2) then bb, [] else
let _ = (n_mid * 18 + n_highish * 18 + n_high *36 <= 52) or failwith "set_node_numerics" in
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
let vfields_low = map (fun t -> ("200_218",t)) node_new_low in
let vfields_mid = map(fun t->("218_236",t)) (highish bb) in
let vfields = vfields_low @ vfields_mid in
if vfields = [] then bb, [] else
let new_bb = modify_bb bb false [] vfields in
let info =
if n_high > 0 then
{split_type = "high"; split_face = bb.node_236_252}
else
{split_type = "mid"; split_face = bb.node_218_236 @ highish bb} in
new_bb, [info];;
let rec build hard_flag (hyp_list, bb0) =
let _ = next_build_case() in
let _ = report_build_progress() in
let bb, info_list = if not hard_flag then bb0, []
else
let bb1, list1 = set_face_numerics_info bb0 in
let bb2, list2 = set_node_numerics_info bb1 in
bb2, list1 @ list2 in
let f_hint = if hard_flag then add_hints_include_flat else (fun t -> ()) in
let result = solve_f f_hint bb in
let certificate =
if not (is_feas result) then
let _ = next_terminal_case() in
let infeasible =
match result.lpvalue with
| Lp_infeasible -> true
| _ -> false in
(* terminal case *)
let terminal =
if !test_mode then
empty_terminal
else
build_terminal_case hyp_list infeasible bb in
Lp_terminal terminal
else
(* split case *)
let n = length (hd bb.std_faces_not_super) in
if n = 3 then
if hard_flag then
split3_hard hyp_list bb
else
split3 hard_flag hyp_list bb
else if n = 4 then
split4 hard_flag hyp_list bb
else if n = 5 then
split5 hard_flag hyp_list bb
else if n = 6 then
split6 hard_flag hyp_list bb
else
failwith (sprintf "build: incorrect face size - %d" n) in
itlist (fun info c -> Lp_split (info, [c])) info_list certificate
(* split3 - hard cases *)
and split3_hard hyp_list bb =
let _ = assert (bb.hints <> []) in
let bbs, split_type, list =
match hd (bb.hints) with
| Triangle_split d -> switch_std3 d bb, "tri", face_of_dart d bb
| Edge_split d -> switch_edge d bb, "edge", d
| High_low i -> switch_node bb i, (if mem i (highish bb) then "236" else "218"), [i] in
let _ = map clear_hint bbs in
let case_args = zip [hyp_list; hyp_list] bbs in
let cases = map (build true) case_args in
let info = {split_type = split_type; split_face = list} in
Lp_split (info, cases)
(* split3 *)
and split3 hard_flag hyp_list bb =
let split_face = hd (std_tri_prebranch bb) in
let _ = assert (length split_face = 3) in
let bbs = switch3 bb in
let case_args = zip [hyp_list; hyp_list] bbs in
let cases = map (build hard_flag) case_args in
let info = {split_type = "tri"; split_face = split_face} in
Lp_split (info, cases)
(* split4 *)
and split4 hard_flag hyp_list bb =
let split_face = hd bb.std_faces_not_super in
let _ = assert (length split_face = 4) in
let bbs = switch4 bb in
let dart13 = nth split_face 1, nth split_face 2 and
dart24 = nth split_face 0, nth split_face 1 in
let split13 = split_list hyp_list dart13 and
split24 = split_list hyp_list dart24 in
let case_args = zip [split13; split24; split13; split24; hyp_list] bbs in
let cases = map (build hard_flag) case_args in
let info = {split_type = "quad"; split_face = split_face} in
Lp_split (info, cases)
(* split5 *)
and split5 hard_flag hyp_list bb =
let split_face = hd bb.std_faces_not_super in
let _ = assert (length split_face = 5) in
let bbs = switch5 bb in
let darts = rotateL 1 (list_pairs split_face) in
let splits_one = map (split_list hyp_list) darts in
let splits_two = map2 split_list splits_one (rotateL 2 darts) in
let case_args = zip (hyp_list :: (splits_one @ splits_two)) bbs in
let cases = map (build hard_flag) case_args in
let info = {split_type = "pent"; split_face = split_face} in
Lp_split (info, cases)
(* split6 *)
and split6 hard_flag hyp_list bb =
let split_face = hd bb.std_faces_not_super in
let _ = assert (length split_face = 6) in
let bbs = switch6 bb in
let darts = Glpk_link.rotateL 1 (list_pairs split_face) in
let splits = map (split_list hyp_list) darts in
let case_args = zip (hyp_list :: splits) bbs in
let cases = map (build hard_flag) case_args in
let info = {split_type = "hex"; split_face = split_face} in
Lp_split (info, cases);;
(* Moves all hex faces to std56_flat_free *)
(* This operation can be done on all hypermaps since branching on hexes is not required: *)
(* there are no inequalities for std56_flat_free INTER std6 *)
let modify_hex_cases bb =
let faces6 = filter (fun f -> length f = 6) bb.std_faces_not_super in
itlist (fun fc bb -> modify_bb bb true ["flat_free", fc] []) faces6 bb;;
(* Builds an lp certificate *)
let build_certificate modify_hex bb =
let _ = reset_build_counters() in
let bb = if modify_hex then modify_hex_cases bb else bb in
let hyp_list = (snd o convert_to_list) bb.string_rep in
let hard_flag = bb.hints <> [] in
let root = build hard_flag (hyp_list, bb) in
let _ = report_build_progress() in
let certificate = {hypermap_string = bb.string_rep; root_case = root} in
certificate;;
(************************)
(* Builds and saves certificates for all given cases *)
let build_and_save_all =
let counter = ref 0 and
total = ref 0 in
let buf = ref [] and
file_counter = ref 0 and
buf_terminals = ref 0 in
let save_buf name_prefix =
if length !buf > 0 then
let _ = file_counter := !file_counter + 1 in
let fname = sprintf "%s/%s_%d.dat" !output_dir name_prefix !file_counter in
let _ = write_lp_certificates fname !buf in
buf := [];
buf_terminals := 0
else
() in
let process name_prefix max modify_hex bb =
let _ = counter := !counter + 1 in
let _ = report (sprintf "%d/%d" !counter !total) in
let certificate = build_certificate modify_hex bb in
let _ = buf := certificate :: !buf in
let _ = buf_terminals := !buf_terminals + count_terminals certificate.root_case in
if !buf_terminals >= max then
save_buf name_prefix
else
() in
fun name_prefix max modify_hex bbs ->
let _ = total := length bbs and
_ = counter := 0 and
_ = file_counter := 0 and
_ = buf := [] and
_ = buf_terminals := 0 in
let _ = map (process name_prefix max modify_hex) bbs in
save_buf name_prefix;;
end;;