--- 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
--- /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 \<le> A"
+ "A \<le> sparse_row_matrix A2"
+ "sparse_row_matrix c1 \<le> c"
+ "c \<le> sparse_row_matrix c2"
+ "sparse_row_matrix r1 \<le> x"
+ "x \<le> sparse_row_matrix r2"
+ "A * x \<le> sparse_row_matrix (b::('a::lattice_ring) spmat)"
+ shows
+ "c * x \<le> 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 \<le> A"
+ "A \<le> sparse_row_matrix A2"
+ "sparse_row_matrix c1 \<le> c"
+ "c \<le> sparse_row_matrix c2"
+ "sparse_row_matrix r1 \<le> x"
+ "x \<le> sparse_row_matrix r2"
+ "A * x \<le> sparse_row_matrix (b::('a::lattice_ring) spmat)"
+ shows
+ "c * x \<le> 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
--- /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);
--- /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;*)
--- /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;
--- 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";
--- 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 \<le> A"
- "A \<le> sparse_row_matrix A2"
- "sparse_row_matrix c1 \<le> c"
- "c \<le> sparse_row_matrix c2"
- "sparse_row_matrix r1 \<le> x"
- "x \<le> sparse_row_matrix r2"
- "A * x \<le> sparse_row_matrix (b::('a::lattice_ring) spmat)"
- shows
- "c * x \<le> 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 \<le> A"
- "A \<le> sparse_row_matrix A2"
- "sparse_row_matrix c1 \<le> c"
- "c \<le> sparse_row_matrix c2"
- "sparse_row_matrix r1 \<le> x"
- "x \<le> sparse_row_matrix r2"
- "A * x \<le> sparse_row_matrix (b::('a::lattice_ring) spmat)"
- shows
- "c * x \<le> 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
--- 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);
--- 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;*)
--- 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;
--- 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
--- 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
--- /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
--- /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