src/HOL/Matrix_LP/Cplex_tools.ML
author wenzelm
Sun May 12 14:25:16 2013 +0200 (2013-05-12)
changeset 51940 958d439b3013
parent 51930 52fd62618631
child 62505 9e2a65912111
permissions -rw-r--r--
decentralized historic settings;
     1 (*  Title:      HOL/Matrix_LP/Cplex_tools.ML
     2     Author:     Steven Obua
     3 
     4 Relevant Isabelle environment settings:
     5 
     6   # LP_SOLVER is the default solver. It can be changed during runtime via Cplex.set_solver.
     7   # First option: use the commercial cplex solver
     8   #LP_SOLVER=CPLEX
     9   #CPLEX_PATH=cplex
    10   # Second option: use the open source glpk solver
    11   #LP_SOLVER=GLPK
    12   #GLPK_PATH=glpsol
    13 *)
    14 
    15 signature CPLEX =
    16 sig
    17 
    18     datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf
    19                        | cplexNeg of cplexTerm
    20                        | cplexProd of cplexTerm * cplexTerm
    21                        | cplexSum of (cplexTerm list)
    22 
    23     datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
    24 
    25     datatype cplexGoal = cplexMinimize of cplexTerm
    26                | cplexMaximize of cplexTerm
    27 
    28     datatype cplexConstr = cplexConstr of cplexComp *
    29                       (cplexTerm * cplexTerm)
    30 
    31     datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
    32                       * cplexComp * cplexTerm
    33              | cplexBound of cplexTerm * cplexComp * cplexTerm
    34 
    35     datatype cplexProg = cplexProg of string
    36                       * cplexGoal
    37                       * ((string option * cplexConstr)
    38                          list)
    39                       * cplexBounds list
    40 
    41     datatype cplexResult = Unbounded
    42              | Infeasible
    43              | Undefined
    44              | Optimal of string *
    45                       (((* name *) string *
    46                     (* value *) string) list)
    47 
    48     datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK
    49 
    50     exception Load_cplexFile of string
    51     exception Load_cplexResult of string
    52     exception Save_cplexFile of string
    53     exception Execute of string
    54 
    55     val load_cplexFile : string -> cplexProg
    56 
    57     val save_cplexFile : string -> cplexProg -> unit
    58 
    59     val elim_nonfree_bounds : cplexProg -> cplexProg
    60 
    61     val relax_strict_ineqs : cplexProg -> cplexProg
    62 
    63     val is_normed_cplexProg : cplexProg -> bool
    64 
    65     val get_solver : unit -> cplexSolver
    66     val set_solver : cplexSolver -> unit
    67     val solve : cplexProg -> cplexResult
    68 end;
    69 
    70 structure Cplex  : CPLEX =
    71 struct
    72 
    73 datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK
    74 
    75 val cplexsolver = Unsynchronized.ref SOLVER_DEFAULT;
    76 fun get_solver () = !cplexsolver;
    77 fun set_solver s = (cplexsolver := s);
    78 
    79 exception Load_cplexFile of string;
    80 exception Load_cplexResult of string;
    81 exception Save_cplexFile of string;
    82 
    83 datatype cplexTerm = cplexVar of string
    84            | cplexNum of string
    85            | cplexInf
    86                    | cplexNeg of cplexTerm
    87                    | cplexProd of cplexTerm * cplexTerm
    88                    | cplexSum of (cplexTerm list)
    89 datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
    90 datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm
    91 datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm)
    92 datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
    93                       * cplexComp * cplexTerm
    94                      | cplexBound of cplexTerm * cplexComp * cplexTerm
    95 datatype cplexProg = cplexProg of string
    96                   * cplexGoal
    97                   * ((string option * cplexConstr) list)
    98                   * cplexBounds list
    99 
   100 fun rev_cmp cplexLe = cplexGe
   101   | rev_cmp cplexLeq = cplexGeq
   102   | rev_cmp cplexGe = cplexLe
   103   | rev_cmp cplexGeq = cplexLeq
   104   | rev_cmp cplexEq = cplexEq
   105 
   106 fun the NONE = raise (Load_cplexFile "SOME expected")
   107   | the (SOME x) = x;
   108 
   109 fun modulo_signed is_something (cplexNeg u) = is_something u
   110   | modulo_signed is_something u = is_something u
   111 
   112 fun is_Num (cplexNum _) = true
   113   | is_Num _ = false
   114 
   115 fun is_Inf cplexInf = true
   116   | is_Inf _ = false
   117 
   118 fun is_Var (cplexVar _) = true
   119   | is_Var _ = false
   120 
   121 fun is_Neg (cplexNeg _) = true
   122   | is_Neg _ = false
   123 
   124 fun is_normed_Prod (cplexProd (t1, t2)) =
   125     (is_Num t1) andalso (is_Var t2)
   126   | is_normed_Prod x = is_Var x
   127 
   128 fun is_normed_Sum (cplexSum ts) =
   129     (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts
   130   | is_normed_Sum x = modulo_signed is_normed_Prod x
   131 
   132 fun is_normed_Constr (cplexConstr (_, (t1, t2))) =
   133     (is_normed_Sum t1) andalso (modulo_signed is_Num t2)
   134 
   135 fun is_Num_or_Inf x = is_Inf x orelse is_Num x
   136 
   137 fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) =
   138     (c1 = cplexLe orelse c1 = cplexLeq) andalso
   139     (c2 = cplexLe orelse c2 = cplexLeq) andalso
   140     is_Var t2 andalso
   141     modulo_signed is_Num_or_Inf t1 andalso
   142     modulo_signed is_Num_or_Inf t3
   143   | is_normed_Bounds (cplexBound (t1, c, t2)) =
   144     (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2))
   145     orelse
   146     (c <> cplexEq andalso
   147      is_Var t2 andalso (modulo_signed is_Num_or_Inf t1))
   148 
   149 fun term_of_goal (cplexMinimize x) = x
   150   | term_of_goal (cplexMaximize x) = x
   151 
   152 fun is_normed_cplexProg (cplexProg (_, goal, constraints, bounds)) =
   153     is_normed_Sum (term_of_goal goal) andalso
   154     forall (fn (_,x) => is_normed_Constr x) constraints andalso
   155     forall is_normed_Bounds bounds
   156 
   157 fun is_NL s = s = "\n"
   158 
   159 fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s)
   160 
   161 fun is_num a =
   162     let
   163     val b = String.explode a
   164     fun num4 cs = forall Char.isDigit cs
   165     fun num3 [] = true
   166       | num3 (ds as (c::cs)) =
   167         if c = #"+" orelse c = #"-" then
   168         num4 cs
   169         else
   170         num4 ds
   171     fun num2 [] = true
   172       | num2 (c::cs) =
   173         if c = #"e" orelse c = #"E" then num3 cs
   174         else (Char.isDigit c) andalso num2 cs
   175     fun num1 [] = true
   176       | num1 (c::cs) =
   177         if c = #"." then num2 cs
   178         else if c = #"e" orelse c = #"E" then num3 cs
   179         else (Char.isDigit c) andalso num1 cs
   180     fun num [] = true
   181       | num (c::cs) =
   182         if c = #"." then num2 cs
   183         else (Char.isDigit c) andalso num1 cs
   184     in
   185     num b
   186     end
   187 
   188 fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":"
   189 
   190 fun is_cmp s = s = "<" orelse s = ">" orelse s = "<="
   191              orelse s = ">=" orelse s = "="
   192 
   193 fun is_symbol a =
   194     let
   195     val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~"
   196     fun is_symbol_char c = Char.isAlphaNum c orelse
   197                    exists (fn d => d=c) symbol_char
   198     fun is_symbol_start c = is_symbol_char c andalso
   199                 not (Char.isDigit c) andalso
   200                 not (c= #".")
   201     val b = String.explode a
   202     in
   203     b <> [] andalso is_symbol_start (hd b) andalso
   204     forall is_symbol_char b
   205     end
   206 
   207 fun to_upper s = String.implode (map Char.toUpper (String.explode s))
   208 
   209 fun keyword x =
   210     let
   211     val a = to_upper x
   212     in
   213     if a = "BOUNDS" orelse a = "BOUND" then
   214         SOME "BOUNDS"
   215     else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then
   216         SOME "MINIMIZE"
   217     else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then
   218         SOME "MAXIMIZE"
   219     else if a = "ST" orelse a = "S.T." orelse a = "ST." then
   220         SOME "ST"
   221     else if a = "FREE" orelse a = "END" then
   222         SOME a
   223     else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then
   224         SOME "GENERAL"
   225     else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then
   226         SOME "INTEGER"
   227     else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then
   228         SOME "BINARY"
   229     else if a = "INF" orelse a = "INFINITY" then
   230         SOME "INF"
   231     else
   232         NONE
   233     end
   234 
   235 val TOKEN_ERROR = ~1
   236 val TOKEN_BLANK = 0
   237 val TOKEN_NUM = 1
   238 val TOKEN_DELIMITER = 2
   239 val TOKEN_SYMBOL = 3
   240 val TOKEN_LABEL = 4
   241 val TOKEN_CMP = 5
   242 val TOKEN_KEYWORD = 6
   243 val TOKEN_NL = 7
   244 
   245 (* tokenize takes a list of chars as argument and returns a list of
   246    int * string pairs, each string representing a "cplex token",
   247    and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP
   248    or TOKEN_SYMBOL *)
   249 fun tokenize s =
   250     let
   251     val flist = [(is_NL, TOKEN_NL),
   252              (is_blank, TOKEN_BLANK),
   253              (is_num, TOKEN_NUM),
   254                      (is_delimiter, TOKEN_DELIMITER),
   255              (is_cmp, TOKEN_CMP),
   256              (is_symbol, TOKEN_SYMBOL)]
   257     fun match_helper [] s = (fn _ => false, TOKEN_ERROR)
   258       | match_helper (f::fs) s =
   259         if ((fst f) s) then f else match_helper fs s
   260     fun match s = match_helper flist s
   261     fun tok s =
   262         if s = "" then [] else
   263         let
   264         val h = String.substring (s,0,1)
   265         val (f, j) = match h
   266         fun len i =
   267             if size s = i then i
   268             else if f (String.substring (s,0,i+1)) then
   269             len (i+1)
   270             else i
   271         in
   272         if j < 0 then
   273             (if h = "\\" then []
   274              else raise (Load_cplexFile ("token expected, found: "
   275                          ^s)))
   276         else
   277             let
   278             val l = len 1
   279             val u = String.substring (s,0,l)
   280             val v = String.extract (s,l,NONE)
   281             in
   282             if j = 0 then tok v else (j, u) :: tok v
   283             end
   284         end
   285     in
   286     tok s
   287     end
   288 
   289 exception Tokenize of string;
   290 
   291 fun tokenize_general flist s =
   292     let
   293     fun match_helper [] s = raise (Tokenize s)
   294       | match_helper (f::fs) s =
   295         if ((fst f) s) then f else match_helper fs s
   296     fun match s = match_helper flist s
   297     fun tok s =
   298         if s = "" then [] else
   299         let
   300         val h = String.substring (s,0,1)
   301         val (f, j) = match h
   302         fun len i =
   303             if size s = i then i
   304             else if f (String.substring (s,0,i+1)) then
   305             len (i+1)
   306             else i
   307         val l = len 1
   308         in
   309         (j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE))
   310         end
   311     in
   312     tok s
   313     end
   314 
   315 fun load_cplexFile name =
   316     let
   317     val f = TextIO.openIn name
   318         val ignore_NL = Unsynchronized.ref true
   319     val rest = Unsynchronized.ref []
   320 
   321     fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s
   322 
   323     fun readToken_helper () =
   324         if length (!rest) > 0 then
   325         let val u = hd (!rest) in
   326             (
   327              rest := tl (!rest);
   328              SOME u
   329             )
   330         end
   331         else
   332           (case TextIO.inputLine f of
   333             NONE => NONE
   334           | SOME s =>
   335             let val t = tokenize s in
   336             if (length t >= 2 andalso
   337                 snd(hd (tl t)) = ":")
   338             then
   339                 rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t))
   340             else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t))
   341                 andalso is_symbol "TO" (hd (tl t))
   342             then
   343                 rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t))
   344             else
   345                 rest := t;
   346             readToken_helper ()
   347             end)
   348 
   349     fun readToken_helper2 () =
   350         let val c = readToken_helper () in
   351             if c = NONE then NONE
   352                     else if !ignore_NL andalso fst (the c) = TOKEN_NL then
   353             readToken_helper2 ()
   354             else if fst (the c) = TOKEN_SYMBOL
   355                 andalso keyword (snd (the c)) <> NONE
   356             then SOME (TOKEN_KEYWORD, the (keyword (snd (the c))))
   357             else c
   358         end
   359 
   360     fun readToken () = readToken_helper2 ()
   361 
   362     fun pushToken a = rest := (a::(!rest))
   363 
   364     fun is_value token =
   365         fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD
   366                       andalso snd token = "INF")
   367 
   368         fun get_value token =
   369         if fst token = TOKEN_NUM then
   370         cplexNum (snd token)
   371         else if fst token = TOKEN_KEYWORD andalso snd token = "INF"
   372         then
   373         cplexInf
   374         else
   375         raise (Load_cplexFile "num expected")
   376 
   377     fun readTerm_Product only_num =
   378         let val c = readToken () in
   379         if c = NONE then NONE
   380         else if fst (the c) = TOKEN_SYMBOL
   381         then (
   382             if only_num then (pushToken (the c); NONE)
   383             else SOME (cplexVar (snd (the c)))
   384             )
   385         else if only_num andalso is_value (the c) then
   386             SOME (get_value (the c))
   387         else if is_value (the c) then
   388             let val t1 = get_value (the c)
   389             val d = readToken ()
   390             in
   391             if d = NONE then SOME t1
   392             else if fst (the d) = TOKEN_SYMBOL then
   393                 SOME (cplexProd (t1, cplexVar (snd (the d))))
   394             else
   395                 (pushToken (the d); SOME t1)
   396             end
   397         else (pushToken (the c); NONE)
   398         end
   399 
   400     fun readTerm_Signed only_signed only_num =
   401         let
   402         val c = readToken ()
   403         in
   404         if c = NONE then NONE
   405         else
   406             let val d = the c in
   407             if d = (TOKEN_DELIMITER, "+") then
   408                 readTerm_Product only_num
   409              else if d = (TOKEN_DELIMITER, "-") then
   410                  SOME (cplexNeg (the (readTerm_Product
   411                               only_num)))
   412              else (pushToken d;
   413                    if only_signed then NONE
   414                    else readTerm_Product only_num)
   415             end
   416         end
   417 
   418     fun readTerm_Sum first_signed =
   419         let val c = readTerm_Signed first_signed false in
   420         if c = NONE then [] else (the c)::(readTerm_Sum true)
   421         end
   422 
   423     fun readTerm () =
   424         let val c = readTerm_Sum false in
   425         if c = [] then NONE
   426         else if tl c = [] then SOME (hd c)
   427         else SOME (cplexSum c)
   428         end
   429 
   430     fun readLabeledTerm () =
   431         let val c = readToken () in
   432         if c = NONE then (NONE, NONE)
   433         else if fst (the c) = TOKEN_LABEL then
   434             let val t = readTerm () in
   435             if t = NONE then
   436                 raise (Load_cplexFile ("term after label "^
   437                            (snd (the c))^
   438                            " expected"))
   439             else (SOME (snd (the c)), t)
   440             end
   441         else (pushToken (the c); (NONE, readTerm ()))
   442         end
   443 
   444     fun readGoal () =
   445         let
   446         val g = readToken ()
   447         in
   448             if g = SOME (TOKEN_KEYWORD, "MAXIMIZE") then
   449             cplexMaximize (the (snd (readLabeledTerm ())))
   450         else if g = SOME (TOKEN_KEYWORD, "MINIMIZE") then
   451             cplexMinimize (the (snd (readLabeledTerm ())))
   452         else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected")
   453         end
   454 
   455     fun str2cmp b =
   456         (case b of
   457          "<" => cplexLe
   458            | "<=" => cplexLeq
   459            | ">" => cplexGe
   460            | ">=" => cplexGeq
   461                | "=" => cplexEq
   462            | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP")))
   463 
   464     fun readConstraint () =
   465             let
   466         val t = readLabeledTerm ()
   467         fun make_constraint b t1 t2 =
   468                     cplexConstr
   469             (str2cmp b,
   470              (t1, t2))
   471         in
   472         if snd t = NONE then NONE
   473         else
   474             let val c = readToken () in
   475             if c = NONE orelse fst (the c) <> TOKEN_CMP
   476             then raise (Load_cplexFile "TOKEN_CMP expected")
   477             else
   478                 let val n = readTerm_Signed false true in
   479                 if n = NONE then
   480                     raise (Load_cplexFile "num expected")
   481                 else
   482                     SOME (fst t,
   483                       make_constraint (snd (the c))
   484                               (the (snd t))
   485                               (the n))
   486                 end
   487             end
   488         end
   489 
   490         fun readST () =
   491         let
   492         fun readbody () =
   493             let val t = readConstraint () in
   494             if t = NONE then []
   495             else if (is_normed_Constr (snd (the t))) then
   496                 (the t)::(readbody ())
   497             else if (fst (the t) <> NONE) then
   498                 raise (Load_cplexFile
   499                        ("constraint '"^(the (fst (the t)))^
   500                     "'is not normed"))
   501             else
   502                 raise (Load_cplexFile
   503                        "constraint is not normed")
   504             end
   505         in
   506         if readToken () = SOME (TOKEN_KEYWORD, "ST")
   507         then
   508             readbody ()
   509         else
   510             raise (Load_cplexFile "ST expected")
   511         end
   512 
   513     fun readCmp () =
   514         let val c = readToken () in
   515         if c = NONE then NONE
   516         else if fst (the c) = TOKEN_CMP then
   517             SOME (str2cmp (snd (the c)))
   518         else (pushToken (the c); NONE)
   519         end
   520 
   521     fun skip_NL () =
   522         let val c = readToken () in
   523         if c <> NONE andalso fst (the c) = TOKEN_NL then
   524             skip_NL ()
   525         else
   526             (pushToken (the c); ())
   527         end
   528 
   529     fun make_bounds c t1 t2 =
   530         cplexBound (t1, c, t2)
   531 
   532     fun readBound () =
   533         let
   534         val _ = skip_NL ()
   535         val t1 = readTerm ()
   536         in
   537         if t1 = NONE then NONE
   538         else
   539             let
   540             val c1 = readCmp ()
   541             in
   542             if c1 = NONE then
   543                 let
   544                 val c = readToken ()
   545                 in
   546                 if c = SOME (TOKEN_KEYWORD, "FREE") then
   547                     SOME (
   548                     cplexBounds (cplexNeg cplexInf,
   549                          cplexLeq,
   550                          the t1,
   551                          cplexLeq,
   552                          cplexInf))
   553                 else
   554                     raise (Load_cplexFile "FREE expected")
   555                 end
   556             else
   557                 let
   558                 val t2 = readTerm ()
   559                 in
   560                 if t2 = NONE then
   561                     raise (Load_cplexFile "term expected")
   562                 else
   563                     let val c2 = readCmp () in
   564                     if c2 = NONE then
   565                         SOME (make_bounds (the c1)
   566                                   (the t1)
   567                                   (the t2))
   568                     else
   569                         SOME (
   570                         cplexBounds (the t1,
   571                              the c1,
   572                              the t2,
   573                              the c2,
   574                              the (readTerm())))
   575                     end
   576                 end
   577             end
   578         end
   579 
   580     fun readBounds () =
   581         let
   582         fun makestring _ = "?"
   583         fun readbody () =
   584             let
   585             val b = readBound ()
   586             in
   587             if b = NONE then []
   588             else if (is_normed_Bounds (the b)) then
   589                 (the b)::(readbody())
   590             else (
   591                 raise (Load_cplexFile
   592                        ("bounds are not normed in: "^
   593                     (makestring (the b)))))
   594             end
   595         in
   596         if readToken () = SOME (TOKEN_KEYWORD, "BOUNDS") then
   597             readbody ()
   598         else raise (Load_cplexFile "BOUNDS expected")
   599         end
   600 
   601         fun readEnd () =
   602         if readToken () = SOME (TOKEN_KEYWORD, "END") then ()
   603         else raise (Load_cplexFile "END expected")
   604 
   605     val result_Goal = readGoal ()
   606     val result_ST = readST ()
   607     val _ =    ignore_NL := false
   608         val result_Bounds = readBounds ()
   609         val _ = ignore_NL := true
   610         val _ = readEnd ()
   611     val _ = TextIO.closeIn f
   612     in
   613     cplexProg (name, result_Goal, result_ST, result_Bounds)
   614     end
   615 
   616 fun save_cplexFile filename (cplexProg (_, goal, constraints, bounds)) =
   617     let
   618     val f = TextIO.openOut filename
   619 
   620     fun basic_write s = TextIO.output(f, s)
   621 
   622     val linebuf = Unsynchronized.ref ""
   623     fun buf_flushline s =
   624         (basic_write (!linebuf);
   625          basic_write "\n";
   626          linebuf := s)
   627     fun buf_add s = linebuf := (!linebuf) ^ s
   628 
   629     fun write s =
   630         if (String.size s) + (String.size (!linebuf)) >= 250 then
   631         buf_flushline ("    "^s)
   632         else
   633         buf_add s
   634 
   635         fun writeln s = (buf_add s; buf_flushline "")
   636 
   637     fun write_term (cplexVar x) = write x
   638       | write_term (cplexNum x) = write x
   639       | write_term cplexInf = write "inf"
   640       | write_term (cplexProd (cplexNum "1", b)) = write_term b
   641       | write_term (cplexProd (a, b)) =
   642         (write_term a; write " "; write_term b)
   643           | write_term (cplexNeg x) = (write " - "; write_term x)
   644           | write_term (cplexSum ts) = write_terms ts
   645     and write_terms [] = ()
   646       | write_terms (t::ts) =
   647         (if (not (is_Neg t)) then write " + " else ();
   648          write_term t; write_terms ts)
   649 
   650     fun write_goal (cplexMaximize term) =
   651         (writeln "MAXIMIZE"; write_term term; writeln "")
   652       | write_goal (cplexMinimize term) =
   653         (writeln "MINIMIZE"; write_term term; writeln "")
   654 
   655     fun write_cmp cplexLe = write "<"
   656       | write_cmp cplexLeq = write "<="
   657       | write_cmp cplexEq = write "="
   658       | write_cmp cplexGe = write ">"
   659       | write_cmp cplexGeq = write ">="
   660 
   661     fun write_constr (cplexConstr (cmp, (a,b))) =
   662         (write_term a;
   663          write " ";
   664          write_cmp cmp;
   665          write " ";
   666          write_term b)
   667 
   668     fun write_constraints [] = ()
   669       | write_constraints (c::cs) =
   670         (if (fst c <> NONE)
   671          then
   672          (write (the (fst c)); write ": ")
   673          else
   674          ();
   675          write_constr (snd c);
   676          writeln "";
   677          write_constraints cs)
   678 
   679     fun write_bounds [] = ()
   680       | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) =
   681         ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf
   682          andalso (c1 = cplexLeq orelse c1 = cplexLe)
   683          andalso (c2 = cplexLeq orelse c2 = cplexLe)
   684           then
   685           (write_term t2; write " free")
   686           else
   687           (write_term t1; write " "; write_cmp c1; write " ";
   688            write_term t2; write " "; write_cmp c2; write " ";
   689            write_term t3)
   690          ); writeln ""; write_bounds bs)
   691       | write_bounds ((cplexBound (t1, c, t2)) :: bs) =
   692         (write_term t1; write " ";
   693          write_cmp c; write " ";
   694          write_term t2; writeln ""; write_bounds bs)
   695 
   696     val _ = write_goal goal
   697         val _ = (writeln ""; writeln "ST")
   698     val _ = write_constraints constraints
   699         val _ = (writeln ""; writeln "BOUNDS")
   700     val _ = write_bounds bounds
   701         val _ = (writeln ""; writeln "END")
   702         val _ = TextIO.closeOut f
   703     in
   704     ()
   705     end
   706 
   707 fun norm_Constr (constr as cplexConstr (c, (t1, t2))) =
   708     if not (modulo_signed is_Num t2) andalso
   709        modulo_signed is_Num t1
   710     then
   711     [cplexConstr (rev_cmp c, (t2, t1))]
   712     else if (c = cplexLe orelse c = cplexLeq) andalso
   713         (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf)
   714     then
   715     []
   716     else if (c = cplexGe orelse c = cplexGeq) andalso
   717         (t1 = cplexInf orelse t2 = cplexNeg cplexInf)
   718     then
   719     []
   720     else
   721     [constr]
   722 
   723 fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) =
   724     (norm_Constr(cplexConstr (c1, (t1, t2))))
   725     @ (norm_Constr(cplexConstr (c2, (t2, t3))))
   726   | bound2constr (cplexBound (t1, cplexEq, t2)) =
   727     (norm_Constr(cplexConstr (cplexLeq, (t1, t2))))
   728     @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1))))
   729   | bound2constr (cplexBound (t1, c1, t2)) =
   730     norm_Constr(cplexConstr (c1, (t1,t2)))
   731 
   732 val emptyset = Symtab.empty
   733 
   734 fun singleton v = Symtab.update (v, ()) emptyset
   735 
   736 fun merge a b = Symtab.merge (op =) (a, b)
   737 
   738 fun mergemap f ts = fold (fn x => fn table => merge table (f x)) ts Symtab.empty
   739 
   740 fun diff a b = Symtab.fold (Symtab.delete_safe o fst) b a
   741 
   742 fun collect_vars (cplexVar v) = singleton v
   743   | collect_vars (cplexNeg t) = collect_vars t
   744   | collect_vars (cplexProd (t1, t2)) =
   745     merge (collect_vars t1) (collect_vars t2)
   746   | collect_vars (cplexSum ts) = mergemap collect_vars ts
   747   | collect_vars _ = emptyset
   748 
   749 (* Eliminates all nonfree bounds from the linear program and produces an
   750    equivalent program with only free bounds
   751    IF for the input program P holds: is_normed_cplexProg P *)
   752 fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) =
   753     let
   754     fun collect_constr_vars (_, cplexConstr (_, (t1,_))) =
   755         (collect_vars t1)
   756 
   757     val cvars = merge (collect_vars (term_of_goal goal))
   758               (mergemap collect_constr_vars constraints)
   759 
   760     fun collect_lower_bounded_vars
   761         (cplexBounds (_, _, cplexVar v, _, _)) =
   762         singleton v
   763       |  collect_lower_bounded_vars
   764          (cplexBound (_, cplexLe, cplexVar v)) =
   765          singleton v
   766       |  collect_lower_bounded_vars
   767          (cplexBound (_, cplexLeq, cplexVar v)) =
   768          singleton v
   769       |  collect_lower_bounded_vars
   770          (cplexBound (cplexVar v, cplexGe,_)) =
   771          singleton v
   772       |  collect_lower_bounded_vars
   773          (cplexBound (cplexVar v, cplexGeq, _)) =
   774          singleton v
   775       | collect_lower_bounded_vars
   776         (cplexBound (cplexVar v, cplexEq, _)) =
   777         singleton v
   778       |  collect_lower_bounded_vars _ = emptyset
   779 
   780     val lvars = mergemap collect_lower_bounded_vars bounds
   781     val positive_vars = diff cvars lvars
   782     val zero = cplexNum "0"
   783 
   784     fun make_pos_constr v =
   785         (NONE, cplexConstr (cplexGeq, ((cplexVar v), zero)))
   786 
   787     fun make_free_bound v =
   788         cplexBounds (cplexNeg cplexInf, cplexLeq,
   789              cplexVar v, cplexLeq,
   790              cplexInf)
   791 
   792     val pos_constrs = rev (Symtab.fold
   793                   (fn (k, _) => cons (make_pos_constr k))
   794                   positive_vars [])
   795         val bound_constrs = map (pair NONE)
   796                 (maps bound2constr bounds)
   797     val constraints' = constraints @ pos_constrs @ bound_constrs
   798     val bounds' = rev (Symtab.fold (fn (v, _) => cons (make_free_bound v)) cvars []);
   799     in
   800     cplexProg (name, goal, constraints', bounds')
   801     end
   802 
   803 fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) =
   804     let
   805     fun relax cplexLe = cplexLeq
   806       | relax cplexGe = cplexGeq
   807       | relax x = x
   808 
   809     fun relax_constr (n, cplexConstr(c, (t1, t2))) =
   810         (n, cplexConstr(relax c, (t1, t2)))
   811 
   812     fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) =
   813         cplexBounds (t1, relax c1, t2, relax c2, t3)
   814       | relax_bounds (cplexBound (t1, c, t2)) =
   815         cplexBound (t1, relax c, t2)
   816     in
   817     cplexProg (name,
   818            goals,
   819            map relax_constr constrs,
   820            map relax_bounds bounds)
   821     end
   822 
   823 datatype cplexResult = Unbounded
   824              | Infeasible
   825              | Undefined
   826              | Optimal of string * ((string * string) list)
   827 
   828 fun is_separator x = forall (fn c => c = #"-") (String.explode x)
   829 
   830 fun is_sign x = (x = "+" orelse x = "-")
   831 
   832 fun is_colon x = (x = ":")
   833 
   834 fun is_resultsymbol a =
   835     let
   836     val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-"
   837     fun is_symbol_char c = Char.isAlphaNum c orelse
   838                    exists (fn d => d=c) symbol_char
   839     fun is_symbol_start c = is_symbol_char c andalso
   840                 not (Char.isDigit c) andalso
   841                 not (c= #".") andalso
   842                 not (c= #"-")
   843     val b = String.explode a
   844     in
   845     b <> [] andalso is_symbol_start (hd b) andalso
   846     forall is_symbol_char b
   847     end
   848 
   849 val TOKEN_SIGN = 100
   850 val TOKEN_COLON = 101
   851 val TOKEN_SEPARATOR = 102
   852 
   853 fun load_glpkResult name =
   854     let
   855     val flist = [(is_NL, TOKEN_NL),
   856              (is_blank, TOKEN_BLANK),
   857              (is_num, TOKEN_NUM),
   858              (is_sign, TOKEN_SIGN),
   859                      (is_colon, TOKEN_COLON),
   860              (is_cmp, TOKEN_CMP),
   861              (is_resultsymbol, TOKEN_SYMBOL),
   862              (is_separator, TOKEN_SEPARATOR)]
   863 
   864     val tokenize = tokenize_general flist
   865 
   866     val f = TextIO.openIn name
   867 
   868     val rest = Unsynchronized.ref []
   869 
   870     fun readToken_helper () =
   871         if length (!rest) > 0 then
   872         let val u = hd (!rest) in
   873             (
   874              rest := tl (!rest);
   875              SOME u
   876             )
   877         end
   878         else
   879         (case TextIO.inputLine f of
   880           NONE => NONE
   881         | SOME s => (rest := tokenize s; readToken_helper()))
   882 
   883     fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty)
   884 
   885     fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest)))
   886 
   887     fun readToken () =
   888         let val t = readToken_helper () in
   889         if is_tt t TOKEN_BLANK then
   890             readToken ()
   891         else if is_tt t TOKEN_NL then
   892             let val t2 = readToken_helper () in
   893             if is_tt t2 TOKEN_SIGN then
   894                 (pushToken (SOME (TOKEN_SEPARATOR, "-")); t)
   895             else
   896                 (pushToken t2; t)
   897             end
   898         else if is_tt t TOKEN_SIGN then
   899             let val t2 = readToken_helper () in
   900             if is_tt t2 TOKEN_NUM then
   901                 (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2))))
   902             else
   903                 (pushToken t2; t)
   904             end
   905         else
   906             t
   907         end
   908 
   909         fun readRestOfLine P =
   910         let
   911         val t = readToken ()
   912         in
   913         if is_tt t TOKEN_NL orelse t = NONE
   914         then P
   915         else readRestOfLine P
   916         end
   917 
   918     fun readHeader () =
   919         let
   920         fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
   921         fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ())))
   922         val t1 = readToken ()
   923         val t2 = readToken ()
   924         in
   925         if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON
   926         then
   927             case to_upper (snd (the t1)) of
   928             "STATUS" => (readStatus ())::(readHeader ())
   929               | "OBJECTIVE" => (readObjective())::(readHeader ())
   930               | _ => (readRestOfLine (); readHeader ())
   931         else
   932             (pushToken t2; pushToken t1; [])
   933         end
   934 
   935     fun skip_until_sep () =
   936         let val x = readToken () in
   937         if is_tt x TOKEN_SEPARATOR then
   938             readRestOfLine ()
   939         else
   940             skip_until_sep ()
   941         end
   942 
   943     fun load_value () =
   944         let
   945         val t1 = readToken ()
   946         val t2 = readToken ()
   947         in
   948         if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then
   949             let
   950             val t = readToken ()
   951             val state = if is_tt t TOKEN_NL then readToken () else t
   952             val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected")
   953             val k = readToken ()
   954             in
   955             if is_tt k TOKEN_NUM then
   956                 readRestOfLine (SOME (snd (the t2), snd (the k)))
   957             else
   958                 raise (Load_cplexResult "number expected")
   959             end
   960         else
   961             (pushToken t2; pushToken t1; NONE)
   962         end
   963 
   964     fun load_values () =
   965         let val v = load_value () in
   966         if v = NONE then [] else (the v)::(load_values ())
   967         end
   968 
   969     val header = readHeader ()
   970 
   971     val result =
   972         case AList.lookup (op =) header "STATUS" of
   973         SOME "INFEASIBLE" => Infeasible
   974           | SOME "UNBOUNDED" => Unbounded
   975           | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) header "OBJECTIVE"),
   976                        (skip_until_sep ();
   977                         skip_until_sep ();
   978                         load_values ()))
   979           | _ => Undefined
   980 
   981     val _ = TextIO.closeIn f
   982     in
   983     result
   984     end
   985     handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
   986      | Option.Option => raise (Load_cplexResult "Option")
   987 
   988 fun load_cplexResult name =
   989     let
   990     val flist = [(is_NL, TOKEN_NL),
   991              (is_blank, TOKEN_BLANK),
   992              (is_num, TOKEN_NUM),
   993              (is_sign, TOKEN_SIGN),
   994                      (is_colon, TOKEN_COLON),
   995              (is_cmp, TOKEN_CMP),
   996              (is_resultsymbol, TOKEN_SYMBOL)]
   997 
   998     val tokenize = tokenize_general flist
   999 
  1000     val f = TextIO.openIn name
  1001 
  1002     val rest = Unsynchronized.ref []
  1003 
  1004     fun readToken_helper () =
  1005         if length (!rest) > 0 then
  1006         let val u = hd (!rest) in
  1007             (
  1008              rest := tl (!rest);
  1009              SOME u
  1010             )
  1011         end
  1012         else
  1013         (case TextIO.inputLine f of
  1014           NONE => NONE
  1015         | SOME s => (rest := tokenize s; readToken_helper()))
  1016 
  1017     fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty)
  1018 
  1019     fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest)))
  1020 
  1021     fun readToken () =
  1022         let val t = readToken_helper () in
  1023         if is_tt t TOKEN_BLANK then
  1024             readToken ()
  1025         else if is_tt t TOKEN_SIGN then
  1026             let val t2 = readToken_helper () in
  1027             if is_tt t2 TOKEN_NUM then
  1028                 (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2))))
  1029             else
  1030                 (pushToken t2; t)
  1031             end
  1032         else
  1033             t
  1034         end
  1035 
  1036         fun readRestOfLine P =
  1037         let
  1038         val t = readToken ()
  1039         in
  1040         if is_tt t TOKEN_NL orelse t = NONE
  1041         then P
  1042         else readRestOfLine P
  1043         end
  1044 
  1045     fun readHeader () =
  1046         let
  1047         fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
  1048         fun readObjective () =
  1049             let
  1050             val t = readToken ()
  1051             in
  1052             if is_tt t TOKEN_SYMBOL andalso to_upper (snd (the t)) = "VALUE" then
  1053                 readRestOfLine ("OBJECTIVE", snd (the (readToken())))
  1054             else
  1055                 readRestOfLine ("OBJECTIVE_NAME", snd (the t))
  1056             end
  1057 
  1058         val t = readToken ()
  1059         in
  1060         if is_tt t TOKEN_SYMBOL then
  1061             case to_upper (snd (the t)) of
  1062             "STATUS" => (readStatus ())::(readHeader ())
  1063               | "OBJECTIVE" => (readObjective ())::(readHeader ())
  1064               | "SECTION" => (pushToken t; [])
  1065               | _ => (readRestOfLine (); readHeader ())
  1066         else
  1067             (readRestOfLine (); readHeader ())
  1068         end
  1069 
  1070     fun skip_nls () =
  1071         let val x = readToken () in
  1072         if is_tt x TOKEN_NL then
  1073             skip_nls ()
  1074         else
  1075             (pushToken x; ())
  1076         end
  1077 
  1078     fun skip_paragraph () =
  1079         if is_tt (readToken ()) TOKEN_NL then
  1080         (if is_tt (readToken ()) TOKEN_NL then
  1081              skip_nls ()
  1082          else
  1083              skip_paragraph ())
  1084         else
  1085         skip_paragraph ()
  1086 
  1087     fun load_value () =
  1088         let
  1089         val t1 = readToken ()
  1090         val t1 = if is_tt t1 TOKEN_SYMBOL andalso snd (the t1) = "A" then readToken () else t1
  1091         in
  1092         if is_tt t1 TOKEN_NUM then
  1093             let
  1094             val name = readToken ()
  1095             val status = readToken ()
  1096             val value = readToken ()
  1097             in
  1098             if is_tt name TOKEN_SYMBOL andalso
  1099                is_tt status TOKEN_SYMBOL andalso
  1100                is_tt value TOKEN_NUM
  1101             then
  1102                 readRestOfLine (SOME (snd (the name), snd (the value)))
  1103             else
  1104                 raise (Load_cplexResult "column line expected")
  1105             end
  1106         else
  1107             (pushToken t1; NONE)
  1108         end
  1109 
  1110     fun load_values () =
  1111         let val v = load_value () in
  1112         if v = NONE then [] else (the v)::(load_values ())
  1113         end
  1114 
  1115     val header = readHeader ()
  1116 
  1117     val result =
  1118         case AList.lookup (op =) header "STATUS" of
  1119         SOME "INFEASIBLE" => Infeasible
  1120           | SOME "NONOPTIMAL" => Unbounded
  1121           | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) header "OBJECTIVE"),
  1122                        (skip_paragraph ();
  1123                         skip_paragraph ();
  1124                         skip_paragraph ();
  1125                         skip_paragraph ();
  1126                         skip_paragraph ();
  1127                         load_values ()))
  1128           | _ => Undefined
  1129 
  1130     val _ = TextIO.closeIn f
  1131     in
  1132     result
  1133     end
  1134     handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
  1135      | Option.Option => raise (Load_cplexResult "Option")
  1136 
  1137 exception Execute of string;
  1138 
  1139 fun tmp_file s = Path.implode (Path.expand (File.tmp_path (Path.basic s)));
  1140 fun wrap s = "\""^s^"\"";
  1141 
  1142 fun solve_glpk prog =
  1143     let
  1144     val name = string_of_int (Time.toMicroseconds (Time.now ()))
  1145     val lpname = tmp_file (name^".lp")
  1146     val resultname = tmp_file (name^".txt")
  1147     val _ = save_cplexFile lpname prog
  1148     val cplex_path = getenv "GLPK_PATH"
  1149     val cplex = if cplex_path = "" then "glpsol" else cplex_path
  1150     val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname)
  1151     val answer = #1 (Isabelle_System.bash_output command)
  1152     in
  1153     let
  1154         val result = load_glpkResult resultname
  1155         val _ = OS.FileSys.remove lpname
  1156         val _ = OS.FileSys.remove resultname
  1157     in
  1158         result
  1159     end
  1160     handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer))
  1161          | exn => if Exn.is_interrupt exn then reraise exn else raise (Execute answer)
  1162     end
  1163 
  1164 fun solve_cplex prog =
  1165     let
  1166     fun write_script s lp r =
  1167         let
  1168         val f = TextIO.openOut s
  1169         val _ = TextIO.output (f, "read\n"^lp^"\noptimize\nwrite\n"^r^"\nquit")
  1170         val _ = TextIO.closeOut f
  1171         in
  1172         ()
  1173         end
  1174 
  1175     val name = string_of_int (Time.toMicroseconds (Time.now ()))
  1176     val lpname = tmp_file (name^".lp")
  1177     val resultname = tmp_file (name^".txt")
  1178     val scriptname = tmp_file (name^".script")
  1179     val _ = save_cplexFile lpname prog
  1180     val _ = write_script scriptname lpname resultname
  1181     in
  1182     let
  1183         val result = load_cplexResult resultname
  1184         val _ = OS.FileSys.remove lpname
  1185         val _ = OS.FileSys.remove resultname
  1186         val _ = OS.FileSys.remove scriptname
  1187     in
  1188         result
  1189     end
  1190     end
  1191 
  1192 fun solve prog =
  1193     case get_solver () of
  1194       SOLVER_DEFAULT =>
  1195         (case getenv "LP_SOLVER" of
  1196        "CPLEX" => solve_cplex prog
  1197          | "GLPK" => solve_glpk prog
  1198          | _ => raise (Execute ("LP_SOLVER must be set to CPLEX or to GLPK")))
  1199     | SOLVER_CPLEX => solve_cplex prog
  1200     | SOLVER_GLPK => solve_glpk prog
  1201 
  1202 end;