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