# HG changeset patch # User obua # Date 1121241230 -7200 # Node ID 92ff7c90358574852075eb750e0e46871e5f0f86 # Parent 26fccaaf9cb458b344bac162b62dd6f6fa321b47 - added cplex package to HOL/Matrix diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Complex/ROOT.ML --- a/src/HOL/Complex/ROOT.ML Tue Jul 12 23:42:11 2005 +0200 +++ b/src/HOL/Complex/ROOT.ML Wed Jul 13 09:53:50 2005 +0200 @@ -5,6 +5,6 @@ The Complex Numbers *) -with_path "../Real" use_thy "Real"; +with_path "../Real" use_thy "Float"; with_path "../Hyperreal" use_thy "Hyperreal"; time_use_thy "Complex_Main"; diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/IsaMakefile --- a/src/HOL/IsaMakefile Tue Jul 12 23:42:11 2005 +0200 +++ b/src/HOL/IsaMakefile Wed Jul 13 09:53:50 2005 +0200 @@ -627,7 +627,10 @@ $(LOG)/HOL-Complex-Matrix.gz: $(OUT)/HOL-Complex \ Matrix/MatrixGeneral.thy Matrix/Matrix.thy Matrix/SparseMatrix.thy \ - Matrix/document/root.tex Matrix/ROOT.ML + Matrix/document/root.tex Matrix/ROOT.ML \ + Matrix/cplex/ROOT.ML Matrix/cplex/Cplex.thy Matrix/cplex/CplexMatrixConverter.ML \ + Matrix/cplex/Cplex_tools.ML Matrix/cplex/FloatSparseMatrix.thy \ + Matrix/cplex/FloatSparseMatrixBuilder.ML Matrix/cplex/fspmlp.ML @$(ISATOOL) usedir $(OUT)/HOL-Complex Matrix diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/ROOT.ML --- a/src/HOL/Matrix/ROOT.ML Tue Jul 12 23:42:11 2005 +0200 +++ b/src/HOL/Matrix/ROOT.ML Wed Jul 13 09:53:50 2005 +0200 @@ -3,4 +3,6 @@ *) use_thy "SparseMatrix"; +cd "cplex"; use_thy "Cplex"; cd ".."; + diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/Cplex.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/Cplex.thy Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,14 @@ +(* Title: HOL/Matrix/cplex/Cplex.thy + ID: $Id$ + Author: Steven Obua +*) + +theory Cplex +imports FloatSparseMatrix +files "Cplex_tools.ML" "CplexMatrixConverter.ML" "FloatSparseMatrixBuilder.ML" "fspmlp.ML" +begin + +end + + + diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/CplexMatrixConverter.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/CplexMatrixConverter.ML Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,137 @@ +(* Title: HOL/Matrix/cplex/CplexMatrixConverter.ML + ID: $Id$ + Author: Steven Obua +*) + +signature MATRIX_BUILDER = +sig + type vector + type matrix + + val empty_vector : vector + val empty_matrix : matrix + + exception Nat_expected of int + val set_elem : vector -> int -> string -> vector + val set_vector : matrix -> int -> vector -> matrix +end; + +signature CPLEX_MATRIX_CONVERTER = +sig + structure cplex : CPLEX + structure matrix_builder : MATRIX_BUILDER + type vector = matrix_builder.vector + type matrix = matrix_builder.matrix + type naming = int * (int -> string) * (string -> int) + + exception Converter of string + + (* program must fulfill is_normed_cplexProg and must be an element of the image of elim_nonfree_bounds *) + (* convert_prog maximize c A b naming *) + val convert_prog : cplex.cplexProg -> bool * vector * matrix * vector * naming + + (* results must be optimal, converts_results returns the optimal value as string and the solution as vector *) + (* convert_results results name2index *) + val convert_results : cplex.cplexResult -> (string -> int) -> string * vector +end; + +functor MAKE_CPLEX_MATRIX_CONVERTER (structure cplex: CPLEX and matrix_builder: MATRIX_BUILDER) : CPLEX_MATRIX_CONVERTER = +struct + +structure cplex = cplex +structure matrix_builder = matrix_builder +type matrix = matrix_builder.matrix +type vector = matrix_builder.vector +type naming = int * (int -> string) * (string -> int) + +open matrix_builder +open cplex + +exception Converter of string; + +structure Inttab = TableFun(type key = int val ord = int_ord); + +fun neg_term (cplexNeg t) = t + | neg_term (cplexSum ts) = cplexSum (map neg_term ts) + | neg_term t = cplexNeg t + +fun convert_prog (cplexProg (s, goal, constrs, bounds)) = + let + fun build_naming index i2s s2i [] = (index, i2s, s2i) + | build_naming index i2s s2i (cplexBounds (cplexNeg cplexInf, cplexLeq, cplexVar v, cplexLeq, cplexInf)::bounds) + = build_naming (index+1) (Inttab.update ((index, v), i2s)) (Symtab.update_new ((v, index), s2i)) bounds + | build_naming _ _ _ _ = raise (Converter "nonfree bound") + + val (varcount, i2s_tab, s2i_tab) = build_naming 0 Inttab.empty Symtab.empty bounds + + fun i2s i = case Inttab.lookup (i2s_tab, i) of NONE => raise (Converter "index not found") + | SOME n => n + fun s2i s = case Symtab.lookup (s2i_tab, s) of NONE => raise (Converter ("name not found: "^s)) + | SOME i => i + fun num2str positive (cplexNeg t) = num2str (not positive) t + | num2str positive (cplexNum num) = if positive then num else "-"^num + | num2str _ _ = raise (Converter "term is not a (possibly signed) number") + + fun setprod vec positive (cplexNeg t) = setprod vec (not positive) t + | setprod vec positive (cplexVar v) = set_elem vec (s2i v) (if positive then "1" else "-1") + | setprod vec positive (cplexProd (cplexNum num, cplexVar v)) = + set_elem vec (s2i v) (if positive then num else "-"^num) + | setprod _ _ _ = raise (Converter "term is not a normed product") + + fun sum2vec (cplexSum ts) = Library.foldl (fn (vec, t) => setprod vec true t) (empty_vector, ts) + | sum2vec t = setprod empty_vector true t + + fun constrs2Ab j A b [] = (A, b) + | constrs2Ab j A b ((_, cplexConstr (cplexLeq, (t1,t2)))::cs) = + constrs2Ab (j+1) (set_vector A j (sum2vec t1)) (set_elem b j (num2str true t2)) cs + | constrs2Ab j A b ((_, cplexConstr (cplexGeq, (t1,t2)))::cs) = + constrs2Ab (j+1) (set_vector A j (sum2vec (neg_term t1))) (set_elem b j (num2str true (neg_term t2))) cs + | constrs2Ab j A b ((_, cplexConstr (cplexEq, (t1,t2)))::cs) = + constrs2Ab j A b ((NONE, cplexConstr (cplexLeq, (t1,t2))):: + (NONE, cplexConstr (cplexGeq, (t1, t2)))::cs) + | constrs2Ab _ _ _ _ = raise (Converter "no strict constraints allowed") + + val (A, b) = constrs2Ab 0 empty_matrix empty_vector constrs + + val (goal_maximize, goal_term) = + case goal of + (cplexMaximize t) => (true, t) + | (cplexMinimize t) => (false, t) + in + (goal_maximize, sum2vec goal_term, A, b, (varcount, i2s, s2i)) + end + +fun convert_results (cplex.Optimal (opt, entries)) name2index = + let + fun setv (v, (name, value)) = (matrix_builder.set_elem v (name2index name) value) + in + (opt, Library.foldl setv (matrix_builder.empty_vector, entries)) + end + | convert_results _ _ = raise (Converter "No optimal result") + + +end; + +structure SimpleMatrixBuilder : MATRIX_BUILDER = +struct +type vector = (int * string) list +type matrix = (int * vector) list + +val empty_matrix = [] +val empty_vector = [] + +exception Nat_expected of int; + +fun set_elem v i s = v @ [(i, s)] + +fun set_vector m i v = m @ [(i, v)] + +end; + + + +structure SimpleCplexMatrixConverter = MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = SimpleMatrixBuilder); + + + + diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/Cplex_tools.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/Cplex_tools.ML Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,1220 @@ +(* Title: HOL/Matrix/cplex/Cplex_tools.ML + ID: $Id$ + Author: Steven Obua +*) + +signature CPLEX = +sig + + datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf + | cplexNeg of cplexTerm + | cplexProd of cplexTerm * cplexTerm + | cplexSum of (cplexTerm list) + + datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq + + datatype cplexGoal = cplexMinimize of cplexTerm + | cplexMaximize of cplexTerm + + datatype cplexConstr = cplexConstr of cplexComp * + (cplexTerm * cplexTerm) + + datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm + * cplexComp * cplexTerm + | cplexBound of cplexTerm * cplexComp * cplexTerm + + datatype cplexProg = cplexProg of string + * cplexGoal + * ((string option * cplexConstr) + list) + * cplexBounds list + + datatype cplexResult = Unbounded + | Infeasible + | Undefined + | Optimal of string * + (((* name *) string * + (* value *) string) list) + + exception Load_cplexFile of string + exception Load_cplexResult of string + exception Save_cplexFile of string + exception Execute of string + + val load_cplexFile : string -> cplexProg + + val save_cplexFile : string -> cplexProg -> unit + + val elim_nonfree_bounds : cplexProg -> cplexProg + + val relax_strict_ineqs : cplexProg -> cplexProg + + val is_normed_cplexProg : cplexProg -> bool + + val solve : cplexProg -> cplexResult +end; + +structure Cplex : CPLEX = +struct + +exception Load_cplexFile of string; +exception Load_cplexResult of string; +exception Save_cplexFile of string; + +datatype cplexTerm = cplexVar of string + | cplexNum of string + | cplexInf + | cplexNeg of cplexTerm + | cplexProd of cplexTerm * cplexTerm + | cplexSum of (cplexTerm list) +datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq +datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm +datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm) +datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm + * cplexComp * cplexTerm + | cplexBound of cplexTerm * cplexComp * cplexTerm +datatype cplexProg = cplexProg of string + * cplexGoal + * ((string option * cplexConstr) list) + * cplexBounds list + +fun rev_cmp cplexLe = cplexGe + | rev_cmp cplexLeq = cplexGeq + | rev_cmp cplexGe = cplexLe + | rev_cmp cplexGeq = cplexLeq + | rev_cmp cplexEq = cplexEq + +fun the NONE = raise (Load_cplexFile "SOME expected") + | the (SOME x) = x; + +fun modulo_signed is_something (cplexNeg u) = is_something u + | modulo_signed is_something u = is_something u + +fun is_Num (cplexNum _) = true + | is_Num _ = false + +fun is_Inf cplexInf = true + | is_Inf _ = false + +fun is_Var (cplexVar _) = true + | is_Var _ = false + +fun is_Neg (cplexNeg x ) = true + | is_Neg _ = false + +fun is_normed_Prod (cplexProd (t1, t2)) = + (is_Num t1) andalso (is_Var t2) + | is_normed_Prod x = is_Var x + +fun is_normed_Sum (cplexSum ts) = + (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts + | is_normed_Sum x = modulo_signed is_normed_Prod x + +fun is_normed_Constr (cplexConstr (c, (t1, t2))) = + (is_normed_Sum t1) andalso (modulo_signed is_Num t2) + +fun is_Num_or_Inf x = is_Inf x orelse is_Num x + +fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) = + (c1 = cplexLe orelse c1 = cplexLeq) andalso + (c2 = cplexLe orelse c2 = cplexLeq) andalso + is_Var t2 andalso + modulo_signed is_Num_or_Inf t1 andalso + modulo_signed is_Num_or_Inf t3 + | is_normed_Bounds (cplexBound (t1, c, t2)) = + (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2)) + orelse + (c <> cplexEq andalso + is_Var t2 andalso (modulo_signed is_Num_or_Inf t1)) + +fun term_of_goal (cplexMinimize x) = x + | term_of_goal (cplexMaximize x) = x + +fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) = + is_normed_Sum (term_of_goal goal) andalso + forall (fn (_,x) => is_normed_Constr x) constraints andalso + forall is_normed_Bounds bounds + +fun is_NL s = s = "\n" + +fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s) + +fun is_num a = + let + val b = String.explode a + fun num4 cs = forall Char.isDigit cs + fun num3 [] = true + | num3 (ds as (c::cs)) = + if c = #"+" orelse c = #"-" then + num4 cs + else + num4 ds + fun num2 [] = true + | num2 (c::cs) = + if c = #"e" orelse c = #"E" then num3 cs + else (Char.isDigit c) andalso num2 cs + fun num1 [] = true + | num1 (c::cs) = + if c = #"." then num2 cs + else if c = #"e" orelse c = #"E" then num3 cs + else (Char.isDigit c) andalso num1 cs + fun num [] = true + | num (c::cs) = + if c = #"." then num2 cs + else (Char.isDigit c) andalso num1 cs + in + num b + end + +fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":" + +fun is_cmp s = s = "<" orelse s = ">" orelse s = "<=" + orelse s = ">=" orelse s = "=" + +fun is_symbol a = + let + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~" + fun is_symbol_char c = Char.isAlphaNum c orelse + exists (fn d => d=c) symbol_char + fun is_symbol_start c = is_symbol_char c andalso + not (Char.isDigit c) andalso + not (c= #".") + val b = String.explode a + in + b <> [] andalso is_symbol_start (hd b) andalso + forall is_symbol_char b + end + +fun to_upper s = String.implode (map Char.toUpper (String.explode s)) + +fun keyword x = + let + val a = to_upper x + in + if a = "BOUNDS" orelse a = "BOUND" then + SOME "BOUNDS" + else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then + SOME "MINIMIZE" + else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then + SOME "MAXIMIZE" + else if a = "ST" orelse a = "S.T." orelse a = "ST." then + SOME "ST" + else if a = "FREE" orelse a = "END" then + SOME a + else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then + SOME "GENERAL" + else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then + SOME "INTEGER" + else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then + SOME "BINARY" + else if a = "INF" orelse a = "INFINITY" then + SOME "INF" + else + NONE + end + +val TOKEN_ERROR = ~1 +val TOKEN_BLANK = 0 +val TOKEN_NUM = 1 +val TOKEN_DELIMITER = 2 +val TOKEN_SYMBOL = 3 +val TOKEN_LABEL = 4 +val TOKEN_CMP = 5 +val TOKEN_KEYWORD = 6 +val TOKEN_NL = 7 + +(* tokenize takes a list of chars as argument and returns a list of + int * string pairs, each string representing a "cplex token", + and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP + or TOKEN_SYMBOL *) +fun tokenize s = + let + val flist = [(is_NL, TOKEN_NL), + (is_blank, TOKEN_BLANK), + (is_num, TOKEN_NUM), + (is_delimiter, TOKEN_DELIMITER), + (is_cmp, TOKEN_CMP), + (is_symbol, TOKEN_SYMBOL)] + fun match_helper [] s = (fn x => false, TOKEN_ERROR) + | match_helper (f::fs) s = + if ((fst f) s) then f else match_helper fs s + fun match s = match_helper flist s + fun tok s = + if s = "" then [] else + let + val h = String.substring (s,0,1) + val (f, j) = match h + fun len i = + if size s = i then i + else if f (String.substring (s,0,i+1)) then + len (i+1) + else i + in + if j < 0 then + (if h = "\\" then [] + else raise (Load_cplexFile ("token expected, found: " + ^s))) + else + let + val l = len 1 + val u = String.substring (s,0,l) + val v = String.extract (s,l,NONE) + in + if j = 0 then tok v else (j, u) :: tok v + end + end + in + tok s + end + +exception Tokenize of string; + +fun tokenize_general flist s = + let + fun match_helper [] s = raise (Tokenize s) + | match_helper (f::fs) s = + if ((fst f) s) then f else match_helper fs s + fun match s = match_helper flist s + fun tok s = + if s = "" then [] else + let + val h = String.substring (s,0,1) + val (f, j) = match h + fun len i = + if size s = i then i + else if f (String.substring (s,0,i+1)) then + len (i+1) + else i + val l = len 1 + in + (j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE)) + end + in + tok s + end + +fun load_cplexFile name = + let + val f = TextIO.openIn name + val ignore_NL = ref true + val rest = ref [] + + fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + SOME u + ) + end + else + let val s = TextIO.inputLine f in + if s = "" then NONE else + let val t = tokenize s in + if (length t >= 2 andalso + snd(hd (tl t)) = ":") + then + rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t)) + else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t)) + andalso is_symbol "TO" (hd (tl t)) + then + rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t)) + else + rest := t; + readToken_helper () + end + end + + fun readToken_helper2 () = + let val c = readToken_helper () in + if c = NONE then NONE + else if !ignore_NL andalso fst (the c) = TOKEN_NL then + readToken_helper2 () + else if fst (the c) = TOKEN_SYMBOL + andalso keyword (snd (the c)) <> NONE + then SOME (TOKEN_KEYWORD, the (keyword (snd (the c)))) + else c + end + + fun readToken () = readToken_helper2 () + + fun pushToken a = rest := (a::(!rest)) + + fun is_value token = + fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD + andalso snd token = "INF") + + fun get_value token = + if fst token = TOKEN_NUM then + cplexNum (snd token) + else if fst token = TOKEN_KEYWORD andalso snd token = "INF" + then + cplexInf + else + raise (Load_cplexFile "num expected") + + fun readTerm_Product only_num = + let val c = readToken () in + if c = NONE then NONE + else if fst (the c) = TOKEN_SYMBOL + then ( + if only_num then (pushToken (the c); NONE) + else SOME (cplexVar (snd (the c))) + ) + else if only_num andalso is_value (the c) then + SOME (get_value (the c)) + else if is_value (the c) then + let val t1 = get_value (the c) + val d = readToken () + in + if d = NONE then SOME t1 + else if fst (the d) = TOKEN_SYMBOL then + SOME (cplexProd (t1, cplexVar (snd (the d)))) + else + (pushToken (the d); SOME t1) + end + else (pushToken (the c); NONE) + end + + fun readTerm_Signed only_signed only_num = + let + val c = readToken () + in + if c = NONE then NONE + else + let val d = the c in + if d = (TOKEN_DELIMITER, "+") then + readTerm_Product only_num + else if d = (TOKEN_DELIMITER, "-") then + SOME (cplexNeg (the (readTerm_Product + only_num))) + else (pushToken d; + if only_signed then NONE + else readTerm_Product only_num) + end + end + + fun readTerm_Sum first_signed = + let val c = readTerm_Signed first_signed false in + if c = NONE then [] else (the c)::(readTerm_Sum true) + end + + fun readTerm () = + let val c = readTerm_Sum false in + if c = [] then NONE + else if tl c = [] then SOME (hd c) + else SOME (cplexSum c) + end + + fun readLabeledTerm () = + let val c = readToken () in + if c = NONE then (NONE, NONE) + else if fst (the c) = TOKEN_LABEL then + let val t = readTerm () in + if t = NONE then + raise (Load_cplexFile ("term after label "^ + (snd (the c))^ + " expected")) + else (SOME (snd (the c)), t) + end + else (pushToken (the c); (NONE, readTerm ())) + end + + fun readGoal () = + let + val g = readToken () + in + if g = SOME (TOKEN_KEYWORD, "MAXIMIZE") then + cplexMaximize (the (snd (readLabeledTerm ()))) + else if g = SOME (TOKEN_KEYWORD, "MINIMIZE") then + cplexMinimize (the (snd (readLabeledTerm ()))) + else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected") + end + + fun str2cmp b = + (case b of + "<" => cplexLe + | "<=" => cplexLeq + | ">" => cplexGe + | ">=" => cplexGeq + | "=" => cplexEq + | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP"))) + + fun readConstraint () = + let + val t = readLabeledTerm () + fun make_constraint b t1 t2 = + cplexConstr + (str2cmp b, + (t1, t2)) + in + if snd t = NONE then NONE + else + let val c = readToken () in + if c = NONE orelse fst (the c) <> TOKEN_CMP + then raise (Load_cplexFile "TOKEN_CMP expected") + else + let val n = readTerm_Signed false true in + if n = NONE then + raise (Load_cplexFile "num expected") + else + SOME (fst t, + make_constraint (snd (the c)) + (the (snd t)) + (the n)) + end + end + end + + fun readST () = + let + fun readbody () = + let val t = readConstraint () in + if t = NONE then [] + else if (is_normed_Constr (snd (the t))) then + (the t)::(readbody ()) + else if (fst (the t) <> NONE) then + raise (Load_cplexFile + ("constraint '"^(the (fst (the t)))^ + "'is not normed")) + else + raise (Load_cplexFile + "constraint is not normed") + end + in + if readToken () = SOME (TOKEN_KEYWORD, "ST") + then + readbody () + else + raise (Load_cplexFile "ST expected") + end + + fun readCmp () = + let val c = readToken () in + if c = NONE then NONE + else if fst (the c) = TOKEN_CMP then + SOME (str2cmp (snd (the c))) + else (pushToken (the c); NONE) + end + + fun skip_NL () = + let val c = readToken () in + if c <> NONE andalso fst (the c) = TOKEN_NL then + skip_NL () + else + (pushToken (the c); ()) + end + + fun is_var (cplexVar _) = true + | is_var _ = false + + fun make_bounds c t1 t2 = + cplexBound (t1, c, t2) + + fun readBound () = + let + val _ = skip_NL () + val t1 = readTerm () + in + if t1 = NONE then NONE + else + let + val c1 = readCmp () + in + if c1 = NONE then + let + val c = readToken () + in + if c = SOME (TOKEN_KEYWORD, "FREE") then + SOME ( + cplexBounds (cplexNeg cplexInf, + cplexLeq, + the t1, + cplexLeq, + cplexInf)) + else + raise (Load_cplexFile "FREE expected") + end + else + let + val t2 = readTerm () + in + if t2 = NONE then + raise (Load_cplexFile "term expected") + else + let val c2 = readCmp () in + if c2 = NONE then + SOME (make_bounds (the c1) + (the t1) + (the t2)) + else + SOME ( + cplexBounds (the t1, + the c1, + the t2, + the c2, + the (readTerm()))) + end + end + end + end + + fun readBounds () = + let + fun makestring b = "?" + fun readbody () = + let + val b = readBound () + in + if b = NONE then [] + else if (is_normed_Bounds (the b)) then + (the b)::(readbody()) + else ( + raise (Load_cplexFile + ("bounds are not normed in: "^ + (makestring (the b))))) + end + in + if readToken () = SOME (TOKEN_KEYWORD, "BOUNDS") then + readbody () + else raise (Load_cplexFile "BOUNDS expected") + end + + fun readEnd () = + if readToken () = SOME (TOKEN_KEYWORD, "END") then () + else raise (Load_cplexFile "END expected") + + val result_Goal = readGoal () + val result_ST = readST () + val _ = ignore_NL := false + val result_Bounds = readBounds () + val _ = ignore_NL := true + val _ = readEnd () + val _ = TextIO.closeIn f + in + cplexProg (name, result_Goal, result_ST, result_Bounds) + end + +fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) = + let + val f = TextIO.openOut filename + + fun basic_write s = TextIO.output(f, s) + + val linebuf = ref "" + fun buf_flushline s = + (basic_write (!linebuf); + basic_write "\n"; + linebuf := s) + fun buf_add s = linebuf := (!linebuf) ^ s + + fun write s = + if (String.size s) + (String.size (!linebuf)) >= 250 then + buf_flushline (" "^s) + else + buf_add s + + fun writeln s = (buf_add s; buf_flushline "") + + fun write_term (cplexVar x) = write x + | write_term (cplexNum x) = write x + | write_term cplexInf = write "inf" + | write_term (cplexProd (cplexNum "1", b)) = write_term b + | write_term (cplexProd (a, b)) = + (write_term a; write " "; write_term b) + | write_term (cplexNeg x) = (write " - "; write_term x) + | write_term (cplexSum ts) = write_terms ts + and write_terms [] = () + | write_terms (t::ts) = + (if (not (is_Neg t)) then write " + " else (); + write_term t; write_terms ts) + + fun write_goal (cplexMaximize term) = + (writeln "MAXIMIZE"; write_term term; writeln "") + | write_goal (cplexMinimize term) = + (writeln "MINIMIZE"; write_term term; writeln "") + + fun write_cmp cplexLe = write "<" + | write_cmp cplexLeq = write "<=" + | write_cmp cplexEq = write "=" + | write_cmp cplexGe = write ">" + | write_cmp cplexGeq = write ">=" + + fun write_constr (cplexConstr (cmp, (a,b))) = + (write_term a; + write " "; + write_cmp cmp; + write " "; + write_term b) + + fun write_constraints [] = () + | write_constraints (c::cs) = + (if (fst c <> NONE) + then + (write (the (fst c)); write ": ") + else + (); + write_constr (snd c); + writeln ""; + write_constraints cs) + + fun write_bounds [] = () + | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) = + ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf + andalso (c1 = cplexLeq orelse c1 = cplexLe) + andalso (c2 = cplexLeq orelse c2 = cplexLe) + then + (write_term t2; write " free") + else + (write_term t1; write " "; write_cmp c1; write " "; + write_term t2; write " "; write_cmp c2; write " "; + write_term t3) + ); writeln ""; write_bounds bs) + | write_bounds ((cplexBound (t1, c, t2)) :: bs) = + (write_term t1; write " "; + write_cmp c; write " "; + write_term t2; writeln ""; write_bounds bs) + + val _ = write_goal goal + val _ = (writeln ""; writeln "ST") + val _ = write_constraints constraints + val _ = (writeln ""; writeln "BOUNDS") + val _ = write_bounds bounds + val _ = (writeln ""; writeln "END") + val _ = TextIO.closeOut f + in + () + end + +fun norm_Constr (constr as cplexConstr (c, (t1, t2))) = + if not (modulo_signed is_Num t2) andalso + modulo_signed is_Num t1 + then + [cplexConstr (rev_cmp c, (t2, t1))] + else if (c = cplexLe orelse c = cplexLeq) andalso + (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf) + then + [] + else if (c = cplexGe orelse c = cplexGeq) andalso + (t1 = cplexInf orelse t2 = cplexNeg cplexInf) + then + [] + else + [constr] + +fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) = + (norm_Constr(cplexConstr (c1, (t1, t2)))) + @ (norm_Constr(cplexConstr (c2, (t2, t3)))) + | bound2constr (cplexBound (t1, cplexEq, t2)) = + (norm_Constr(cplexConstr (cplexLeq, (t1, t2)))) + @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1)))) + | bound2constr (cplexBound (t1, c1, t2)) = + norm_Constr(cplexConstr (c1, (t1,t2))) + +val emptyset = Symtab.empty + +fun singleton v = Symtab.update ((v, ()), emptyset) + +fun merge a b = Symtab.merge (op =) (a, b) + +fun mergemap f ts = Library.foldl + (fn (table, x) => merge table (f x)) + (Symtab.empty, ts) + +fun diff a b = Symtab.foldl (fn (a, (k,v)) => + (Symtab.delete k a) handle UNDEF => a) + (a, b) + +fun collect_vars (cplexVar v) = singleton v + | collect_vars (cplexNeg t) = collect_vars t + | collect_vars (cplexProd (t1, t2)) = + merge (collect_vars t1) (collect_vars t2) + | collect_vars (cplexSum ts) = mergemap collect_vars ts + | collect_vars _ = emptyset + +(* Eliminates all nonfree bounds from the linear program and produces an + equivalent program with only free bounds + IF for the input program P holds: is_normed_cplexProg P *) +fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) = + let + fun collect_constr_vars (_, cplexConstr (c, (t1,_))) = + (collect_vars t1) + + val cvars = merge (collect_vars (term_of_goal goal)) + (mergemap collect_constr_vars constraints) + + fun collect_lower_bounded_vars + (cplexBounds (t1, c1, cplexVar v, c2, t3)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (_, cplexLe, cplexVar v)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (_, cplexLeq, cplexVar v)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexGe,_)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexGeq, _)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexEq, _)) = + singleton v + | collect_lower_bounded_vars _ = emptyset + + val lvars = mergemap collect_lower_bounded_vars bounds + val positive_vars = diff cvars lvars + val zero = cplexNum "0" + + fun make_pos_constr v = + (NONE, cplexConstr (cplexGeq, ((cplexVar v), zero))) + + fun make_free_bound v = + cplexBounds (cplexNeg cplexInf, cplexLeq, + cplexVar v, cplexLeq, + cplexInf) + + val pos_constrs = rev(Symtab.foldl + (fn (l, (k,v)) => (make_pos_constr k) :: l) + ([], positive_vars)) + val bound_constrs = map (fn x => (NONE, x)) + (Library.flat (map bound2constr bounds)) + val constraints' = constraints @ pos_constrs @ bound_constrs + val bounds' = rev(Symtab.foldl (fn (l, (v,_)) => + (make_free_bound v)::l) + ([], cvars)) + in + cplexProg (name, goal, constraints', bounds') + end + +fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) = + let + fun relax cplexLe = cplexLeq + | relax cplexGe = cplexGeq + | relax x = x + + fun relax_constr (n, cplexConstr(c, (t1, t2))) = + (n, cplexConstr(relax c, (t1, t2))) + + fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) = + cplexBounds (t1, relax c1, t2, relax c2, t3) + | relax_bounds (cplexBound (t1, c, t2)) = + cplexBound (t1, relax c, t2) + in + cplexProg (name, + goals, + map relax_constr constrs, + map relax_bounds bounds) + end + +datatype cplexResult = Unbounded + | Infeasible + | Undefined + | Optimal of string * ((string * string) list) + +fun is_separator x = forall (fn c => c = #"-") (String.explode x) + +fun is_sign x = (x = "+" orelse x = "-") + +fun is_colon x = (x = ":") + +fun is_resultsymbol a = + let + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-" + fun is_symbol_char c = Char.isAlphaNum c orelse + exists (fn d => d=c) symbol_char + fun is_symbol_start c = is_symbol_char c andalso + not (Char.isDigit c) andalso + not (c= #".") andalso + not (c= #"-") + val b = String.explode a + in + b <> [] andalso is_symbol_start (hd b) andalso + forall is_symbol_char b + end + +val TOKEN_SIGN = 100 +val TOKEN_COLON = 101 +val TOKEN_SEPARATOR = 102 + +fun load_glpkResult name = + let + val flist = [(is_NL, TOKEN_NL), + (is_blank, TOKEN_BLANK), + (is_num, TOKEN_NUM), + (is_sign, TOKEN_SIGN), + (is_colon, TOKEN_COLON), + (is_cmp, TOKEN_CMP), + (is_resultsymbol, TOKEN_SYMBOL), + (is_separator, TOKEN_SEPARATOR)] + + val tokenize = tokenize_general flist + + val f = TextIO.openIn name + + val rest = ref [] + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + SOME u + ) + end + else + let val s = TextIO.inputLine f in + if s = "" then NONE else (rest := tokenize s; readToken_helper()) + end + + fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty) + + fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest))) + + fun readToken () = + let val t = readToken_helper () in + if is_tt t TOKEN_BLANK then + readToken () + else if is_tt t TOKEN_NL then + let val t2 = readToken_helper () in + if is_tt t2 TOKEN_SIGN then + (pushToken (SOME (TOKEN_SEPARATOR, "-")); t) + else + (pushToken t2; t) + end + else if is_tt t TOKEN_SIGN then + let val t2 = readToken_helper () in + if is_tt t2 TOKEN_NUM then + (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2)))) + else + (pushToken t2; t) + end + else + t + end + + fun readRestOfLine P = + let + val t = readToken () + in + if is_tt t TOKEN_NL orelse t = NONE + then P + else readRestOfLine P + end + + fun readHeader () = + let + fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ()))) + fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ()))) + val t1 = readToken () + val t2 = readToken () + in + if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON + then + case to_upper (snd (the t1)) of + "STATUS" => (readStatus ())::(readHeader ()) + | "OBJECTIVE" => (readObjective())::(readHeader ()) + | _ => (readRestOfLine (); readHeader ()) + else + (pushToken t2; pushToken t1; []) + end + + fun skip_until_sep () = + let val x = readToken () in + if is_tt x TOKEN_SEPARATOR then + readRestOfLine () + else + skip_until_sep () + end + + fun load_value () = + let + val t1 = readToken () + val t2 = readToken () + in + if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then + let + val t = readToken () + val state = if is_tt t TOKEN_NL then readToken () else t + val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected") + val k = readToken () + in + if is_tt k TOKEN_NUM then + readRestOfLine (SOME (snd (the t2), snd (the k))) + else + raise (Load_cplexResult "number expected") + end + else + (pushToken t2; pushToken t1; NONE) + end + + fun load_values () = + let val v = load_value () in + if v = NONE then [] else (the v)::(load_values ()) + end + + val header = readHeader () + + val result = + case assoc (header, "STATUS") of + SOME "INFEASIBLE" => Infeasible + | SOME "UNBOUNDED" => Unbounded + | SOME "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")), + (skip_until_sep (); + skip_until_sep (); + load_values ())) + | _ => Undefined + + val _ = TextIO.closeIn f + in + result + end + handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s)) + | Option => raise (Load_cplexResult "Option") + | x => raise x + +fun load_cplexResult name = + let + val flist = [(is_NL, TOKEN_NL), + (is_blank, TOKEN_BLANK), + (is_num, TOKEN_NUM), + (is_sign, TOKEN_SIGN), + (is_colon, TOKEN_COLON), + (is_cmp, TOKEN_CMP), + (is_resultsymbol, TOKEN_SYMBOL)] + + val tokenize = tokenize_general flist + + val f = TextIO.openIn name + + val rest = ref [] + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + SOME u + ) + end + else + let + val s = TextIO.inputLine f + in + if s = "" then NONE else (rest := tokenize s; readToken_helper()) + end + + fun is_tt tok ty = (tok <> NONE andalso (fst (the tok)) = ty) + + fun pushToken a = if a = NONE then () else (rest := ((the a)::(!rest))) + + fun readToken () = + let val t = readToken_helper () in + if is_tt t TOKEN_BLANK then + readToken () + else if is_tt t TOKEN_SIGN then + let val t2 = readToken_helper () in + if is_tt t2 TOKEN_NUM then + (SOME (TOKEN_NUM, (snd (the t))^(snd (the t2)))) + else + (pushToken t2; t) + end + else + t + end + + fun readRestOfLine P = + let + val t = readToken () + in + if is_tt t TOKEN_NL orelse t = NONE + then P + else readRestOfLine P + end + + fun readHeader () = + let + fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ()))) + fun readObjective () = + let + val t = readToken () + in + if is_tt t TOKEN_SYMBOL andalso to_upper (snd (the t)) = "VALUE" then + readRestOfLine ("OBJECTIVE", snd (the (readToken()))) + else + readRestOfLine ("OBJECTIVE_NAME", snd (the t)) + end + + val t = readToken () + in + if is_tt t TOKEN_SYMBOL then + case to_upper (snd (the t)) of + "STATUS" => (readStatus ())::(readHeader ()) + | "OBJECTIVE" => (readObjective ())::(readHeader ()) + | "SECTION" => (pushToken t; []) + | _ => (readRestOfLine (); readHeader ()) + else + (readRestOfLine (); readHeader ()) + end + + fun skip_nls () = + let val x = readToken () in + if is_tt x TOKEN_NL then + skip_nls () + else + (pushToken x; ()) + end + + fun skip_paragraph () = + if is_tt (readToken ()) TOKEN_NL then + (if is_tt (readToken ()) TOKEN_NL then + skip_nls () + else + skip_paragraph ()) + else + skip_paragraph () + + fun load_value () = + let + val t1 = readToken () + val t1 = if is_tt t1 TOKEN_SYMBOL andalso snd (the t1) = "A" then readToken () else t1 + in + if is_tt t1 TOKEN_NUM then + let + val name = readToken () + val status = readToken () + val value = readToken () + in + if is_tt name TOKEN_SYMBOL andalso + is_tt status TOKEN_SYMBOL andalso + is_tt value TOKEN_NUM + then + readRestOfLine (SOME (snd (the name), snd (the value))) + else + raise (Load_cplexResult "column line expected") + end + else + (pushToken t1; NONE) + end + + fun load_values () = + let val v = load_value () in + if v = NONE then [] else (the v)::(load_values ()) + end + + val header = readHeader () + + val result = + case assoc (header, "STATUS") of + SOME "INFEASIBLE" => Infeasible + | SOME "NONOPTIMAL" => Unbounded + | SOME "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")), + (skip_paragraph (); + skip_paragraph (); + skip_paragraph (); + skip_paragraph (); + skip_paragraph (); + load_values ())) + | _ => Undefined + + val _ = TextIO.closeIn f + in + result + end + handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s)) + | Option => raise (Load_cplexResult "Option") + | x => raise x + +exception Execute of string; + +fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s]))); +fun wrap s = "\""^s^"\""; + +fun solve_glpk prog = + let + val name = LargeInt.toString (Time.toMicroseconds (Time.now ())) + val lpname = tmp_file (name^".lp") + val resultname = tmp_file (name^".txt") + val _ = save_cplexFile lpname prog + val cplex_path = getenv "LP_SOLVER_PATH" + val cplex = if cplex_path = "" then "glpsol" else cplex_path + val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname) + val answer = execute command + in + let + val result = load_glpkResult resultname + val _ = OS.FileSys.remove lpname + val _ = OS.FileSys.remove resultname + in + result + end + handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer)) + | _ => raise (Execute answer) + end + +fun solve_cplex prog = + let + fun write_script s lp r = + let + val f = TextIO.openOut s + val _ = TextIO.output (f, "read\n"^lp^"\noptimize\nwrite\n"^r^"\nquit") + val _ = TextIO.closeOut f + in + () + end + + val name = LargeInt.toString (Time.toMicroseconds (Time.now ())) + val lpname = tmp_file (name^".lp") + val resultname = tmp_file (name^".txt") + val scriptname = tmp_file (name^".script") + val _ = save_cplexFile lpname prog + val cplex_path = getenv "LP_SOLVER_PATH" + val cplex = if cplex_path = "" then "cplex" else cplex_path + val _ = write_script scriptname lpname resultname + val command = (wrap cplex)^" < "^(wrap scriptname)^" > /dev/null" + val answer = "return code "^(Int.toString (system command)) + in + let + val result = load_cplexResult resultname + val _ = OS.FileSys.remove lpname + val _ = OS.FileSys.remove resultname + val _ = OS.FileSys.remove scriptname + in + result + end + end + +fun solve prog = + case getenv "LP_SOLVER_NAME" of + "CPLEX" => solve_cplex prog + | "GLPK" => solve_glpk prog + | _ => raise (Execute ("LP_SOLVER.NAME must be set to CPLEX or to GLPK")); + +end; + +(* +val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45" +val demoout = "/home/obua/flyspeck/kepler/LP/test.out" +val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol" + +fun loadcplex () = Cplex.relax_strict_ineqs + (Cplex.load_cplexFile demofile) + +fun writecplex lp = Cplex.save_cplexFile demoout lp + +fun test () = + let + val lp = loadcplex () + val lp2 = Cplex.elim_nonfree_bounds lp + in + writecplex lp2 + end + +fun loadresult () = Cplex.load_cplexResult demoresult; +*) + +(*val prog = Cplex.load_cplexFile "/home/obua/tmp/pent/graph_0.lpt"; +val _ = Cplex.solve prog;*) diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/FloatSparseMatrix.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/FloatSparseMatrix.thy Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,8 @@ +(* Title: HOL/Matrix/cplex/FloatSparseMatrix.thy + ID: $Id$ + Author: Steven Obua +*) + +theory FloatSparseMatrix = Float + SparseMatrix: + +end diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,393 @@ +(* Title: HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML + ID: $Id$ + Author: Steven Obua +*) + +structure FloatSparseMatrixBuilder : +sig + include MATRIX_BUILDER + + structure cplex : CPLEX + + type float = IntInf.int*IntInf.int + type floatfunc = float -> float + + + val float2cterm : IntInf.int * IntInf.int -> cterm + + val approx_value : int -> floatfunc -> string -> cterm * cterm + val approx_vector : int -> floatfunc -> vector -> cterm * cterm + val approx_matrix : int -> floatfunc -> matrix -> cterm * cterm + + val mk_spvec_entry : int -> float -> term + val empty_spvec : term + val cons_spvec : term -> term -> term + val empty_spmat : term + val mk_spmat_entry : int -> term -> term + val cons_spmat : term -> term -> term + val sign_term : term -> cterm + + val v_elem_at : vector -> int -> string option + val m_elem_at : matrix -> int -> vector option + val v_only_elem : vector -> int option + val v_fold : ('a * (int * string) -> 'a) -> 'a -> vector -> 'a + val m_fold : ('a * (int * vector) -> 'a) -> 'a -> matrix -> 'a + + val transpose_matrix : matrix -> matrix + + val cut_vector : int -> vector -> vector + val cut_matrix : vector -> (int option) -> matrix -> matrix + + (* cplexProg c A b *) + val cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int)) + (* dual_cplexProg c A b *) + val dual_cplexProg : vector -> matrix -> vector -> (cplex.cplexProg * (string -> int)) + + val real_spmatT : typ + val real_spvecT : typ +end += +struct + + +structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord)); + +type vector = string Inttab.table +type matrix = vector Inttab.table +type float = IntInf.int*IntInf.int +type floatfunc = float -> float + +val th = theory "FloatSparseMatrix" +val sg = sign_of th + +fun readtype s = Sign.intern_type sg s +fun readterm s = Sign.intern_const sg s + +val ty_list = readtype "list" +val term_Nil = readterm "Nil" +val term_Cons = readterm "Cons" + +val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT) +val spvecT = Type (ty_list, [spvec_elemT]) +val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT) +val spmatT = Type (ty_list, [spmat_elemT]) + +val real_spmatT = spmatT +val real_spvecT = spvecT + +val empty_matrix_const = Const (term_Nil, spmatT) +val empty_vector_const = Const (term_Nil, spvecT) + +val Cons_spvec_const = Const (term_Cons, spvec_elemT --> spvecT --> spvecT) +val Cons_spmat_const = Const (term_Cons, spmat_elemT --> spmatT --> spmatT) + +val float_const = Const (readterm "float", HOLogic.mk_prodT (HOLogic.intT, HOLogic.intT) --> HOLogic.realT) + +val zero = IntInf.fromInt 0 +val minus_one = IntInf.fromInt ~1 +val two = IntInf.fromInt 2 + +fun mk_intinf ty n = + let + fun mk_bit n = if n = zero then HOLogic.false_const else HOLogic.true_const + + fun bin_of n = + if n = zero then HOLogic.pls_const + else if n = minus_one then HOLogic.min_const + else + let + (*val (q,r) = IntInf.divMod (n, two): doesn't work in SML 10.0.7, but in newer versions!!!*) + val q = IntInf.div (n, two) + val r = IntInf.mod (n, two) + in + HOLogic.bit_const $ bin_of q $ mk_bit r + end + in + HOLogic.number_of_const ty $ (bin_of n) + end + +fun mk_float (a,b) = + float_const $ (HOLogic.mk_prod ((mk_intinf HOLogic.intT a), (mk_intinf HOLogic.intT b))) + +fun float2cterm (a,b) = cterm_of sg (mk_float (a,b)) + +fun approx_value_term prec f value = + let + val (flower, fupper) = ExactFloatingPoint.approx_decstr_by_bin prec value + val (flower, fupper) = (f flower, f fupper) + in + (mk_float flower, mk_float fupper) + end + +fun approx_value prec pprt value = + let + val (flower, fupper) = approx_value_term prec pprt value + in + (cterm_of sg flower, cterm_of sg fupper) + end + +fun sign_term t = cterm_of sg t + +val empty_spvec = empty_vector_const + +val empty_spmat = empty_matrix_const + +fun mk_spvec_entry i f = + let + val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i) + val term_f = mk_float f + in + HOLogic.mk_prod (term_i, term_f) + end + +fun mk_spmat_entry i e = + let + val term_i = mk_intinf HOLogic.natT (IntInf.fromInt i) + in + HOLogic.mk_prod (term_i, e) + end + +fun cons_spvec h t = Cons_spvec_const $ h $ t + +fun cons_spmat h t = Cons_spmat_const $ h $ t + +fun approx_vector_term prec pprt vector = + let + fun app ((vlower, vupper), (index, s)) = + let + val (flower, fupper) = approx_value_term prec pprt s + val index = mk_intinf HOLogic.natT (IntInf.fromInt index) + val elower = HOLogic.mk_prod (index, flower) + val eupper = HOLogic.mk_prod (index, fupper) + in + (Cons_spvec_const $ elower $ vlower, + Cons_spvec_const $ eupper $ vupper) + end + in + Inttab.foldl app ((empty_vector_const, empty_vector_const), vector) + end + +fun approx_matrix_term prec pprt matrix = + let + fun app ((mlower, mupper), (index, vector)) = + let + val (vlower, vupper) = approx_vector_term prec pprt vector + val index = mk_intinf HOLogic.natT (IntInf.fromInt index) + val elower = HOLogic.mk_prod (index, vlower) + val eupper = HOLogic.mk_prod (index, vupper) + in + (Cons_spmat_const $ elower $ mlower, + Cons_spmat_const $ eupper $ mupper) + end + + val (mlower, mupper) = Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix) + in + Inttab.foldl app ((empty_matrix_const, empty_matrix_const), matrix) + end + +fun approx_vector prec pprt vector = + let + val (l, u) = approx_vector_term prec pprt vector + in + (cterm_of sg l, cterm_of sg u) + end + +fun approx_matrix prec pprt matrix = + let + val (l, u) = approx_matrix_term prec pprt matrix + in + (cterm_of sg l, cterm_of sg u) + end + + +exception Nat_expected of int; + +val zero_interval = approx_value_term 1 I "0" + +fun set_elem vector index str = + if index < 0 then + raise (Nat_expected index) + else if (approx_value_term 1 I str) = zero_interval then + vector + else + Inttab.update ((index, str), vector) + +fun set_vector matrix index vector = + if index < 0 then + raise (Nat_expected index) + else if Inttab.is_empty vector then + matrix + else + Inttab.update ((index, vector), matrix) + +val empty_matrix = Inttab.empty +val empty_vector = Inttab.empty + +(* dual stuff *) + +structure cplex = Cplex + +fun transpose_matrix matrix = + let + fun upd m j i x = + case Inttab.lookup (m, j) of + SOME v => Inttab.update ((j, Inttab.update ((i, x), v)), m) + | NONE => Inttab.update ((j, Inttab.update ((i, x), Inttab.empty)), m) + + fun updv j (m, (i, s)) = upd m i j s + + fun updm (m, (j, v)) = Inttab.foldl (updv j) (m, v) + in + Inttab.foldl updm (empty_matrix, matrix) + end + +exception No_name of string; + +exception Superfluous_constr_right_hand_sides + +fun cplexProg c A b = + let + val ytable = ref Inttab.empty + fun indexof s = + if String.size s = 0 then raise (No_name s) + else case Int.fromString (String.extract(s, 1, NONE)) of + SOME i => i | NONE => raise (No_name s) + + fun nameof i = + let + val s = "x"^(Int.toString i) + val _ = ytable := (Inttab.update ((i, s), !ytable)) + in + s + end + + fun split_numstr s = + if String.isPrefix "-" s then (false,String.extract(s, 1, NONE)) + else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) + else (true, s) + + fun mk_term index s = + let + val (p, s) = split_numstr s + val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) + in + if p then prod else cplex.cplexNeg prod + end + + fun vec2sum vector = + cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector)) + + fun mk_constr index vector c = + let + val s = case Inttab.lookup (c, index) of SOME s => s | NONE => "0" + val (p, s) = split_numstr s + val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) + in + (NONE, cplex.cplexConstr (cplex.cplexLeq, (vec2sum vector, num))) + end + + fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c + + val (list, b) = Inttab.foldl + (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c)) + (([], b), A) + val _ = if Inttab.is_empty b then () else raise Superfluous_constr_right_hand_sides + + fun mk_free y = cplex.cplexBounds (cplex.cplexNeg cplex.cplexInf, cplex.cplexLeq, + cplex.cplexVar y, cplex.cplexLeq, + cplex.cplexInf) + + val yvars = Inttab.foldl (fn (l, (i, y)) => (mk_free y)::l) ([], !ytable) + + val prog = cplex.cplexProg ("original", cplex.cplexMaximize (vec2sum c), list, yvars) + in + (prog, indexof) + end + + +fun dual_cplexProg c A b = + let + fun indexof s = + if String.size s = 0 then raise (No_name s) + else case Int.fromString (String.extract(s, 1, NONE)) of + SOME i => i | NONE => raise (No_name s) + + fun nameof i = "y"^(Int.toString i) + + fun split_numstr s = + if String.isPrefix "-" s then (false,String.extract(s, 1, NONE)) + else if String.isPrefix "+" s then (true, String.extract(s, 1, NONE)) + else (true, s) + + fun mk_term index s = + let + val (p, s) = split_numstr s + val prod = cplex.cplexProd (cplex.cplexNum s, cplex.cplexVar (nameof index)) + in + if p then prod else cplex.cplexNeg prod + end + + fun vec2sum vector = + cplex.cplexSum (Inttab.foldl (fn (list, (index, s)) => (mk_term index s)::list) ([], vector)) + + fun mk_constr index vector c = + let + val s = case Inttab.lookup (c, index) of SOME s => s | NONE => "0" + val (p, s) = split_numstr s + val num = if p then cplex.cplexNum s else cplex.cplexNeg (cplex.cplexNum s) + in + (NONE, cplex.cplexConstr (cplex.cplexEq, (vec2sum vector, num))) + end + + fun delete index c = Inttab.delete index c handle Inttab.UNDEF _ => c + + val (list, c) = Inttab.foldl + (fn ((list, c), (index, v)) => ((mk_constr index v c)::list, delete index c)) + (([], c), transpose_matrix A) + val _ = if Inttab.is_empty c then () else raise Superfluous_constr_right_hand_sides + + val prog = cplex.cplexProg ("dual", cplex.cplexMinimize (vec2sum b), list, []) + in + (prog, indexof) + end + +fun cut_vector size v = + let + val count = ref 0 + fun app (v, (i, s)) = + if (!count < size) then + (count := !count +1 ; Inttab.update ((i,s),v)) + else + v + in + Inttab.foldl app (empty_vector, v) + end + +fun cut_matrix vfilter vsize m = + let + fun app (m, (i, v)) = + if (Inttab.lookup (vfilter, i) = NONE) then + m + else + case vsize of + NONE => Inttab.update ((i,v), m) + | SOME s => Inttab.update((i, cut_vector s v),m) + in + Inttab.foldl app (empty_matrix, m) + end + +fun v_elem_at v i = Inttab.lookup (v,i) +fun m_elem_at m i = Inttab.lookup (m,i) + +fun v_only_elem v = + case Inttab.min_key v of + NONE => NONE + | SOME vmin => (case Inttab.max_key v of + NONE => SOME vmin + | SOME vmax => if vmin = vmax then SOME vmin else NONE) + +fun v_fold f a v = Inttab.foldl f (a,v) + +fun m_fold f a m = Inttab.foldl f (a,m) + +end; diff -r 26fccaaf9cb4 -r 92ff7c903585 src/HOL/Matrix/cplex/fspmlp.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/cplex/fspmlp.ML Wed Jul 13 09:53:50 2005 +0200 @@ -0,0 +1,319 @@ +(* Title: HOL/Matrix/cplex/fspmlp.ML + ID: $Id$ + Author: Steven Obua +*) + +signature FSPMLP = +sig + type linprog + + val y : linprog -> cterm + val A : linprog -> cterm * cterm + val b : linprog -> cterm + val c : linprog -> cterm * cterm + val r : linprog -> cterm + val r12 : linprog -> cterm * cterm + + exception Load of string + + val load : string -> int -> bool -> linprog +end + +structure fspmlp : FSPMLP = +struct + +type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm * (cterm * cterm) + +fun y (c1, c2, c3, c4, c5, _) = c1 +fun A (c1, c2, c3, c4, c5, _) = c2 +fun b (c1, c2, c3, c4, c5, _) = c3 +fun c (c1, c2, c3, c4, c5, _) = c4 +fun r (c1, c2, c3, c4, c5, _) = c5 +fun r12 (c1, c2, c3, c4, c5, c6) = c6 + +structure CplexFloatSparseMatrixConverter = +MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder); + +datatype bound_type = LOWER | UPPER + +fun intbound_ord ((i1, b1),(i2,b2)) = + if i1 < i2 then LESS + else if i1 = i2 then + (if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER) + else GREATER + +structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord)); + +structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord); +(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *) +(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *) + +exception Internal of string; + +fun add_row_bound g dest_key row_index row_bound = + let + val x = + case VarGraph.lookup (g, dest_key) of + NONE => (NONE, Inttab.update ((row_index, (row_bound, [])), Inttab.empty)) + | SOME (sure_bound, f) => + (sure_bound, + case Inttab.lookup (f, row_index) of + NONE => Inttab.update ((row_index, (row_bound, [])), f) + | SOME _ => raise (Internal "add_row_bound")) + in + VarGraph.update ((dest_key, x), g) + end + +fun update_sure_bound g (key as (_, btype)) bound = + let + val x = + case VarGraph.lookup (g, key) of + NONE => (SOME bound, Inttab.empty) + | SOME (NONE, f) => (SOME bound, f) + | SOME (SOME old_bound, f) => + (SOME ((case btype of + UPPER => FloatArith.min + | LOWER => FloatArith.max) + old_bound bound), f) + in + VarGraph.update ((key, x), g) + end + +fun get_sure_bound g key = + case VarGraph.lookup (g, key) of + NONE => NONE + | SOME (sure_bound, _) => sure_bound + +(*fun get_row_bound g key row_index = + case VarGraph.lookup (g, key) of + NONE => NONE + | SOME (sure_bound, f) => + (case Inttab.lookup (f, row_index) of + NONE => NONE + | SOME (row_bound, _) => (sure_bound, row_bound))*) + +fun add_edge g src_key dest_key row_index coeff = + case VarGraph.lookup (g, dest_key) of + NONE => raise (Internal "add_edge: dest_key not found") + | SOME (sure_bound, f) => + (case Inttab.lookup (f, row_index) of + NONE => raise (Internal "add_edge: row_index not found") + | SOME (row_bound, sources) => + VarGraph.update ((dest_key, (sure_bound, Inttab.update ((row_index, (row_bound, (coeff, src_key) :: sources)), f))), g)) + +fun split_graph g = + let + fun split ((r1, r2), (key, (sure_bound, _))) = + case sure_bound of + NONE => (r1, r2) + | SOME bound => + (case key of + (u, UPPER) => (r1, Inttab.update ((u, bound), r2)) + | (u, LOWER) => (Inttab.update ((u, bound), r1), r2)) + in + VarGraph.foldl split ((Inttab.empty, Inttab.empty), g) + end + +fun it2list t = + let + fun tolist (l, a) = a::l + in + Inttab.foldl tolist ([], t) + end + +(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm). + If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *) +fun propagate_sure_bounds safe names g = + let + (* returns NONE if no new sure bound could be calculated, otherwise the new sure bound is returned *) + fun calc_sure_bound_from_sources g (key as (_, btype)) = + let + fun mult_upper x (lower, upper) = + if FloatArith.is_negative x then + FloatArith.mul x lower + else + FloatArith.mul x upper + + fun mult_lower x (lower, upper) = + if FloatArith.is_negative x then + FloatArith.mul x upper + else + FloatArith.mul x lower + + val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower + + fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) = + let + fun add_src_bound (sum, (coeff, src_key)) = + case sum of + NONE => NONE + | SOME x => + (case get_sure_bound g src_key of + NONE => NONE + | SOME src_sure_bound => SOME (FloatArith.add x (mult_btype src_sure_bound coeff))) + in + case Library.foldl add_src_bound (SOME row_bound, sources) of + NONE => sure_bound + | new_sure_bound as (SOME new_bound) => + (case sure_bound of + NONE => new_sure_bound + | SOME old_bound => + SOME (case btype of + UPPER => FloatArith.min old_bound new_bound + | LOWER => FloatArith.max old_bound new_bound)) + end + in + case VarGraph.lookup (g, key) of + NONE => NONE + | SOME (sure_bound, f) => + let + val x = Inttab.foldl calc_sure_bound (sure_bound, f) + in + if x = sure_bound then NONE else x + end + end + + fun propagate ((g, b), (key, _)) = + case calc_sure_bound_from_sources g key of + NONE => (g,b) + | SOME bound => (update_sure_bound g key bound, + if safe then + case get_sure_bound g key of + NONE => true + | _ => b + else + true) + + val (g, b) = VarGraph.foldl propagate ((g, false), g) + in + if b then propagate_sure_bounds safe names g else g + end + +exception Load of string; + +fun calcr safe_propagation xlen names prec A b = + let + val empty = Inttab.empty + + fun instab t i x = Inttab.update ((i,x), t) + + fun isnegstr x = String.isPrefix "-" x + fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x + + fun test_1 (lower, upper) = + if lower = upper then + (if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1 + else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1 + else 0) + else 0 + + fun calcr (g, (row_index, a)) = + let + val b = FloatSparseMatrixBuilder.v_elem_at b row_index + val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of NONE => "0" | SOME b => b) + val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) => + (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a + + fun fold_dest_nodes (g, (dest_index, dest_value)) = + let + val dest_test = test_1 dest_value + in + if dest_test = 0 then + g + else let + val (dest_key as (_, dest_btype), row_bound) = + if dest_test = ~1 then + ((dest_index, LOWER), FloatArith.neg b2) + else + ((dest_index, UPPER), b2) + + fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) = + if src_index = dest_index then g + else + let + val coeff = case dest_btype of + UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower) + | LOWER => src_value + in + if FloatArith.is_negative src_lower then + add_edge g (src_index, UPPER) dest_key row_index coeff + else + add_edge g (src_index, LOWER) dest_key row_index coeff + end + in + Library.foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a) + end + end + in + case approx_a of + [] => g + | [(u, a)] => + let + val atest = test_1 a + in + if atest = ~1 then + update_sure_bound g (u, LOWER) (FloatArith.neg b2) + else if atest = 1 then + update_sure_bound g (u, UPPER) b2 + else + g + end + | _ => Library.foldl fold_dest_nodes (g, approx_a) + end + + val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A + val g = propagate_sure_bounds safe_propagation names g + + val (r1, r2) = split_graph g + + fun add_row_entry m index value = + let + val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 value) FloatSparseMatrixBuilder.empty_spvec + in + FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m + end + + fun abs_estimate i r1 r2 = + if i = 0 then + let val e = FloatSparseMatrixBuilder.empty_spmat in (e, (e, e)) end + else + let + val index = xlen-i + val (r, (r12_1, r12_2)) = abs_estimate (i-1) r1 r2 + val b1 = case Inttab.lookup (r1, index) of NONE => raise (Load ("x-value not bounded from below: "^(names index))) | SOME x => x + val b2 = case Inttab.lookup (r2, index) of NONE => raise (Load ("x-value not bounded from above: "^(names index))) | SOME x => x + val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2) + in + (add_row_entry r index abs_max, (add_row_entry r12_1 index b1, add_row_entry r12_2 index b2)) + end + val sign = FloatSparseMatrixBuilder.sign_term + val (r, (r1, r2)) = abs_estimate xlen r1 r2 + in + (sign r, (sign r1, sign r2)) + end + +fun load filename prec safe_propagation = + let + val prog = Cplex.load_cplexFile filename + val prog = Cplex.elim_nonfree_bounds prog + val prog = Cplex.relax_strict_ineqs prog + val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog + val (r, (r1, r2)) = calcr safe_propagation xlen names prec A b + val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems" + val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b + val results = Cplex.solve dualprog + val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof + val A = FloatSparseMatrixBuilder.cut_matrix v NONE A + fun id x = x + val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v + val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b) + val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c + val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v + val A = FloatSparseMatrixBuilder.approx_matrix prec id A + val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b + val c = FloatSparseMatrixBuilder.approx_matrix prec id c + in + (y1, A, b2, c, r, (r1, r2)) + end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s))) + +end