# HG changeset patch # User haftmann # Date 1278917892 -7200 # Node ID 3489daf839d51edc4a0df05f1032ec49f1bce222 # Parent 00ff97087ab54925821f8ceb5f0c81c578a5ff4d more regular session structure diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/IsaMakefile --- a/src/HOL/IsaMakefile Fri Jul 09 17:00:42 2010 +0200 +++ b/src/HOL/IsaMakefile Mon Jul 12 08:58:12 2010 +0200 @@ -1057,10 +1057,10 @@ $(SRC)/Tools/Compute_Oracle/compute.ML Matrix/ComputeFloat.thy \ Matrix/ComputeHOL.thy Matrix/ComputeNumeral.thy Tools/float_arith.ML \ Matrix/Matrix.thy Matrix/SparseMatrix.thy Matrix/LP.thy \ - Matrix/document/root.tex Matrix/ROOT.ML Matrix/cplex/Cplex.thy \ - Matrix/cplex/CplexMatrixConverter.ML Matrix/cplex/Cplex_tools.ML \ - Matrix/cplex/FloatSparseMatrixBuilder.ML Matrix/cplex/fspmlp.ML \ - Matrix/cplex/matrixlp.ML + Matrix/document/root.tex Matrix/ROOT.ML Matrix/Cplex.thy \ + Matrix/CplexMatrixConverter.ML Matrix/Cplex_tools.ML \ + Matrix/FloatSparseMatrixBuilder.ML Matrix/fspmlp.ML \ + Matrix/matrixlp.ML @$(ISABELLE_TOOL) usedir -g true $(OUT)/HOL Matrix diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/Cplex.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/Cplex.thy Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,66 @@ +(* Title: HOL/Matrix/cplex/Cplex.thy + Author: Steven Obua +*) + +theory Cplex +imports SparseMatrix LP ComputeFloat ComputeNumeral +uses "Cplex_tools.ML" "CplexMatrixConverter.ML" "FloatSparseMatrixBuilder.ML" + "fspmlp.ML" ("matrixlp.ML") +begin + +lemma spm_mult_le_dual_prts: + assumes + "sorted_sparse_matrix A1" + "sorted_sparse_matrix A2" + "sorted_sparse_matrix c1" + "sorted_sparse_matrix c2" + "sorted_sparse_matrix y" + "sorted_sparse_matrix r1" + "sorted_sparse_matrix r2" + "sorted_spvec b" + "le_spmat [] y" + "sparse_row_matrix A1 \ A" + "A \ sparse_row_matrix A2" + "sparse_row_matrix c1 \ c" + "c \ sparse_row_matrix c2" + "sparse_row_matrix r1 \ x" + "x \ sparse_row_matrix r2" + "A * x \ sparse_row_matrix (b::('a::lattice_ring) spmat)" + shows + "c * x \ sparse_row_matrix (add_spmat (mult_spmat y b) + (let s1 = diff_spmat c1 (mult_spmat y A2); s2 = diff_spmat c2 (mult_spmat y A1) in + add_spmat (mult_spmat (pprt_spmat s2) (pprt_spmat r2)) (add_spmat (mult_spmat (pprt_spmat s1) (nprt_spmat r2)) + (add_spmat (mult_spmat (nprt_spmat s2) (pprt_spmat r1)) (mult_spmat (nprt_spmat s1) (nprt_spmat r1))))))" + apply (simp add: Let_def) + apply (insert assms) + apply (simp add: sparse_row_matrix_op_simps algebra_simps) + apply (rule mult_le_dual_prts[where A=A, simplified Let_def algebra_simps]) + apply (auto) + done + +lemma spm_mult_le_dual_prts_no_let: + assumes + "sorted_sparse_matrix A1" + "sorted_sparse_matrix A2" + "sorted_sparse_matrix c1" + "sorted_sparse_matrix c2" + "sorted_sparse_matrix y" + "sorted_sparse_matrix r1" + "sorted_sparse_matrix r2" + "sorted_spvec b" + "le_spmat [] y" + "sparse_row_matrix A1 \ A" + "A \ sparse_row_matrix A2" + "sparse_row_matrix c1 \ c" + "c \ sparse_row_matrix c2" + "sparse_row_matrix r1 \ x" + "x \ sparse_row_matrix r2" + "A * x \ sparse_row_matrix (b::('a::lattice_ring) spmat)" + shows + "c * x \ sparse_row_matrix (add_spmat (mult_spmat y b) + (mult_est_spmat r1 r2 (diff_spmat c1 (mult_spmat y A2)) (diff_spmat c2 (mult_spmat y A1))))" + by (simp add: assms mult_est_spmat_def spm_mult_le_dual_prts[where A=A, simplified Let_def]) + +use "matrixlp.ML" + +end diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/CplexMatrixConverter.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/CplexMatrixConverter.ML Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,128 @@ +(* Title: HOL/Matrix/cplex/CplexMatrixConverter.ML + 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; + +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) = fold (fn t => fn vec => setprod vec true t) ts empty_vector + | 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 (name, value) v = matrix_builder.set_elem v (name2index name) value + in + (opt, fold setv entries (matrix_builder.empty_vector)) + 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 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/Cplex_tools.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/Cplex_tools.ML Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,1225 @@ +(* Title: HOL/Matrix/cplex/Cplex_tools.ML + 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) + + datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK + + 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 get_solver : unit -> cplexSolver + val set_solver : cplexSolver -> unit + val solve : cplexProg -> cplexResult +end; + +structure Cplex : CPLEX = +struct + +datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK + +val cplexsolver = Unsynchronized.ref SOLVER_DEFAULT; +fun get_solver () = !cplexsolver; +fun set_solver s = (cplexsolver := s); + +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 = Unsynchronized.ref true + val rest = Unsynchronized.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 + (case TextIO.inputLine f of + NONE => NONE + | SOME s => + 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) + + 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 = Unsynchronized.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 = fold (fn x => fn table => merge table (f x)) ts Symtab.empty + +fun diff a b = Symtab.fold (Symtab.delete_safe o fst) b a + +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.fold + (fn (k, v) => cons (make_pos_constr k)) + positive_vars []) + val bound_constrs = map (pair NONE) + (maps bound2constr bounds) + val constraints' = constraints @ pos_constrs @ bound_constrs + val bounds' = rev (Symtab.fold (fn (v, _) => cons (make_free_bound v)) 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 = Unsynchronized.ref [] + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + SOME u + ) + end + else + (case TextIO.inputLine f of + NONE => NONE + | SOME s => (rest := tokenize s; readToken_helper())) + + 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 AList.lookup (op =) header "STATUS" of + SOME "INFEASIBLE" => Infeasible + | SOME "UNBOUNDED" => Unbounded + | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) 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 = Unsynchronized.ref [] + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + SOME u + ) + end + else + (case TextIO.inputLine f of + NONE => NONE + | SOME s => (rest := tokenize s; readToken_helper())) + + 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 AList.lookup (op =) header "STATUS" of + SOME "INFEASIBLE" => Infeasible + | SOME "NONOPTIMAL" => Unbounded + | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) 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.implode (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 "GLPK_PATH" + val cplex = if cplex_path = "" then "glpsol" else cplex_path + val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname) + val answer = #1 (bash_output 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 "CPLEX_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 (bash 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 get_solver () of + SOLVER_DEFAULT => + (case getenv "LP_SOLVER" of + "CPLEX" => solve_cplex prog + | "GLPK" => solve_glpk prog + | _ => raise (Execute ("LP_SOLVER must be set to CPLEX or to GLPK"))) + | SOLVER_CPLEX => solve_cplex prog + | SOLVER_GLPK => solve_glpk prog + +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 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/FloatSparseMatrixBuilder.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/FloatSparseMatrixBuilder.ML Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,284 @@ +(* Title: HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML + Author: Steven Obua +*) + +signature FLOAT_SPARSE_MATIRX_BUILDER = +sig + include MATRIX_BUILDER + + structure cplex : CPLEX + + type float = Float.float + val approx_value : int -> (float -> float) -> string -> term * term + val approx_vector : int -> (float -> float) -> vector -> term * term + val approx_matrix : int -> (float -> float) -> matrix -> term * term + + val mk_spvec_entry : int -> float -> term + val mk_spvec_entry' : int -> term -> term + val mk_spmat_entry : int -> term -> term + val spvecT: typ + val spmatT: typ + + 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 : (int * string -> 'a -> 'a) -> vector -> 'a -> 'a + val m_fold : (int * vector -> 'a -> 'a) -> matrix -> 'a -> 'a + + val transpose_matrix : matrix -> matrix + + val cut_vector : int -> vector -> vector + val cut_matrix : vector -> int option -> matrix -> matrix + + val delete_matrix : int list -> matrix -> matrix + val cut_matrix' : int list -> matrix -> matrix + val delete_vector : int list -> vector -> vector + val cut_vector' : int list -> vector -> vector + + val indices_of_matrix : matrix -> int list + val indices_of_vector : vector -> int list + + (* 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) +end; + +structure FloatSparseMatrixBuilder : FLOAT_SPARSE_MATIRX_BUILDER = +struct + +type float = Float.float +structure Inttab = Table(type key = int val ord = rev_order o int_ord); + +type vector = string Inttab.table +type matrix = vector Inttab.table + +val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT); +val spvecT = HOLogic.listT spvec_elemT; +val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT); +val spmatT = HOLogic.listT spmat_elemT; + +fun approx_value prec f = + FloatArith.approx_float prec (fn (x, y) => (f x, f y)); + +fun mk_spvec_entry i f = + HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, FloatArith.mk_float f); + +fun mk_spvec_entry' i x = + HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, x); + +fun mk_spmat_entry i e = + HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, e); + +fun approx_vector prec pprt vector = + let + fun app (index, s) (lower, upper) = + let + val (flower, fupper) = approx_value prec pprt s + val index = HOLogic.mk_number HOLogic.natT index + val elower = HOLogic.mk_prod (index, flower) + val eupper = HOLogic.mk_prod (index, fupper) + in (elower :: lower, eupper :: upper) end; + in + pairself (HOLogic.mk_list spvec_elemT) (Inttab.fold app vector ([], [])) + end; + +fun approx_matrix prec pprt vector = + let + fun app (index, v) (lower, upper) = + let + val (flower, fupper) = approx_vector prec pprt v + val index = HOLogic.mk_number HOLogic.natT index + val elower = HOLogic.mk_prod (index, flower) + val eupper = HOLogic.mk_prod (index, fupper) + in (elower :: lower, eupper :: upper) end; + in + pairself (HOLogic.mk_list spmat_elemT) (Inttab.fold app vector ([], [])) + end; + +exception Nat_expected of int; + +val zero_interval = approx_value 1 I "0" + +fun set_elem vector index str = + if index < 0 then + raise (Nat_expected index) + else if (approx_value 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 j (i, s) = + Inttab.map_default (i, Inttab.empty) (Inttab.update (j, s)); + fun updm (j, v) = Inttab.fold (upd j) v; + in Inttab.fold updm matrix empty_matrix end; + +exception No_name of string; + +exception Superfluous_constr_right_hand_sides + +fun cplexProg c A b = + let + val ytable = Unsynchronized.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 _ = Unsynchronized.change ytable (Inttab.update (i, s)) + 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.fold (fn (index, s) => fn list => (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.fold + (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) + A ([], b) + 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.fold (fn (i, y) => fn l => (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.fold (fn (index, s) => fn list => (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.fold + (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) + (transpose_matrix A) ([], c) + 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 = Unsynchronized.ref 0; + fun app (i, s) = if (!count < size) then + (count := !count +1 ; Inttab.update (i, s)) + else I + in + Inttab.fold app v empty_vector + end + +fun cut_matrix vfilter vsize m = + let + fun app (i, v) = + if is_none (Inttab.lookup vfilter i) then I + else case vsize + of NONE => Inttab.update (i, v) + | SOME s => Inttab.update (i, cut_vector s v) + in Inttab.fold app m empty_matrix 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 = Inttab.fold f; +fun m_fold f = Inttab.fold f; + +fun indices_of_vector v = Inttab.keys v +fun indices_of_matrix m = Inttab.keys m +fun delete_vector indices v = fold Inttab.delete indices v +fun delete_matrix indices m = fold Inttab.delete indices m +fun cut_matrix' indices m = fold (fn i => fn m => (case Inttab.lookup m i of NONE => m | SOME v => Inttab.update (i, v) m)) indices Inttab.empty +fun cut_vector' indices v = fold (fn i => fn v => (case Inttab.lookup v i of NONE => v | SOME x => Inttab.update (i, x) v)) indices Inttab.empty + + + +end; diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/ROOT.ML --- a/src/HOL/Matrix/ROOT.ML Fri Jul 09 17:00:42 2010 +0200 +++ b/src/HOL/Matrix/ROOT.ML Mon Jul 12 08:58:12 2010 +0200 @@ -1,1 +1,2 @@ -use_thys ["SparseMatrix", "cplex/Cplex"]; + +use_thy "Cplex"; diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/Cplex.thy --- a/src/HOL/Matrix/cplex/Cplex.thy Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,66 +0,0 @@ -(* Title: HOL/Matrix/cplex/Cplex.thy - Author: Steven Obua -*) - -theory Cplex -imports SparseMatrix LP ComputeFloat ComputeNumeral -uses "Cplex_tools.ML" "CplexMatrixConverter.ML" "FloatSparseMatrixBuilder.ML" - "fspmlp.ML" ("matrixlp.ML") -begin - -lemma spm_mult_le_dual_prts: - assumes - "sorted_sparse_matrix A1" - "sorted_sparse_matrix A2" - "sorted_sparse_matrix c1" - "sorted_sparse_matrix c2" - "sorted_sparse_matrix y" - "sorted_sparse_matrix r1" - "sorted_sparse_matrix r2" - "sorted_spvec b" - "le_spmat [] y" - "sparse_row_matrix A1 \ A" - "A \ sparse_row_matrix A2" - "sparse_row_matrix c1 \ c" - "c \ sparse_row_matrix c2" - "sparse_row_matrix r1 \ x" - "x \ sparse_row_matrix r2" - "A * x \ sparse_row_matrix (b::('a::lattice_ring) spmat)" - shows - "c * x \ sparse_row_matrix (add_spmat (mult_spmat y b) - (let s1 = diff_spmat c1 (mult_spmat y A2); s2 = diff_spmat c2 (mult_spmat y A1) in - add_spmat (mult_spmat (pprt_spmat s2) (pprt_spmat r2)) (add_spmat (mult_spmat (pprt_spmat s1) (nprt_spmat r2)) - (add_spmat (mult_spmat (nprt_spmat s2) (pprt_spmat r1)) (mult_spmat (nprt_spmat s1) (nprt_spmat r1))))))" - apply (simp add: Let_def) - apply (insert assms) - apply (simp add: sparse_row_matrix_op_simps algebra_simps) - apply (rule mult_le_dual_prts[where A=A, simplified Let_def algebra_simps]) - apply (auto) - done - -lemma spm_mult_le_dual_prts_no_let: - assumes - "sorted_sparse_matrix A1" - "sorted_sparse_matrix A2" - "sorted_sparse_matrix c1" - "sorted_sparse_matrix c2" - "sorted_sparse_matrix y" - "sorted_sparse_matrix r1" - "sorted_sparse_matrix r2" - "sorted_spvec b" - "le_spmat [] y" - "sparse_row_matrix A1 \ A" - "A \ sparse_row_matrix A2" - "sparse_row_matrix c1 \ c" - "c \ sparse_row_matrix c2" - "sparse_row_matrix r1 \ x" - "x \ sparse_row_matrix r2" - "A * x \ sparse_row_matrix (b::('a::lattice_ring) spmat)" - shows - "c * x \ sparse_row_matrix (add_spmat (mult_spmat y b) - (mult_est_spmat r1 r2 (diff_spmat c1 (mult_spmat y A2)) (diff_spmat c2 (mult_spmat y A1))))" - by (simp add: assms mult_est_spmat_def spm_mult_le_dual_prts[where A=A, simplified Let_def]) - -use "matrixlp.ML" - -end diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/CplexMatrixConverter.ML --- a/src/HOL/Matrix/cplex/CplexMatrixConverter.ML Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,128 +0,0 @@ -(* Title: HOL/Matrix/cplex/CplexMatrixConverter.ML - 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; - -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) = fold (fn t => fn vec => setprod vec true t) ts empty_vector - | 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 (name, value) v = matrix_builder.set_elem v (name2index name) value - in - (opt, fold setv entries (matrix_builder.empty_vector)) - 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 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/Cplex_tools.ML --- a/src/HOL/Matrix/cplex/Cplex_tools.ML Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1225 +0,0 @@ -(* Title: HOL/Matrix/cplex/Cplex_tools.ML - 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) - - datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK - - 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 get_solver : unit -> cplexSolver - val set_solver : cplexSolver -> unit - val solve : cplexProg -> cplexResult -end; - -structure Cplex : CPLEX = -struct - -datatype cplexSolver = SOLVER_DEFAULT | SOLVER_CPLEX | SOLVER_GLPK - -val cplexsolver = Unsynchronized.ref SOLVER_DEFAULT; -fun get_solver () = !cplexsolver; -fun set_solver s = (cplexsolver := s); - -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 = Unsynchronized.ref true - val rest = Unsynchronized.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 - (case TextIO.inputLine f of - NONE => NONE - | SOME s => - 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) - - 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 = Unsynchronized.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 = fold (fn x => fn table => merge table (f x)) ts Symtab.empty - -fun diff a b = Symtab.fold (Symtab.delete_safe o fst) b a - -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.fold - (fn (k, v) => cons (make_pos_constr k)) - positive_vars []) - val bound_constrs = map (pair NONE) - (maps bound2constr bounds) - val constraints' = constraints @ pos_constrs @ bound_constrs - val bounds' = rev (Symtab.fold (fn (v, _) => cons (make_free_bound v)) 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 = Unsynchronized.ref [] - - fun readToken_helper () = - if length (!rest) > 0 then - let val u = hd (!rest) in - ( - rest := tl (!rest); - SOME u - ) - end - else - (case TextIO.inputLine f of - NONE => NONE - | SOME s => (rest := tokenize s; readToken_helper())) - - 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 AList.lookup (op =) header "STATUS" of - SOME "INFEASIBLE" => Infeasible - | SOME "UNBOUNDED" => Unbounded - | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) 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 = Unsynchronized.ref [] - - fun readToken_helper () = - if length (!rest) > 0 then - let val u = hd (!rest) in - ( - rest := tl (!rest); - SOME u - ) - end - else - (case TextIO.inputLine f of - NONE => NONE - | SOME s => (rest := tokenize s; readToken_helper())) - - 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 AList.lookup (op =) header "STATUS" of - SOME "INFEASIBLE" => Infeasible - | SOME "NONOPTIMAL" => Unbounded - | SOME "OPTIMAL" => Optimal (the (AList.lookup (op =) 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.implode (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 "GLPK_PATH" - val cplex = if cplex_path = "" then "glpsol" else cplex_path - val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname) - val answer = #1 (bash_output 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 "CPLEX_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 (bash 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 get_solver () of - SOLVER_DEFAULT => - (case getenv "LP_SOLVER" of - "CPLEX" => solve_cplex prog - | "GLPK" => solve_glpk prog - | _ => raise (Execute ("LP_SOLVER must be set to CPLEX or to GLPK"))) - | SOLVER_CPLEX => solve_cplex prog - | SOLVER_GLPK => solve_glpk prog - -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 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML --- a/src/HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,284 +0,0 @@ -(* Title: HOL/Matrix/cplex/FloatSparseMatrixBuilder.ML - Author: Steven Obua -*) - -signature FLOAT_SPARSE_MATIRX_BUILDER = -sig - include MATRIX_BUILDER - - structure cplex : CPLEX - - type float = Float.float - val approx_value : int -> (float -> float) -> string -> term * term - val approx_vector : int -> (float -> float) -> vector -> term * term - val approx_matrix : int -> (float -> float) -> matrix -> term * term - - val mk_spvec_entry : int -> float -> term - val mk_spvec_entry' : int -> term -> term - val mk_spmat_entry : int -> term -> term - val spvecT: typ - val spmatT: typ - - 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 : (int * string -> 'a -> 'a) -> vector -> 'a -> 'a - val m_fold : (int * vector -> 'a -> 'a) -> matrix -> 'a -> 'a - - val transpose_matrix : matrix -> matrix - - val cut_vector : int -> vector -> vector - val cut_matrix : vector -> int option -> matrix -> matrix - - val delete_matrix : int list -> matrix -> matrix - val cut_matrix' : int list -> matrix -> matrix - val delete_vector : int list -> vector -> vector - val cut_vector' : int list -> vector -> vector - - val indices_of_matrix : matrix -> int list - val indices_of_vector : vector -> int list - - (* 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) -end; - -structure FloatSparseMatrixBuilder : FLOAT_SPARSE_MATIRX_BUILDER = -struct - -type float = Float.float -structure Inttab = Table(type key = int val ord = rev_order o int_ord); - -type vector = string Inttab.table -type matrix = vector Inttab.table - -val spvec_elemT = HOLogic.mk_prodT (HOLogic.natT, HOLogic.realT); -val spvecT = HOLogic.listT spvec_elemT; -val spmat_elemT = HOLogic.mk_prodT (HOLogic.natT, spvecT); -val spmatT = HOLogic.listT spmat_elemT; - -fun approx_value prec f = - FloatArith.approx_float prec (fn (x, y) => (f x, f y)); - -fun mk_spvec_entry i f = - HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, FloatArith.mk_float f); - -fun mk_spvec_entry' i x = - HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, x); - -fun mk_spmat_entry i e = - HOLogic.mk_prod (HOLogic.mk_number HOLogic.natT i, e); - -fun approx_vector prec pprt vector = - let - fun app (index, s) (lower, upper) = - let - val (flower, fupper) = approx_value prec pprt s - val index = HOLogic.mk_number HOLogic.natT index - val elower = HOLogic.mk_prod (index, flower) - val eupper = HOLogic.mk_prod (index, fupper) - in (elower :: lower, eupper :: upper) end; - in - pairself (HOLogic.mk_list spvec_elemT) (Inttab.fold app vector ([], [])) - end; - -fun approx_matrix prec pprt vector = - let - fun app (index, v) (lower, upper) = - let - val (flower, fupper) = approx_vector prec pprt v - val index = HOLogic.mk_number HOLogic.natT index - val elower = HOLogic.mk_prod (index, flower) - val eupper = HOLogic.mk_prod (index, fupper) - in (elower :: lower, eupper :: upper) end; - in - pairself (HOLogic.mk_list spmat_elemT) (Inttab.fold app vector ([], [])) - end; - -exception Nat_expected of int; - -val zero_interval = approx_value 1 I "0" - -fun set_elem vector index str = - if index < 0 then - raise (Nat_expected index) - else if (approx_value 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 j (i, s) = - Inttab.map_default (i, Inttab.empty) (Inttab.update (j, s)); - fun updm (j, v) = Inttab.fold (upd j) v; - in Inttab.fold updm matrix empty_matrix end; - -exception No_name of string; - -exception Superfluous_constr_right_hand_sides - -fun cplexProg c A b = - let - val ytable = Unsynchronized.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 _ = Unsynchronized.change ytable (Inttab.update (i, s)) - 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.fold (fn (index, s) => fn list => (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.fold - (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) - A ([], b) - 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.fold (fn (i, y) => fn l => (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.fold (fn (index, s) => fn list => (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.fold - (fn (index, v) => fn (list, c) => ((mk_constr index v c)::list, delete index c)) - (transpose_matrix A) ([], c) - 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 = Unsynchronized.ref 0; - fun app (i, s) = if (!count < size) then - (count := !count +1 ; Inttab.update (i, s)) - else I - in - Inttab.fold app v empty_vector - end - -fun cut_matrix vfilter vsize m = - let - fun app (i, v) = - if is_none (Inttab.lookup vfilter i) then I - else case vsize - of NONE => Inttab.update (i, v) - | SOME s => Inttab.update (i, cut_vector s v) - in Inttab.fold app m empty_matrix 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 = Inttab.fold f; -fun m_fold f = Inttab.fold f; - -fun indices_of_vector v = Inttab.keys v -fun indices_of_matrix m = Inttab.keys m -fun delete_vector indices v = fold Inttab.delete indices v -fun delete_matrix indices m = fold Inttab.delete indices m -fun cut_matrix' indices m = fold (fn i => fn m => (case Inttab.lookup m i of NONE => m | SOME v => Inttab.update (i, v) m)) indices Inttab.empty -fun cut_vector' indices v = fold (fn i => fn v => (case Inttab.lookup v i of NONE => v | SOME x => Inttab.update (i, x) v)) indices Inttab.empty - - - -end; diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/fspmlp.ML --- a/src/HOL/Matrix/cplex/fspmlp.ML Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,322 +0,0 @@ -(* Title: HOL/Matrix/cplex/fspmlp.ML - Author: Steven Obua -*) - -signature FSPMLP = -sig - type linprog - type vector = FloatSparseMatrixBuilder.vector - type matrix = FloatSparseMatrixBuilder.matrix - - val y : linprog -> term - val A : linprog -> term * term - val b : linprog -> term - val c : linprog -> term * term - val r12 : linprog -> term * term - - exception Load of string - - val load : string -> int -> bool -> linprog -end - -structure Fspmlp : FSPMLP = -struct - -type vector = FloatSparseMatrixBuilder.vector -type matrix = FloatSparseMatrixBuilder.matrix - -type linprog = term * (term * term) * term * (term * term) * (term * term) - -fun y (c1, c2, c3, c4, _) = c1 -fun A (c1, c2, c3, c4, _) = c2 -fun b (c1, c2, c3, c4, _) = c3 -fun c (c1, c2, c3, c4, _) = c4 -fun r12 (c1, c2, c3, c4, c6) = c6 - -structure CplexFloatSparseMatrixConverter = -MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder); - -datatype bound_type = LOWER | UPPER - -fun intbound_ord ((i1: int, 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 = Table(type key = int val ord = (rev_order o int_ord)); - -structure VarGraph = Table(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 => Float.min - | LOWER => Float.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 (key, (sure_bound, _)) (r1, r2) = 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.fold split g (Inttab.empty, Inttab.empty) end - -fun it2list t = Inttab.fold cons t []; - -(* 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 Float.sign x = LESS then - Float.mult x lower - else - Float.mult x upper - - fun mult_lower x (lower, upper) = - if Float.sign x = LESS then - Float.mult x upper - else - Float.mult x lower - - val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower - - fun calc_sure_bound (row_index, (row_bound, sources)) sure_bound = - let - fun add_src_bound (coeff, src_key) sum = - case sum of - NONE => NONE - | SOME x => - (case get_sure_bound g src_key of - NONE => NONE - | SOME src_sure_bound => SOME (Float.add x (mult_btype src_sure_bound coeff))) - in - case fold add_src_bound sources (SOME row_bound) 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 => Float.min old_bound new_bound - | LOWER => Float.max old_bound new_bound)) - end - in - case VarGraph.lookup g key of - NONE => NONE - | SOME (sure_bound, f) => - let - val x = Inttab.fold calc_sure_bound f sure_bound - in - if x = sure_bound then NONE else x - end - end - - fun propagate (key, _) (g, b) = - 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.fold propagate g (g, false) - in - if b then propagate_sure_bounds safe names g else g - end - -exception Load of string; - -val empty_spvec = @{term "Nil :: real spvec"}; -fun cons_spvec x xs = @{term "Cons :: nat * real => real spvec => real spvec"} $ x $ xs; -val empty_spmat = @{term "Nil :: real spmat"}; -fun cons_spmat x xs = @{term "Cons :: nat * real spvec => real spmat => real spmat"} $ x $ xs; - -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 Float.eq (lower, (~1, 0)) then ~1 - else if Float.eq (lower, (1, 0)) then 1 - else 0) - else 0 - - fun calcr (row_index, a) g = - let - val b = FloatSparseMatrixBuilder.v_elem_at b row_index - val (_, b2) = FloatArith.approx_decstr_by_bin prec (case b of NONE => "0" | SOME b => b) - val approx_a = FloatSparseMatrixBuilder.v_fold (fn (i, s) => fn l => - (i, FloatArith.approx_decstr_by_bin prec s)::l) a [] - - fun fold_dest_nodes (dest_index, dest_value) g = - 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), Float.neg b2) - else - ((dest_index, UPPER), b2) - - fun fold_src_nodes (src_index, src_value as (src_lower, src_upper)) g = - if src_index = dest_index then g - else - let - val coeff = case dest_btype of - UPPER => (Float.neg src_upper, Float.neg src_lower) - | LOWER => src_value - in - if Float.sign src_lower = LESS 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 - fold fold_src_nodes approx_a (add_row_bound g dest_key row_index row_bound) - 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) (Float.neg b2) - else if atest = 1 then - update_sure_bound g (u, UPPER) b2 - else - g - end - | _ => fold fold_dest_nodes approx_a g - end - - val g = FloatSparseMatrixBuilder.m_fold calcr A VarGraph.empty - - val g = propagate_sure_bounds safe_propagation names g - - val (r1, r2) = split_graph g - - fun add_row_entry m index f vname value = - let - val v = (case value of - SOME value => FloatSparseMatrixBuilder.mk_spvec_entry 0 value - | NONE => FloatSparseMatrixBuilder.mk_spvec_entry' 0 (f $ (Var ((vname,0), HOLogic.realT)))) - val vec = cons_spvec v empty_spvec - in - cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m - end - - fun abs_estimate i r1 r2 = - if i = 0 then - let val e = empty_spmat in (e, e) end - else - let - val index = xlen-i - val (r12_1, r12_2) = abs_estimate (i-1) r1 r2 - val b1 = Inttab.lookup r1 index - val b2 = Inttab.lookup r2 index - in - (add_row_entry r12_1 index @{term "lbound :: real => real"} ((names index)^"l") b1, - add_row_entry r12_2 index @{term "ubound :: real => real"} ((names index)^"u") b2) - end - - val (r1, r2) = abs_estimate xlen r1 r2 - - in - (r1, 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 (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 Float.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, (r1, r2)) - end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s))) - -end diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/cplex/matrixlp.ML --- a/src/HOL/Matrix/cplex/matrixlp.ML Fri Jul 09 17:00:42 2010 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,100 +0,0 @@ -(* Title: HOL/Matrix/cplex/matrixlp.ML - Author: Steven Obua -*) - -signature MATRIX_LP = -sig - val lp_dual_estimate_prt : string -> int -> thm - val lp_dual_estimate_prt_primitive : - cterm * (cterm * cterm) * (cterm * cterm) * cterm * (cterm * cterm) -> thm - val matrix_compute : cterm -> thm - val matrix_simplify : thm -> thm - val prove_bound : string -> int -> thm - val float2real : string * string -> Real.real -end - -structure MatrixLP : MATRIX_LP = -struct - -fun inst_real thm = - let val certT = ctyp_of (Thm.theory_of_thm thm) in - Drule.export_without_context (Thm.instantiate - ([(certT (TVar (hd (OldTerm.term_tvars (prop_of thm)))), certT HOLogic.realT)], []) thm) - end - -local - -val cert = cterm_of @{theory} - -in - -fun lp_dual_estimate_prt_primitive (y, (A1, A2), (c1, c2), b, (r1, r2)) = - let - val th = inst_real @{thm "spm_mult_le_dual_prts_no_let"} - fun var s x = (cert (Var ((s,0), FloatSparseMatrixBuilder.spmatT)), x) - val th = Thm.instantiate ([], [var "A1" A1, var "A2" A2, var "y" y, var "c1" c1, var "c2" c2, - var "r1" r1, var "r2" r2, var "b" b]) th - in th end - -fun lp_dual_estimate_prt lptfile prec = - let - val certificate = - let - open Fspmlp - val l = load lptfile prec false - in - (y l |> cert, A l |> pairself cert, c l |> pairself cert, b l |> cert, r12 l |> pairself cert) - end - in - lp_dual_estimate_prt_primitive certificate - end - -end - -fun prep ths = ComputeHOL.prep_thms ths - -fun inst_tvar ty thm = - let - val certT = Thm.ctyp_of (Thm.theory_of_thm thm); - val ord = prod_ord (prod_ord string_ord int_ord) (list_ord string_ord) - val v = TVar (hd (sort ord (OldTerm.term_tvars (prop_of thm)))) - in - Drule.export_without_context (Thm.instantiate ([(certT v, certT ty)], []) thm) - end - -fun inst_tvars [] thms = thms - | inst_tvars (ty::tys) thms = inst_tvars tys (map (inst_tvar ty) thms) - -local - val ths = ComputeHOL.prep_thms @{thms "ComputeHOL.compute_list_case" "ComputeHOL.compute_let" - "ComputeHOL.compute_if" "ComputeFloat.arith" "SparseMatrix.sparse_row_matrix_arith_simps" - "ComputeHOL.compute_bool" "ComputeHOL.compute_pair" - "SparseMatrix.sorted_sp_simps" "ComputeNumeral.number_norm" - "ComputeNumeral.natnorm"}; - val computer = PCompute.make Compute.SML @{theory} ths [] -in - -fun matrix_compute c = hd (PCompute.rewrite computer [c]) - -end - -fun matrix_simplify th = - let - val simp_th = matrix_compute (cprop_of th) - val th = Thm.strip_shyps (Thm.equal_elim simp_th th) - fun removeTrue th = removeTrue (Thm.implies_elim th TrueI) handle _ => th (* FIXME avoid handle _ *) - in - removeTrue th - end - -fun prove_bound lptfile prec = - let - val th = lp_dual_estimate_prt lptfile prec - in - matrix_simplify th - end - -val realFromStr = the o Real.fromString; -fun float2real (x, y) = realFromStr x * Math.pow (2.0, realFromStr y); - -end diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/fspmlp.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/fspmlp.ML Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,322 @@ +(* Title: HOL/Matrix/cplex/fspmlp.ML + Author: Steven Obua +*) + +signature FSPMLP = +sig + type linprog + type vector = FloatSparseMatrixBuilder.vector + type matrix = FloatSparseMatrixBuilder.matrix + + val y : linprog -> term + val A : linprog -> term * term + val b : linprog -> term + val c : linprog -> term * term + val r12 : linprog -> term * term + + exception Load of string + + val load : string -> int -> bool -> linprog +end + +structure Fspmlp : FSPMLP = +struct + +type vector = FloatSparseMatrixBuilder.vector +type matrix = FloatSparseMatrixBuilder.matrix + +type linprog = term * (term * term) * term * (term * term) * (term * term) + +fun y (c1, c2, c3, c4, _) = c1 +fun A (c1, c2, c3, c4, _) = c2 +fun b (c1, c2, c3, c4, _) = c3 +fun c (c1, c2, c3, c4, _) = c4 +fun r12 (c1, c2, c3, c4, c6) = c6 + +structure CplexFloatSparseMatrixConverter = +MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder); + +datatype bound_type = LOWER | UPPER + +fun intbound_ord ((i1: int, 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 = Table(type key = int val ord = (rev_order o int_ord)); + +structure VarGraph = Table(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 => Float.min + | LOWER => Float.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 (key, (sure_bound, _)) (r1, r2) = 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.fold split g (Inttab.empty, Inttab.empty) end + +fun it2list t = Inttab.fold cons t []; + +(* 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 Float.sign x = LESS then + Float.mult x lower + else + Float.mult x upper + + fun mult_lower x (lower, upper) = + if Float.sign x = LESS then + Float.mult x upper + else + Float.mult x lower + + val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower + + fun calc_sure_bound (row_index, (row_bound, sources)) sure_bound = + let + fun add_src_bound (coeff, src_key) sum = + case sum of + NONE => NONE + | SOME x => + (case get_sure_bound g src_key of + NONE => NONE + | SOME src_sure_bound => SOME (Float.add x (mult_btype src_sure_bound coeff))) + in + case fold add_src_bound sources (SOME row_bound) 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 => Float.min old_bound new_bound + | LOWER => Float.max old_bound new_bound)) + end + in + case VarGraph.lookup g key of + NONE => NONE + | SOME (sure_bound, f) => + let + val x = Inttab.fold calc_sure_bound f sure_bound + in + if x = sure_bound then NONE else x + end + end + + fun propagate (key, _) (g, b) = + 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.fold propagate g (g, false) + in + if b then propagate_sure_bounds safe names g else g + end + +exception Load of string; + +val empty_spvec = @{term "Nil :: real spvec"}; +fun cons_spvec x xs = @{term "Cons :: nat * real => real spvec => real spvec"} $ x $ xs; +val empty_spmat = @{term "Nil :: real spmat"}; +fun cons_spmat x xs = @{term "Cons :: nat * real spvec => real spmat => real spmat"} $ x $ xs; + +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 Float.eq (lower, (~1, 0)) then ~1 + else if Float.eq (lower, (1, 0)) then 1 + else 0) + else 0 + + fun calcr (row_index, a) g = + let + val b = FloatSparseMatrixBuilder.v_elem_at b row_index + val (_, b2) = FloatArith.approx_decstr_by_bin prec (case b of NONE => "0" | SOME b => b) + val approx_a = FloatSparseMatrixBuilder.v_fold (fn (i, s) => fn l => + (i, FloatArith.approx_decstr_by_bin prec s)::l) a [] + + fun fold_dest_nodes (dest_index, dest_value) g = + 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), Float.neg b2) + else + ((dest_index, UPPER), b2) + + fun fold_src_nodes (src_index, src_value as (src_lower, src_upper)) g = + if src_index = dest_index then g + else + let + val coeff = case dest_btype of + UPPER => (Float.neg src_upper, Float.neg src_lower) + | LOWER => src_value + in + if Float.sign src_lower = LESS 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 + fold fold_src_nodes approx_a (add_row_bound g dest_key row_index row_bound) + 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) (Float.neg b2) + else if atest = 1 then + update_sure_bound g (u, UPPER) b2 + else + g + end + | _ => fold fold_dest_nodes approx_a g + end + + val g = FloatSparseMatrixBuilder.m_fold calcr A VarGraph.empty + + val g = propagate_sure_bounds safe_propagation names g + + val (r1, r2) = split_graph g + + fun add_row_entry m index f vname value = + let + val v = (case value of + SOME value => FloatSparseMatrixBuilder.mk_spvec_entry 0 value + | NONE => FloatSparseMatrixBuilder.mk_spvec_entry' 0 (f $ (Var ((vname,0), HOLogic.realT)))) + val vec = cons_spvec v empty_spvec + in + cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) m + end + + fun abs_estimate i r1 r2 = + if i = 0 then + let val e = empty_spmat in (e, e) end + else + let + val index = xlen-i + val (r12_1, r12_2) = abs_estimate (i-1) r1 r2 + val b1 = Inttab.lookup r1 index + val b2 = Inttab.lookup r2 index + in + (add_row_entry r12_1 index @{term "lbound :: real => real"} ((names index)^"l") b1, + add_row_entry r12_2 index @{term "ubound :: real => real"} ((names index)^"u") b2) + end + + val (r1, r2) = abs_estimate xlen r1 r2 + + in + (r1, 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 (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 Float.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, (r1, r2)) + end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s))) + +end diff -r 00ff97087ab5 -r 3489daf839d5 src/HOL/Matrix/matrixlp.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/matrixlp.ML Mon Jul 12 08:58:12 2010 +0200 @@ -0,0 +1,100 @@ +(* Title: HOL/Matrix/cplex/matrixlp.ML + Author: Steven Obua +*) + +signature MATRIX_LP = +sig + val lp_dual_estimate_prt : string -> int -> thm + val lp_dual_estimate_prt_primitive : + cterm * (cterm * cterm) * (cterm * cterm) * cterm * (cterm * cterm) -> thm + val matrix_compute : cterm -> thm + val matrix_simplify : thm -> thm + val prove_bound : string -> int -> thm + val float2real : string * string -> Real.real +end + +structure MatrixLP : MATRIX_LP = +struct + +fun inst_real thm = + let val certT = ctyp_of (Thm.theory_of_thm thm) in + Drule.export_without_context (Thm.instantiate + ([(certT (TVar (hd (OldTerm.term_tvars (prop_of thm)))), certT HOLogic.realT)], []) thm) + end + +local + +val cert = cterm_of @{theory} + +in + +fun lp_dual_estimate_prt_primitive (y, (A1, A2), (c1, c2), b, (r1, r2)) = + let + val th = inst_real @{thm "spm_mult_le_dual_prts_no_let"} + fun var s x = (cert (Var ((s,0), FloatSparseMatrixBuilder.spmatT)), x) + val th = Thm.instantiate ([], [var "A1" A1, var "A2" A2, var "y" y, var "c1" c1, var "c2" c2, + var "r1" r1, var "r2" r2, var "b" b]) th + in th end + +fun lp_dual_estimate_prt lptfile prec = + let + val certificate = + let + open Fspmlp + val l = load lptfile prec false + in + (y l |> cert, A l |> pairself cert, c l |> pairself cert, b l |> cert, r12 l |> pairself cert) + end + in + lp_dual_estimate_prt_primitive certificate + end + +end + +fun prep ths = ComputeHOL.prep_thms ths + +fun inst_tvar ty thm = + let + val certT = Thm.ctyp_of (Thm.theory_of_thm thm); + val ord = prod_ord (prod_ord string_ord int_ord) (list_ord string_ord) + val v = TVar (hd (sort ord (OldTerm.term_tvars (prop_of thm)))) + in + Drule.export_without_context (Thm.instantiate ([(certT v, certT ty)], []) thm) + end + +fun inst_tvars [] thms = thms + | inst_tvars (ty::tys) thms = inst_tvars tys (map (inst_tvar ty) thms) + +local + val ths = ComputeHOL.prep_thms @{thms "ComputeHOL.compute_list_case" "ComputeHOL.compute_let" + "ComputeHOL.compute_if" "ComputeFloat.arith" "SparseMatrix.sparse_row_matrix_arith_simps" + "ComputeHOL.compute_bool" "ComputeHOL.compute_pair" + "SparseMatrix.sorted_sp_simps" "ComputeNumeral.number_norm" + "ComputeNumeral.natnorm"}; + val computer = PCompute.make Compute.SML @{theory} ths [] +in + +fun matrix_compute c = hd (PCompute.rewrite computer [c]) + +end + +fun matrix_simplify th = + let + val simp_th = matrix_compute (cprop_of th) + val th = Thm.strip_shyps (Thm.equal_elim simp_th th) + fun removeTrue th = removeTrue (Thm.implies_elim th TrueI) handle _ => th (* FIXME avoid handle _ *) + in + removeTrue th + end + +fun prove_bound lptfile prec = + let + val th = lp_dual_estimate_prt lptfile prec + in + matrix_simplify th + end + +val realFromStr = the o Real.fromString; +fun float2real (x, y) = realFromStr x * Math.pow (2.0, realFromStr y); + +end