(* ========================================================================== *) (* FLYSPECK - GLPK *) (* *) (* Linear Programming, AMPL format (non-formal) *) (* Chapter: Packing *) (* Lemma: OXLZLEZ *) (* Author: Thomas C. Hales *) (* Date: 2009-09-22, rechecked 2010-06-03, 2012-06-08 *) (* ========================================================================== *) (* This file generates the informal linear programming part of the cluster inequality (OXLZLEZ). These linear programs reduce the 3 and 4 blade cases to a single case: 4 blades with 3 quarters and 1 4-cell of weight 0.5. This final case is handled separately. needs new mktop on platforms that do not support dynamic loading of Str. ocamlmktop unix.cma str.cma -o ocampl ./ocampl *) flyspeck_needs "../glpk/glpk_link.ml";; module Oxlzlez_informal = struct open Str;; open List;; open Glpk_link;; let sprintf = Printf.sprintf;; let glpk_dir = flyspeck_dir ^ "/../glpk";; (* external files *) let model = glpk_dir^ "/minorlp/OXLZLEZ.mod";; let tmpfile = "/tmp/OXLZLEZ_informal.dat";; let dumpfile = "/tmp/OXLZLEZ_informal.out";; type lptype = Lp_unset | Lp_infeasible | Lp_value of (((string * int list) * string * float) list)*float;; let string_of_lptype t = match t with | Lp_infeasible -> "infeasible" | Lp_unset -> "unset" | Lp_value (_,u) -> Printf.sprintf "%3.3f" u;; (* fields of bladerunner are documented in OXLZLEZ.mod. See CBLADE,SBLADE,NONSBLADE,QU,QX,QY,QXD,NONQXD,NEGQU,POSQU, HALFWT,FULLWT,SHORTY4,LONGY4 *) type bladerunner = { mutable lpvalue : lptype; mutable history : string; cblade : int; (* number of blades *) sblade : int list; nonsblade : int list; qu: int list; qx : int list; qy : int list; qxd : int list; nonqxd : int list; negqu : int list; posqu : int list; halfwt : int list; fullwt : int list; shorty4: int list; longy4:int list; };; let next br i = (i+1) mod br.cblade;; let backup = ref [];; (* the initial configuration always has a quarter at 0 *) let mk_br n = let br = {cblade = n; sblade = [0;1]; (* quarter at 0. "raw" blades: face 0 goes with raw blades 0 & 1. *) nonsblade = []; qu = [0]; qx = []; qy = []; qxd = []; nonqxd = []; negqu = []; posqu = []; halfwt = []; fullwt = []; shorty4=[]; longy4=[]; lpvalue = Lp_unset; history=""; } in let _ = (backup := br ::!backup) in br;; let modify_br h br fields = let _ = (br.history <- h) in let add s vs = nub((get_values s fields) @ vs) in let br' = { cblade = br.cblade; sblade = add "sblade" br.sblade; nonsblade = add "nonsblade" br.nonsblade; qu = add "qu" br.qu; qx = add "qx" br.qx; qy = add "qy" br.qy; qxd = add "qxd" br.qxd; nonqxd = add "nonqxd" br.nonqxd; negqu = add "negqu" br.negqu; posqu = add "posqu" br.posqu; halfwt = add "halfwt" br. halfwt; fullwt = add "fullwt" br.fullwt; shorty4 = add "shorty4" br.shorty4; (* y4 <= 2.1 *) longy4 = add "longy4" br.longy4; (* y4 >= 2.1 *) lpvalue = Lp_unset; history = ""; } in let _ = (backup:= br':: !backup) in br' ;; (* Example: move 1 into halfwt let brx = modify_br "" (mk_br 4) ["halfwt",1];; *) let ampl_of_br outs br = let list_of = unsplit " " string_of_int in let p = sprintf in let mk s f = p"set %s := %s;" s (list_of f) in let j = join_lines [ p"param CBLADE := %d;" br.cblade ; mk "SBLADERAW" br.sblade; mk "NONSBLADERAW" br.nonsblade; mk "QU" br.qu; mk "QX" br.qx; mk "QY" br.qy; mk "QXD" br.qxd; mk "NONQXD" br.nonqxd; mk "NEGQU" br.negqu; mk "POSQU" br.posqu; mk "HALFWT" br. halfwt; mk "FULLWT" br.fullwt; mk "LONGY4" br.longy4; mk "SHORTY4" br.shorty4] in Printf.fprintf outs "%s" j;; let test() = let br = mk_br 4 in let br = modify_br "test" br ["qu",1;"qu",2;"negqu",3] in display_ampl tmpfile ampl_of_br br;; (* manipulations of lpvalue field that don't involve glpsol *) let notdone r = (r < 0.0);; let select_notdone f brs = (* selects unset and those satisfying f *) let fil ro = match ro.lpvalue with Lp_unset -> true | Lp_infeasible -> false | Lp_value (_,r) -> f r in filter fil brs;; let is_unset br = match br.lpvalue with Lp_unset -> true | _ -> false;; let calc_min brs = fold_right (fun br x -> match br.lpvalue with |Lp_value (_,y) -> min x y |_ -> x ) brs 10.0;; let find_min brs = let t = calc_min brs in (* let r = Lp_value (t) in *) find (fun br -> match br.lpvalue with | Lp_value (_,y) -> (y=t) | _ -> false) brs;; (* branching. We fail on the branching if the branch has been made already, or if the necessary prior branches have not been followed. *) let branch h br ss i = map (fun s -> modify_br h br [s,i]) ss;; let branch_sblade i br = let _ = not(mem i br.sblade or mem i br.nonsblade) or failwith "sblade" in branch "sblade" br ["sblade"; "nonsblade"] i;; let branch_qxd i br = let _ = mem i br.qx or failwith "qxd1" in let _ = not(mem i br.qxd or mem i br.nonqxd) or failwith "qxd2" in branch "qxd" br ["qxd";"nonqxd"] i;; let branch_negqu i br = let _ = mem i br.qu or failwith "negqu1" in let _ = not(mem i br.negqu or mem i br.posqu) or failwith "negqu2" in branch "negqu" br ["negqu";"posqu"] i;; let branch_wt i br = let _ = mem i br.qx or failwith "wt-qx" in let _ = (mem i br.sblade && mem (next br i) br.sblade) or failwith "wt-blade" in let _ = not(mem i br.halfwt or mem i br.fullwt) or failwith "wt-set" in branch "wt" br ["halfwt";"fullwt"] i;; let branch_y4 i br = let _ = mem i br.qy or failwith "y4-qy" in let _ = (mem i br.sblade && mem (next br i) br.sblade) or failwith "y4-blade" in branch "y4" br ["shorty4";"longy4"] i;; let branch_qu i br = let j = next br i in let _ = not(mem i br.qu or mem i br.qx or mem i br.qy) or failwith "qu-set" in let _ = not(mem i br.nonsblade or mem j br.nonsblade) or failwith "qu-blade" in modify_br "qus" br ["sblade",i;"sblade",j;"qu",i] :: branch "qu" br ["qx";"qy"] i;; (* (* example *) let br = mk_br 3;; let br1 = List.nth (branch_qu 1 br) 1;; let br2 = List.nth (branch_sblade 2 br1) 0;; branch_wt 1 br2;; *) (* Link in glpsol linear programming package *) let strip_id s = let ss = Str.split (Str.regexp "[],[]") s in (hd ss,map int_of_string (tl ss));; let load_dual() = let outputf = Flyspeck_lib.load_file "/tmp/output.out" in let output_split = map (Str.split (Str.regexp " +")) outputf in let output_active = filter (fun xs -> if (List.length xs <3) then false else let n = List.nth xs 0 in let t = List.nth xs 2 in let v = hd (List.rev xs) in (can int_of_string n && t.[0]='N') && not(v="eps")) output_split in let output_format = map (fun xs -> (strip_id (List.nth xs 1),List.nth xs 2, float_of_string (hd (List.rev xs)))) output_active in output_format;; let set_lpvalue nt (dualdata,(f,r)) = (* side effects *) let _ = if (List.length f > 0) then nt.lpvalue <- Lp_infeasible else if (List.length r = 1) then nt.lpvalue <- Lp_value (dualdata, float_of_string(hd r)) else nt.lpvalue <- Lp_unset in nt;; let solve_and_set_lp br = let _ = match br.lpvalue with | Lp_unset -> (* generate dual data as soon as lp is solved *) let fr = (Glpk_link.solve_dual_f model dumpfile "ggsum" ampl_of_br br) in let dualdata = load_dual() in (set_lpvalue br (dualdata,fr)) | _ -> br in let _ = report (string_of_lptype br.lpvalue) in br;; (* control flow *) let ex0 brancher i brs = (flatten (map (brancher i) brs));; let ex brancher i brs = let brs' = ex0 brancher i brs in (* let _ = map set_lpvalue brs' in *) select_notdone notdone (map solve_and_set_lp brs');; let top brancher i (br::rest) = (ex brancher i [br]) @ rest;; let bury (x::xs) = xs @ [x];; (* case of 3 blades: *) let blade3() = let br = mk_br 3 in let br1 = ex branch_qu 2 (branch_qu 1 br) in let br2 = ex branch_negqu 0 br1 in let br3 = ex branch_qxd 2 br2 in let br4 = ex branch_qxd 1 br3 in br4;; (* case of 4 blades *) let blade4() = let cr = mk_br 4 in let cr1 = ex branch_qu 3 (ex0 branch_qu 2 (branch_qu 1 cr)) in let cr2 = bury (top branch_wt 3 cr1) in let cr2' = top branch_y4 3 cr2 in let cr3 = bury (top branch_wt 2 cr2') in let cr4 = top branch_sblade 3 cr3 in let cr5' = top branch_sblade 3 cr4 in let cr5 = top branch_y4 2 cr5' in let cr6 = top branch_sblade 3 cr5 in let cr7 = bury(top branch_wt 1 cr6) in let cr8 = top branch_sblade 2 cr7 in let cr9 = top branch_sblade 2 cr8 in let cr10' = top branch_sblade 2 cr9 in let cr10 = top branch_y4 1 cr10' in let cr11 = top branch_sblade 2 cr10 in cr11;; (* three cases remain, all related by symmetry. 4 blades, 3 quarters, 1 4-cell with weight 0.5 *) let test_structure br = let _ = (sort (<) br.sblade = [0;1;2;3]) or failwith "sblade" in let _ = ((length br.qx = 1) && (br.qx = br.halfwt)) or failwith "weight" in let _ = (length br.qu = 3) or failwith "qu" in true;; let execute() = let b3 = blade3() in let b3_info = if (List.length b3 = 0) then "blade 3 passes" else failwith "blade 3 fails in OXLZLEZ" in let b4 = blade4() in let t0 = nub (map test_structure (b4)) in let b4_info = if (t0=[true]) then "blade 4 passes (ignoring cases with 4 blades, 3 quarters, 1 4-cell with wt 0.5)" else failwith "blade 4 fails in OXLZLEZ" in let s = join_lines [b3_info;b4_info] in let _ = report s in s;; end;;