glpk/minorlp/OXLZLEZ.ml

   1 (* ========================================================================== *)
   2 (* FLYSPECK - GLPK                                                            *)
   3 (*                                                                            *)
   4 (* Linear Programming, AMPL format (non-formal)                               *)
   5 (* Chapter: Packing                                                           *)
   6 (* Lemma: OXLZLEZ                                                             *)
   7 (* Author: Thomas C. Hales                                                    *)
   8 (* Date: 2009-09-22, rechecked 2010-06-03, 2012-06-08                         *)
   9 (* ========================================================================== *)
  10
  11 (*
  12
  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.
  17
  18
  19
  20 needs new mktop on platforms that do not support dynamic loading of Str.
  21
  22 ocamlmktop unix.cma str.cma -o ocampl
  23 ./ocampl
  24
  25 *)
  26
  27
  28 flyspeck_needs "../glpk/glpk_link.ml";;
  29
  30 module Oxlzlez_informal = struct
  31
  32 open Str;;
  33 open List;;
  34 open Glpk_link;;
  35
  36 let sprintf = Printf.sprintf;;
  37
  38 let glpk_dir = flyspeck_dir ^ "/../glpk";;
  39
  40 (* external files *)
  41 let model = glpk_dir^ "/minorlp/OXLZLEZ.mod";;
  42 let tmpfile = "/tmp/OXLZLEZ_informal.dat";;
  43 let dumpfile = "/tmp/OXLZLEZ_informal.out";;
  44
  45 type lptype = Lp_unset
  46               | Lp_infeasible
  47               | Lp_value of (((string * int list) * string * float) list)*float;;
  48
  49
  50 let string_of_lptype t = match t with
  51     | Lp_infeasible -> "infeasible"
  52     | Lp_unset -> "unset"
  53     | Lp_value (_,u) -> Printf.sprintf "%3.3f" u;;
  54
  55
  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
  59 *)
  60
  61 type bladerunner =
  62   {
  63     mutable lpvalue : lptype;
  64     mutable history : string;
  65     cblade : int; (* number of blades *)
  66     sblade : int list;
  67     nonsblade : int list;
  68     qu: int list;
  69     qx : int list;
  70     qy : int list;
  71     qxd : int list;
  72     nonqxd : int list;
  73     negqu : int list;
  74     posqu : int list;
  75     halfwt : int list;
  76     fullwt : int list;
  77     shorty4: int list;
  78     longy4:int list;
  79   };;
  80
  81 let next br i = (i+1) mod br.cblade;;
  82
  83 let backup = ref [];;
  84
  85 (* the initial configuration always has a quarter at 0 *)
  86 let mk_br n =
  87   let br =
  88  {cblade = n;
  89   sblade = [0;1];   (* quarter at 0.  "raw" blades: face 0 goes with raw blades 0 & 1. *)
  90   nonsblade = [];
  91   qu = [0];
  92   qx = [];
  93   qy = [];
  94   qxd = [];
  95   nonqxd = [];
  96   negqu = [];
  97   posqu = [];
  98   halfwt = [];
  99   fullwt = [];
 100   shorty4=[];
 101   longy4=[];
 102   lpvalue = Lp_unset;
 103   history="";
 104  } in
 105   let _ = (backup := br ::!backup) in
 106   br;;
 107
 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
 111   let br' = {
 112 cblade = br.cblade;
 113 sblade = add "sblade" br.sblade;
 114 nonsblade = add "nonsblade" br.nonsblade;
 115 qu = add "qu" br.qu;
 116 qx = add "qx" br.qx;
 117 qy = add "qy" br.qy;
 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 *)
 126 lpvalue = Lp_unset;
 127 history = "";
 128 } in
 129   let _ = (backup:= br':: !backup) in
 130   br'
 131 ;;
 132
 133 (*
 134 Example: move 1 into halfwt
 135 let brx = modify_br "" (mk_br 4)   ["halfwt",1];;
 136 *)
 137
 138 let ampl_of_br outs br =
 139   let list_of = unsplit " " string_of_int in
 140   let p = sprintf in
 141   let mk s f = p"set %s := %s;" s (list_of f) in
 142   let j = join_lines [
 143     p"param CBLADE := %d;" br.cblade ;
 144     mk "SBLADERAW" br.sblade;
 145     mk "NONSBLADERAW" br.nonsblade;
 146     mk "QU" br.qu;
 147     mk "QX" br.qx;
 148     mk "QY" br.qy;
 149     mk "QXD" br.qxd;
 150     mk "NONQXD" br.nonqxd;
 151     mk "NEGQU" br.negqu;
 152     mk "POSQU" br.posqu;
 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;;
 158
 159 let test() =
 160   let br = mk_br 4 in
 161   let br =  modify_br "test" br  ["qu",1;"qu",2;"negqu",3] in
 162     display_ampl tmpfile ampl_of_br br;;
 163
 164 (* manipulations of lpvalue field that don't involve glpsol *)
 165
 166 let notdone r = (r < 0.0);;
 167
 168 let select_notdone f brs = (* selects unset and those satisfying f *)
 169     let fil ro = match ro.lpvalue with
 170         Lp_unset -> true
 171       | Lp_infeasible -> false
 172       | Lp_value (_,r) -> f r in
 173       filter fil brs;;
 174
 175 let is_unset br = match br.lpvalue with
 176     Lp_unset -> true
 177   | _ -> false;;
 178
 179 let calc_min brs = fold_right
 180   (fun br x -> match br.lpvalue with
 181      |Lp_value (_,y) -> min x y
 182      |_ -> x
 183   ) brs 10.0;;
 184
 185 let find_min brs =
 186   let t = calc_min brs in
 187   (* let r = Lp_value (t) in *)
 188     find (fun br ->
 189       match br.lpvalue with
 190         |  Lp_value (_,y) -> (y=t)
 191         | _ -> false) brs;;
 192
 193
 194 (*
 195 branching.
 196 We fail on the branching if the branch has been made already,
 197  or if the necessary prior branches have not been followed.
 198 *)
 199
 200 let branch h br ss i  =
 201    map (fun s -> modify_br h br [s,i]) ss;;
 202
 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;;
 206
 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;;
 211
 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;;
 216
 217 let branch_wt i br  =
 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;;
 222
 223 let branch_y4 i br =
 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;;
 227
 228 let branch_qu i br  =
 229   let j = next br i in
 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;;
 233
 234 (* (* example *)
 235 let br = mk_br 3;;
 236 let br1 = List.nth (branch_qu 1 br) 1;;
 237 let br2 = List.nth (branch_sblade 2 br1) 0;;
 238 branch_wt 1 br2;;
 239 *)
 240
 241 (* Link in glpsol linear programming package *)
 242
 243 let strip_id s =
 244   let ss = Str.split (Str.regexp "[],[]") s in
 245     (hd ss,map int_of_string (tl ss));;
 246
 247
 248 let load_dual() =
 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
 252     (fun xs ->
 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
 259     (fun xs ->
 260        (strip_id (List.nth xs 1),List.nth xs 2,
 261         float_of_string (hd (List.rev xs))))
 262     output_active in
 263     output_format;;
 264
 265 let set_lpvalue nt (dualdata,(f,r)) = (* side effects *)
 266   let _ =
 267     if (List.length f > 0) then nt.lpvalue <- Lp_infeasible
 268     else if (List.length r = 1)
 269     then
 270       nt.lpvalue <- Lp_value  (dualdata, float_of_string(hd r))
 271     else
 272       nt.lpvalue <- Lp_unset in
 273     nt;;
 274
 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))
 281     |  _ -> br in
 282   let _ = report (string_of_lptype br.lpvalue) in
 283   br;;
 284
 285
 286 (* control flow *)
 287
 288 let ex0 brancher i brs = (flatten (map (brancher i) brs));;
 289
 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');;
 294
 295 let top brancher i (br::rest) = (ex brancher i [br]) @ rest;;
 296
 297 let bury (x::xs) = xs @ [x];;
 298
 299 (* case of 3 blades: *)
 300 let blade3() =
 301   let br = mk_br 3 in
 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
 306    br4;;
 307
 308 (* case of 4 blades *)
 309 let blade4() =
 310   let cr = mk_br 4 in
 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
 325    cr11;;
 326
 327 (* three cases remain, all related by symmetry. 4 blades, 3 quarters, 1 4-cell with weight 0.5 *)
 328
 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
 333    true;;
 334
 335
 336 let execute() =
 337   let b3 = blade3() in
 338   let b3_info = if (List.length b3 = 0) then "blade 3 passes" else failwith "blade 3 fails in OXLZLEZ" in
 339   let b4 = blade4() 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
 344   let _ = report s in
 345     s;;
 346
 347 end;;