1 (* ========================================================================== *)
4 (* Linear Programming, AMPL format (non-formal) *)
7 (* Author: Thomas C. Hales *)
8 (* Date: 2009-09-22, rechecked 2010-06-03, 2012-06-08 *)
9 (* ========================================================================== *)
13 This file generates the informal linear programming part of the cluster inequality (OXLZLEZ).
14 These linear programs reduce the 3 and 4 blade cases to a single case:
15 4 blades with 3 quarters and 1 4-cell of weight 0.5.
16 This final case is handled separately.
20 needs new mktop on platforms that do not support dynamic loading of Str.
22 ocamlmktop unix.cma str.cma -o ocampl
28 flyspeck_needs "../glpk/glpk_link.ml";;
30 module Oxlzlez_informal = struct
36 let sprintf = Printf.sprintf;;
38 let glpk_dir = flyspeck_dir ^ "/../glpk";;
41 let model = glpk_dir^ "/minorlp/OXLZLEZ.mod";;
42 let tmpfile = "/tmp/OXLZLEZ_informal.dat";;
43 let dumpfile = "/tmp/OXLZLEZ_informal.out";;
45 type lptype = Lp_unset
47 | Lp_value of (((string * int list) * string * float) list)*float;;
50 let string_of_lptype t = match t with
51 | Lp_infeasible -> "infeasible"
53 | Lp_value (_,u) -> Printf.sprintf "%3.3f" u;;
56 (* fields of bladerunner are documented in OXLZLEZ.mod.
57 See CBLADE,SBLADE,NONSBLADE,QU,QX,QY,QXD,NONQXD,NEGQU,POSQU,
58 HALFWT,FULLWT,SHORTY4,LONGY4
63 mutable lpvalue : lptype;
64 mutable history : string;
65 cblade : int; (* number of blades *)
81 let next br i = (i+1) mod br.cblade;;
85 (* the initial configuration always has a quarter at 0 *)
89 sblade = [0;1]; (* quarter at 0. "raw" blades: face 0 goes with raw blades 0 & 1. *)
105 let _ = (backup := br ::!backup) in
108 let modify_br h br fields =
109 let _ = (br.history <- h) in
110 let add s vs = nub((get_values s fields) @ vs) in
113 sblade = add "sblade" br.sblade;
114 nonsblade = add "nonsblade" br.nonsblade;
118 qxd = add "qxd" br.qxd;
119 nonqxd = add "nonqxd" br.nonqxd;
120 negqu = add "negqu" br.negqu;
121 posqu = add "posqu" br.posqu;
122 halfwt = add "halfwt" br. halfwt;
123 fullwt = add "fullwt" br.fullwt;
124 shorty4 = add "shorty4" br.shorty4; (* y4 <= 2.1 *)
125 longy4 = add "longy4" br.longy4; (* y4 >= 2.1 *)
129 let _ = (backup:= br':: !backup) in
134 Example: move 1 into halfwt
135 let brx = modify_br "" (mk_br 4) ["halfwt",1];;
138 let ampl_of_br outs br =
139 let list_of = unsplit " " string_of_int in
141 let mk s f = p"set %s := %s;" s (list_of f) in
143 p"param CBLADE := %d;" br.cblade ;
144 mk "SBLADERAW" br.sblade;
145 mk "NONSBLADERAW" br.nonsblade;
150 mk "NONQXD" br.nonqxd;
153 mk "HALFWT" br. halfwt;
154 mk "FULLWT" br.fullwt;
155 mk "LONGY4" br.longy4;
156 mk "SHORTY4" br.shorty4] in
157 Printf.fprintf outs "%s" j;;
161 let br = modify_br "test" br ["qu",1;"qu",2;"negqu",3] in
162 display_ampl tmpfile ampl_of_br br;;
164 (* manipulations of lpvalue field that don't involve glpsol *)
166 let notdone r = (r < 0.0);;
168 let select_notdone f brs = (* selects unset and those satisfying f *)
169 let fil ro = match ro.lpvalue with
171 | Lp_infeasible -> false
172 | Lp_value (_,r) -> f r in
175 let is_unset br = match br.lpvalue with
179 let calc_min brs = fold_right
180 (fun br x -> match br.lpvalue with
181 |Lp_value (_,y) -> min x y
186 let t = calc_min brs in
187 (* let r = Lp_value (t) in *)
189 match br.lpvalue with
190 | Lp_value (_,y) -> (y=t)
196 We fail on the branching if the branch has been made already,
197 or if the necessary prior branches have not been followed.
200 let branch h br ss i =
201 map (fun s -> modify_br h br [s,i]) ss;;
203 let branch_sblade i br =
204 let _ = not(mem i br.sblade or mem i br.nonsblade) or failwith "sblade" in
205 branch "sblade" br ["sblade"; "nonsblade"] i;;
207 let branch_qxd i br =
208 let _ = mem i br.qx or failwith "qxd1" in
209 let _ = not(mem i br.qxd or mem i br.nonqxd) or failwith "qxd2" in
210 branch "qxd" br ["qxd";"nonqxd"] i;;
212 let branch_negqu i br =
213 let _ = mem i br.qu or failwith "negqu1" in
214 let _ = not(mem i br.negqu or mem i br.posqu) or failwith "negqu2" in
215 branch "negqu" br ["negqu";"posqu"] i;;
218 let _ = mem i br.qx or failwith "wt-qx" in
219 let _ = (mem i br.sblade && mem (next br i) br.sblade) or failwith "wt-blade" in
220 let _ = not(mem i br.halfwt or mem i br.fullwt) or failwith "wt-set" in
221 branch "wt" br ["halfwt";"fullwt"] i;;
224 let _ = mem i br.qy or failwith "y4-qy" in
225 let _ = (mem i br.sblade && mem (next br i) br.sblade) or failwith "y4-blade" in
226 branch "y4" br ["shorty4";"longy4"] i;;
230 let _ = not(mem i br.qu or mem i br.qx or mem i br.qy) or failwith "qu-set" in
231 let _ = not(mem i br.nonsblade or mem j br.nonsblade) or failwith "qu-blade" in
232 modify_br "qus" br ["sblade",i;"sblade",j;"qu",i] :: branch "qu" br ["qx";"qy"] i;;
236 let br1 = List.nth (branch_qu 1 br) 1;;
237 let br2 = List.nth (branch_sblade 2 br1) 0;;
241 (* Link in glpsol linear programming package *)
244 let ss = Str.split (Str.regexp "[],[]") s in
245 (hd ss,map int_of_string (tl ss));;
249 let outputf = Flyspeck_lib.load_file "/tmp/output.out" in
250 let output_split = map (Str.split (Str.regexp " +")) outputf in
251 let output_active = filter
253 if (List.length xs <3) then false else
254 let n = List.nth xs 0 in
255 let t = List.nth xs 2 in
256 let v = hd (List.rev xs) in
257 (can int_of_string n && t.[0]='N') && not(v="eps")) output_split in
258 let output_format = map
260 (strip_id (List.nth xs 1),List.nth xs 2,
261 float_of_string (hd (List.rev xs))))
265 let set_lpvalue nt (dualdata,(f,r)) = (* side effects *)
267 if (List.length f > 0) then nt.lpvalue <- Lp_infeasible
268 else if (List.length r = 1)
270 nt.lpvalue <- Lp_value (dualdata, float_of_string(hd r))
272 nt.lpvalue <- Lp_unset in
275 let solve_and_set_lp br =
276 let _ = match br.lpvalue with
277 | Lp_unset -> (* generate dual data as soon as lp is solved *)
278 let fr = (Glpk_link.solve_dual_f model dumpfile "ggsum" ampl_of_br br) in
279 let dualdata = load_dual() in
280 (set_lpvalue br (dualdata,fr))
282 let _ = report (string_of_lptype br.lpvalue) in
288 let ex0 brancher i brs = (flatten (map (brancher i) brs));;
290 let ex brancher i brs =
291 let brs' = ex0 brancher i brs in
292 (* let _ = map set_lpvalue brs' in *)
293 select_notdone notdone (map solve_and_set_lp brs');;
295 let top brancher i (br::rest) = (ex brancher i [br]) @ rest;;
297 let bury (x::xs) = xs @ [x];;
299 (* case of 3 blades: *)
302 let br1 = ex branch_qu 2 (branch_qu 1 br) in
303 let br2 = ex branch_negqu 0 br1 in
304 let br3 = ex branch_qxd 2 br2 in
305 let br4 = ex branch_qxd 1 br3 in
308 (* case of 4 blades *)
311 let cr1 = ex branch_qu 3 (ex0 branch_qu 2 (branch_qu 1 cr)) in
312 let cr2 = bury (top branch_wt 3 cr1) in
313 let cr2' = top branch_y4 3 cr2 in
314 let cr3 = bury (top branch_wt 2 cr2') in
315 let cr4 = top branch_sblade 3 cr3 in
316 let cr5' = top branch_sblade 3 cr4 in
317 let cr5 = top branch_y4 2 cr5' in
318 let cr6 = top branch_sblade 3 cr5 in
319 let cr7 = bury(top branch_wt 1 cr6) in
320 let cr8 = top branch_sblade 2 cr7 in
321 let cr9 = top branch_sblade 2 cr8 in
322 let cr10' = top branch_sblade 2 cr9 in
323 let cr10 = top branch_y4 1 cr10' in
324 let cr11 = top branch_sblade 2 cr10 in
327 (* three cases remain, all related by symmetry. 4 blades, 3 quarters, 1 4-cell with weight 0.5 *)
329 let test_structure br =
330 let _ = (sort (<) br.sblade = [0;1;2;3]) or failwith "sblade" in
331 let _ = ((length br.qx = 1) && (br.qx = br.halfwt)) or failwith "weight" in
332 let _ = (length br.qu = 3) or failwith "qu" in
338 let b3_info = if (List.length b3 = 0) then "blade 3 passes" else failwith "blade 3 fails in OXLZLEZ" in
340 let t0 = nub (map test_structure (b4)) in
341 let b4_info = if (t0=[true]) then "blade 4 passes (ignoring cases with 4 blades, 3 quarters, 1 4-cell with wt 0.5)"
342 else failwith "blade 4 fails in OXLZLEZ" in
343 let s = join_lines [b3_info;b4_info] in