1.1 --- a/src/HOL/IsaMakefile Fri Sep 03 10:27:05 2004 +0200
1.2 +++ b/src/HOL/IsaMakefile Fri Sep 03 17:10:36 2004 +0200
1.3 @@ -1,5 +1,3 @@
1.4 -#
1.5 -# $Id$
1.6 #
1.7 # IsaMakefile for HOL
1.8 #
1.9 @@ -80,8 +78,7 @@
1.10 $(SRC)/Provers/quasi.ML \
1.11 $(SRC)/Provers/simplifier.ML $(SRC)/Provers/splitter.ML \
1.12 $(SRC)/Provers/trancl.ML \
1.13 - $(SRC)/TFL/isand.ML $(SRC)/TFL/casesplit.ML $(SRC)/TFL/dcterm.ML\
1.14 - $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
1.15 + $(SRC)/TFL/dcterm.ML $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \
1.16 $(SRC)/TFL/tfl.ML $(SRC)/TFL/thms.ML $(SRC)/TFL/thry.ML \
1.17 $(SRC)/TFL/usyntax.ML $(SRC)/TFL/utils.ML \
1.18 Datatype.thy Datatype_Universe.ML Datatype_Universe.thy \
1.19 @@ -97,8 +94,7 @@
1.20 Nat.thy NatArith.ML NatArith.thy \
1.21 Power.thy PreList.thy Product_Type.ML Product_Type.thy \
1.22 Refute.thy ROOT.ML \
1.23 - Recdef.thy Reconstruction.thy Record.thy\
1.24 - Relation.ML Relation.thy Relation_Power.ML \
1.25 + Recdef.thy Record.thy Relation.ML Relation.thy Relation_Power.ML \
1.26 Relation_Power.thy LOrder.thy OrderedGroup.thy OrderedGroup.ML Ring_and_Field.thy\
1.27 Set.ML Set.thy SetInterval.ML SetInterval.thy \
1.28 Sum_Type.ML Sum_Type.thy Tools/datatype_abs_proofs.ML Tools/datatype_aux.ML \
1.29 @@ -109,7 +105,7 @@
1.30 Tools/primrec_package.ML \
1.31 Tools/prop_logic.ML \
1.32 Tools/recdef_package.ML Tools/recfun_codegen.ML \
1.33 - Tools/record_package.ML Tools/reconstruction.ML\
1.34 + Tools/record_package.ML \
1.35 Tools/refute.ML Tools/refute_isar.ML \
1.36 Tools/rewrite_hol_proof.ML \
1.37 Tools/sat_solver.ML \
1.38 @@ -635,14 +631,18 @@
1.39
1.40 ## HOL-Matrix
1.41
1.42 -HOL-Matrix: HOL $(LOG)/HOL-Matrix.gz
1.43 +HOL-Matrix: HOL HOL-Complex $(LOG)/HOL-Matrix.gz
1.44
1.45 -$(LOG)/HOL-Matrix.gz: $(OUT)/HOL Matrix/ROOT.ML \
1.46 - Matrix/Matrix.thy Matrix/LinProg.thy Matrix/MatrixGeneral.thy \
1.47 - Matrix/SparseMatrix.thy Matrix/document/root.tex
1.48 - @$(ISATOOL) usedir $(OUT)/HOL Matrix
1.49 -
1.50 -
1.51 +$(LOG)/HOL-Matrix.gz: $(OUT)/HOL-Complex \
1.52 + Matrix/MatrixGeneral.thy Matrix/Matrix.thy Matrix/SparseMatrix.thy \
1.53 + Matrix/Float.thy Matrix/FloatArith.ML Matrix/ExactFloatingPoint.ML \
1.54 + Matrix/Cplex.ML Matrix/CplexMatrixConverter.ML \
1.55 + Matrix/FloatSparseMatrixBuilder.ML \
1.56 + Matrix/conv.ML Matrix/eq_codegen.ML Matrix/codegen_prep.ML \
1.57 + Matrix/fspmlp.ML \
1.58 + Matrix/MatrixLP_gensimp.ML Matrix/MatrixLP.ML Matrix/MatrixLP.thy \
1.59 + Matrix/document/root.tex Matrix/ROOT.ML
1.60 + @cd Matrix; $(ISATOOL) usedir -b $(OUT)/HOL-Complex HOL-Matrix
1.61
1.62 ## TLA
1.63
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/HOL/Matrix/Cplex.ML Fri Sep 03 17:10:36 2004 +0200
2.3 @@ -0,0 +1,1027 @@
2.4 +signature CPLEX =
2.5 +sig
2.6 +
2.7 + datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf
2.8 + | cplexNeg of cplexTerm
2.9 + | cplexProd of cplexTerm * cplexTerm
2.10 + | cplexSum of (cplexTerm list)
2.11 +
2.12 + datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
2.13 +
2.14 + datatype cplexGoal = cplexMinimize of cplexTerm
2.15 + | cplexMaximize of cplexTerm
2.16 +
2.17 + datatype cplexConstr = cplexConstr of cplexComp *
2.18 + (cplexTerm * cplexTerm)
2.19 +
2.20 + datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
2.21 + * cplexComp * cplexTerm
2.22 + | cplexBound of cplexTerm * cplexComp * cplexTerm
2.23 +
2.24 + datatype cplexProg = cplexProg of string
2.25 + * cplexGoal
2.26 + * ((string option * cplexConstr)
2.27 + list)
2.28 + * cplexBounds list
2.29 +
2.30 + datatype cplexResult = Unbounded
2.31 + | Infeasible
2.32 + | Undefined
2.33 + | Optimal of string *
2.34 + (((* name *) string *
2.35 + (* value *) string) list)
2.36 +
2.37 + exception Load_cplexFile of string
2.38 + exception Load_cplexResult of string
2.39 + exception Save_cplexFile of string
2.40 + exception Execute of string
2.41 +
2.42 + val load_cplexFile : string -> cplexProg
2.43 +
2.44 + val save_cplexFile : string -> cplexProg -> unit
2.45 +
2.46 + val elim_nonfree_bounds : cplexProg -> cplexProg
2.47 +
2.48 + val relax_strict_ineqs : cplexProg -> cplexProg
2.49 +
2.50 + val is_normed_cplexProg : cplexProg -> bool
2.51 +
2.52 + val load_cplexResult : string -> cplexResult
2.53 +
2.54 + val solve : cplexProg -> cplexResult
2.55 +end;
2.56 +
2.57 +structure Cplex : CPLEX =
2.58 +struct
2.59 +
2.60 +exception Load_cplexFile of string;
2.61 +exception Load_cplexResult of string;
2.62 +exception Save_cplexFile of string;
2.63 +
2.64 +datatype cplexTerm = cplexVar of string
2.65 + | cplexNum of string
2.66 + | cplexInf
2.67 + | cplexNeg of cplexTerm
2.68 + | cplexProd of cplexTerm * cplexTerm
2.69 + | cplexSum of (cplexTerm list)
2.70 +datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq
2.71 +datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm
2.72 +datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm)
2.73 +datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm
2.74 + * cplexComp * cplexTerm
2.75 + | cplexBound of cplexTerm * cplexComp * cplexTerm
2.76 +datatype cplexProg = cplexProg of string
2.77 + * cplexGoal
2.78 + * ((string option * cplexConstr) list)
2.79 + * cplexBounds list
2.80 +
2.81 +fun rev_cmp cplexLe = cplexGe
2.82 + | rev_cmp cplexLeq = cplexGeq
2.83 + | rev_cmp cplexGe = cplexLe
2.84 + | rev_cmp cplexGeq = cplexLeq
2.85 + | rev_cmp cplexEq = cplexEq
2.86 +
2.87 +fun the None = raise (Load_cplexFile "Some expected")
2.88 + | the (Some x) = x;
2.89 +
2.90 +fun modulo_signed is_something (cplexNeg u) = is_something u
2.91 + | modulo_signed is_something u = is_something u
2.92 +
2.93 +fun is_Num (cplexNum _) = true
2.94 + | is_Num _ = false
2.95 +
2.96 +fun is_Inf cplexInf = true
2.97 + | is_Inf _ = false
2.98 +
2.99 +fun is_Var (cplexVar _) = true
2.100 + | is_Var _ = false
2.101 +
2.102 +fun is_Neg (cplexNeg x ) = true
2.103 + | is_Neg _ = false
2.104 +
2.105 +fun is_normed_Prod (cplexProd (t1, t2)) =
2.106 + (is_Num t1) andalso (is_Var t2)
2.107 + | is_normed_Prod x = is_Var x
2.108 +
2.109 +fun is_normed_Sum (cplexSum ts) =
2.110 + (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts
2.111 + | is_normed_Sum x = modulo_signed is_normed_Prod x
2.112 +
2.113 +fun is_normed_Constr (cplexConstr (c, (t1, t2))) =
2.114 + (is_normed_Sum t1) andalso (modulo_signed is_Num t2)
2.115 +
2.116 +fun is_Num_or_Inf x = is_Inf x orelse is_Num x
2.117 +
2.118 +fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) =
2.119 + (c1 = cplexLe orelse c1 = cplexLeq) andalso
2.120 + (c2 = cplexLe orelse c2 = cplexLeq) andalso
2.121 + is_Var t2 andalso
2.122 + modulo_signed is_Num_or_Inf t1 andalso
2.123 + modulo_signed is_Num_or_Inf t3
2.124 + | is_normed_Bounds (cplexBound (t1, c, t2)) =
2.125 + (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2))
2.126 + orelse
2.127 + (c <> cplexEq andalso
2.128 + is_Var t2 andalso (modulo_signed is_Num_or_Inf t1))
2.129 +
2.130 +fun term_of_goal (cplexMinimize x) = x
2.131 + | term_of_goal (cplexMaximize x) = x
2.132 +
2.133 +fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) =
2.134 + is_normed_Sum (term_of_goal goal) andalso
2.135 + forall (fn (_,x) => is_normed_Constr x) constraints andalso
2.136 + forall is_normed_Bounds bounds
2.137 +
2.138 +fun is_NL s = s = "\n"
2.139 +
2.140 +fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s)
2.141 +
2.142 +fun is_num a =
2.143 + let
2.144 + val b = String.explode a
2.145 + fun num4 cs = forall Char.isDigit cs
2.146 + fun num3 [] = true
2.147 + | num3 (ds as (c::cs)) =
2.148 + if c = #"+" orelse c = #"-" then
2.149 + num4 cs
2.150 + else
2.151 + num4 ds
2.152 + fun num2 [] = true
2.153 + | num2 (c::cs) =
2.154 + if c = #"e" orelse c = #"E" then num3 cs
2.155 + else (Char.isDigit c) andalso num2 cs
2.156 + fun num1 [] = true
2.157 + | num1 (c::cs) =
2.158 + if c = #"." then num2 cs
2.159 + else if c = #"e" orelse c = #"E" then num3 cs
2.160 + else (Char.isDigit c) andalso num1 cs
2.161 + fun num [] = true
2.162 + | num (c::cs) =
2.163 + if c = #"." then num2 cs
2.164 + else (Char.isDigit c) andalso num1 cs
2.165 + in
2.166 + num b
2.167 + end
2.168 +
2.169 +fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":"
2.170 +
2.171 +fun is_cmp s = s = "<" orelse s = ">" orelse s = "<="
2.172 + orelse s = ">=" orelse s = "="
2.173 +
2.174 +fun is_symbol a =
2.175 + let
2.176 + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~"
2.177 + fun is_symbol_char c = Char.isAlphaNum c orelse
2.178 + exists (fn d => d=c) symbol_char
2.179 + fun is_symbol_start c = is_symbol_char c andalso
2.180 + not (Char.isDigit c) andalso
2.181 + not (c= #".")
2.182 + val b = String.explode a
2.183 + in
2.184 + b <> [] andalso is_symbol_start (hd b) andalso
2.185 + forall is_symbol_char b
2.186 + end
2.187 +
2.188 +fun to_upper s = String.implode (map Char.toUpper (String.explode s))
2.189 +
2.190 +fun keyword x =
2.191 + let
2.192 + val a = to_upper x
2.193 + in
2.194 + if a = "BOUNDS" orelse a = "BOUND" then
2.195 + Some "BOUNDS"
2.196 + else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then
2.197 + Some "MINIMIZE"
2.198 + else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then
2.199 + Some "MAXIMIZE"
2.200 + else if a = "ST" orelse a = "S.T." orelse a = "ST." then
2.201 + Some "ST"
2.202 + else if a = "FREE" orelse a = "END" then
2.203 + Some a
2.204 + else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then
2.205 + Some "GENERAL"
2.206 + else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then
2.207 + Some "INTEGER"
2.208 + else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then
2.209 + Some "BINARY"
2.210 + else if a = "INF" orelse a = "INFINITY" then
2.211 + Some "INF"
2.212 + else
2.213 + None
2.214 + end
2.215 +
2.216 +val TOKEN_ERROR = ~1
2.217 +val TOKEN_BLANK = 0
2.218 +val TOKEN_NUM = 1
2.219 +val TOKEN_DELIMITER = 2
2.220 +val TOKEN_SYMBOL = 3
2.221 +val TOKEN_LABEL = 4
2.222 +val TOKEN_CMP = 5
2.223 +val TOKEN_KEYWORD = 6
2.224 +val TOKEN_NL = 7
2.225 +
2.226 +(* tokenize takes a list of chars as argument and returns a list of
2.227 + int * string pairs, each string representing a "cplex token",
2.228 + and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP
2.229 + or TOKEN_SYMBOL *)
2.230 +fun tokenize s =
2.231 + let
2.232 + val flist = [(is_NL, TOKEN_NL),
2.233 + (is_blank, TOKEN_BLANK),
2.234 + (is_num, TOKEN_NUM),
2.235 + (is_delimiter, TOKEN_DELIMITER),
2.236 + (is_cmp, TOKEN_CMP),
2.237 + (is_symbol, TOKEN_SYMBOL)]
2.238 + fun match_helper [] s = (fn x => false, TOKEN_ERROR)
2.239 + | match_helper (f::fs) s =
2.240 + if ((fst f) s) then f else match_helper fs s
2.241 + fun match s = match_helper flist s
2.242 + fun tok s =
2.243 + if s = "" then [] else
2.244 + let
2.245 + val h = String.substring (s,0,1)
2.246 + val (f, j) = match h
2.247 + fun len i =
2.248 + if size s = i then i
2.249 + else if f (String.substring (s,0,i+1)) then
2.250 + len (i+1)
2.251 + else i
2.252 + in
2.253 + if j < 0 then
2.254 + (if h = "\\" then []
2.255 + else raise (Load_cplexFile ("token expected, found: "
2.256 + ^s)))
2.257 + else
2.258 + let
2.259 + val l = len 1
2.260 + val u = String.substring (s,0,l)
2.261 + val v = String.extract (s,l,NONE)
2.262 + in
2.263 + if j = 0 then tok v else (j, u) :: tok v
2.264 + end
2.265 + end
2.266 + in
2.267 + tok s
2.268 + end
2.269 +
2.270 +exception Tokenize of string;
2.271 +
2.272 +fun tokenize_general flist s =
2.273 + let
2.274 + fun match_helper [] s = raise (Tokenize s)
2.275 + | match_helper (f::fs) s =
2.276 + if ((fst f) s) then f else match_helper fs s
2.277 + fun match s = match_helper flist s
2.278 + fun tok s =
2.279 + if s = "" then [] else
2.280 + let
2.281 + val h = String.substring (s,0,1)
2.282 + val (f, j) = match h
2.283 + fun len i =
2.284 + if size s = i then i
2.285 + else if f (String.substring (s,0,i+1)) then
2.286 + len (i+1)
2.287 + else i
2.288 + val l = len 1
2.289 + in
2.290 + (j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE))
2.291 + end
2.292 + in
2.293 + tok s
2.294 + end
2.295 +
2.296 +fun load_cplexFile name =
2.297 + let
2.298 + val f = TextIO.openIn name
2.299 + val ignore_NL = ref true
2.300 + val rest = ref []
2.301 +
2.302 + fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s
2.303 +
2.304 + fun readToken_helper () =
2.305 + if length (!rest) > 0 then
2.306 + let val u = hd (!rest) in
2.307 + (
2.308 + rest := tl (!rest);
2.309 + Some u
2.310 + )
2.311 + end
2.312 + else
2.313 + let val s = TextIO.inputLine f in
2.314 + if s = "" then None else
2.315 + let val t = tokenize s in
2.316 + if (length t >= 2 andalso
2.317 + snd(hd (tl t)) = ":")
2.318 + then
2.319 + rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t))
2.320 + else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t))
2.321 + andalso is_symbol "TO" (hd (tl t))
2.322 + then
2.323 + rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t))
2.324 + else
2.325 + rest := t;
2.326 + readToken_helper ()
2.327 + end
2.328 + end
2.329 +
2.330 + fun readToken_helper2 () =
2.331 + let val c = readToken_helper () in
2.332 + if c = None then None
2.333 + else if !ignore_NL andalso fst (the c) = TOKEN_NL then
2.334 + readToken_helper2 ()
2.335 + else if fst (the c) = TOKEN_SYMBOL
2.336 + andalso keyword (snd (the c)) <> None
2.337 + then Some (TOKEN_KEYWORD, the (keyword (snd (the c))))
2.338 + else c
2.339 + end
2.340 +
2.341 + fun readToken () = readToken_helper2 ()
2.342 +
2.343 + fun pushToken a = rest := (a::(!rest))
2.344 +
2.345 + fun is_value token =
2.346 + fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD
2.347 + andalso snd token = "INF")
2.348 +
2.349 + fun get_value token =
2.350 + if fst token = TOKEN_NUM then
2.351 + cplexNum (snd token)
2.352 + else if fst token = TOKEN_KEYWORD andalso snd token = "INF"
2.353 + then
2.354 + cplexInf
2.355 + else
2.356 + raise (Load_cplexFile "num expected")
2.357 +
2.358 + fun readTerm_Product only_num =
2.359 + let val c = readToken () in
2.360 + if c = None then None
2.361 + else if fst (the c) = TOKEN_SYMBOL
2.362 + then (
2.363 + if only_num then (pushToken (the c); None)
2.364 + else Some (cplexVar (snd (the c)))
2.365 + )
2.366 + else if only_num andalso is_value (the c) then
2.367 + Some (get_value (the c))
2.368 + else if is_value (the c) then
2.369 + let val t1 = get_value (the c)
2.370 + val d = readToken ()
2.371 + in
2.372 + if d = None then Some t1
2.373 + else if fst (the d) = TOKEN_SYMBOL then
2.374 + Some (cplexProd (t1, cplexVar (snd (the d))))
2.375 + else
2.376 + (pushToken (the d); Some t1)
2.377 + end
2.378 + else (pushToken (the c); None)
2.379 + end
2.380 +
2.381 + fun readTerm_Signed only_signed only_num =
2.382 + let
2.383 + val c = readToken ()
2.384 + in
2.385 + if c = None then None
2.386 + else
2.387 + let val d = the c in
2.388 + if d = (TOKEN_DELIMITER, "+") then
2.389 + readTerm_Product only_num
2.390 + else if d = (TOKEN_DELIMITER, "-") then
2.391 + Some (cplexNeg (the (readTerm_Product
2.392 + only_num)))
2.393 + else (pushToken d;
2.394 + if only_signed then None
2.395 + else readTerm_Product only_num)
2.396 + end
2.397 + end
2.398 +
2.399 + fun readTerm_Sum first_signed =
2.400 + let val c = readTerm_Signed first_signed false in
2.401 + if c = None then [] else (the c)::(readTerm_Sum true)
2.402 + end
2.403 +
2.404 + fun readTerm () =
2.405 + let val c = readTerm_Sum false in
2.406 + if c = [] then None
2.407 + else if tl c = [] then Some (hd c)
2.408 + else Some (cplexSum c)
2.409 + end
2.410 +
2.411 + fun readLabeledTerm () =
2.412 + let val c = readToken () in
2.413 + if c = None then (None, None)
2.414 + else if fst (the c) = TOKEN_LABEL then
2.415 + let val t = readTerm () in
2.416 + if t = None then
2.417 + raise (Load_cplexFile ("term after label "^
2.418 + (snd (the c))^
2.419 + " expected"))
2.420 + else (Some (snd (the c)), t)
2.421 + end
2.422 + else (pushToken (the c); (None, readTerm ()))
2.423 + end
2.424 +
2.425 + fun readGoal () =
2.426 + let
2.427 + val g = readToken ()
2.428 + in
2.429 + if g = Some (TOKEN_KEYWORD, "MAXIMIZE") then
2.430 + cplexMaximize (the (snd (readLabeledTerm ())))
2.431 + else if g = Some (TOKEN_KEYWORD, "MINIMIZE") then
2.432 + cplexMinimize (the (snd (readLabeledTerm ())))
2.433 + else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected")
2.434 + end
2.435 +
2.436 + fun str2cmp b =
2.437 + (case b of
2.438 + "<" => cplexLe
2.439 + | "<=" => cplexLeq
2.440 + | ">" => cplexGe
2.441 + | ">=" => cplexGeq
2.442 + | "=" => cplexEq
2.443 + | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP")))
2.444 +
2.445 + fun readConstraint () =
2.446 + let
2.447 + val t = readLabeledTerm ()
2.448 + fun make_constraint b t1 t2 =
2.449 + cplexConstr
2.450 + (str2cmp b,
2.451 + (t1, t2))
2.452 + in
2.453 + if snd t = None then None
2.454 + else
2.455 + let val c = readToken () in
2.456 + if c = None orelse fst (the c) <> TOKEN_CMP
2.457 + then raise (Load_cplexFile "TOKEN_CMP expected")
2.458 + else
2.459 + let val n = readTerm_Signed false true in
2.460 + if n = None then
2.461 + raise (Load_cplexFile "num expected")
2.462 + else
2.463 + Some (fst t,
2.464 + make_constraint (snd (the c))
2.465 + (the (snd t))
2.466 + (the n))
2.467 + end
2.468 + end
2.469 + end
2.470 +
2.471 + fun readST () =
2.472 + let
2.473 + fun readbody () =
2.474 + let val t = readConstraint () in
2.475 + if t = None then []
2.476 + else if (is_normed_Constr (snd (the t))) then
2.477 + (the t)::(readbody ())
2.478 + else if (fst (the t) <> None) then
2.479 + raise (Load_cplexFile
2.480 + ("constraint '"^(the (fst (the t)))^
2.481 + "'is not normed"))
2.482 + else
2.483 + raise (Load_cplexFile
2.484 + "constraint is not normed")
2.485 + end
2.486 + in
2.487 + if readToken () = Some (TOKEN_KEYWORD, "ST")
2.488 + then
2.489 + readbody ()
2.490 + else
2.491 + raise (Load_cplexFile "ST expected")
2.492 + end
2.493 +
2.494 + fun readCmp () =
2.495 + let val c = readToken () in
2.496 + if c = None then None
2.497 + else if fst (the c) = TOKEN_CMP then
2.498 + Some (str2cmp (snd (the c)))
2.499 + else (pushToken (the c); None)
2.500 + end
2.501 +
2.502 + fun skip_NL () =
2.503 + let val c = readToken () in
2.504 + if c <> None andalso fst (the c) = TOKEN_NL then
2.505 + skip_NL ()
2.506 + else
2.507 + (pushToken (the c); ())
2.508 + end
2.509 +
2.510 + fun is_var (cplexVar _) = true
2.511 + | is_var _ = false
2.512 +
2.513 + fun make_bounds c t1 t2 =
2.514 + cplexBound (t1, c, t2)
2.515 +
2.516 + fun readBound () =
2.517 + let
2.518 + val _ = skip_NL ()
2.519 + val t1 = readTerm ()
2.520 + in
2.521 + if t1 = None then None
2.522 + else
2.523 + let
2.524 + val c1 = readCmp ()
2.525 + in
2.526 + if c1 = None then
2.527 + let
2.528 + val c = readToken ()
2.529 + in
2.530 + if c = Some (TOKEN_KEYWORD, "FREE") then
2.531 + Some (
2.532 + cplexBounds (cplexNeg cplexInf,
2.533 + cplexLeq,
2.534 + the t1,
2.535 + cplexLeq,
2.536 + cplexInf))
2.537 + else
2.538 + raise (Load_cplexFile "FREE expected")
2.539 + end
2.540 + else
2.541 + let
2.542 + val t2 = readTerm ()
2.543 + in
2.544 + if t2 = None then
2.545 + raise (Load_cplexFile "term expected")
2.546 + else
2.547 + let val c2 = readCmp () in
2.548 + if c2 = None then
2.549 + Some (make_bounds (the c1)
2.550 + (the t1)
2.551 + (the t2))
2.552 + else
2.553 + Some (
2.554 + cplexBounds (the t1,
2.555 + the c1,
2.556 + the t2,
2.557 + the c2,
2.558 + the (readTerm())))
2.559 + end
2.560 + end
2.561 + end
2.562 + end
2.563 +
2.564 + fun readBounds () =
2.565 + let
2.566 + fun makestring b = "?"
2.567 + fun readbody () =
2.568 + let
2.569 + val b = readBound ()
2.570 + in
2.571 + if b = None then []
2.572 + else if (is_normed_Bounds (the b)) then
2.573 + (the b)::(readbody())
2.574 + else (
2.575 + raise (Load_cplexFile
2.576 + ("bounds are not normed in: "^
2.577 + (makestring (the b)))))
2.578 + end
2.579 + in
2.580 + if readToken () = Some (TOKEN_KEYWORD, "BOUNDS") then
2.581 + readbody ()
2.582 + else raise (Load_cplexFile "BOUNDS expected")
2.583 + end
2.584 +
2.585 + fun readEnd () =
2.586 + if readToken () = Some (TOKEN_KEYWORD, "END") then ()
2.587 + else raise (Load_cplexFile "END expected")
2.588 +
2.589 + val result_Goal = readGoal ()
2.590 + val result_ST = readST ()
2.591 + val _ = ignore_NL := false
2.592 + val result_Bounds = readBounds ()
2.593 + val _ = ignore_NL := true
2.594 + val _ = readEnd ()
2.595 + val _ = TextIO.closeIn f
2.596 + in
2.597 + cplexProg (name, result_Goal, result_ST, result_Bounds)
2.598 + end
2.599 +
2.600 +fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) =
2.601 + let
2.602 + val f = TextIO.openOut filename
2.603 +
2.604 + fun basic_write s = TextIO.output(f, s)
2.605 +
2.606 + val linebuf = ref ""
2.607 + fun buf_flushline s =
2.608 + (basic_write (!linebuf);
2.609 + basic_write "\n";
2.610 + linebuf := s)
2.611 + fun buf_add s = linebuf := (!linebuf) ^ s
2.612 +
2.613 + fun write s =
2.614 + if (String.size s) + (String.size (!linebuf)) >= 250 then
2.615 + buf_flushline (" "^s)
2.616 + else
2.617 + buf_add s
2.618 +
2.619 + fun writeln s = (buf_add s; buf_flushline "")
2.620 +
2.621 + fun write_term (cplexVar x) = write x
2.622 + | write_term (cplexNum x) = write x
2.623 + | write_term cplexInf = write "inf"
2.624 + | write_term (cplexProd (cplexNum "1", b)) = write_term b
2.625 + | write_term (cplexProd (a, b)) =
2.626 + (write_term a; write " "; write_term b)
2.627 + | write_term (cplexNeg x) = (write " - "; write_term x)
2.628 + | write_term (cplexSum ts) = write_terms ts
2.629 + and write_terms [] = ()
2.630 + | write_terms (t::ts) =
2.631 + (if (not (is_Neg t)) then write " + " else ();
2.632 + write_term t; write_terms ts)
2.633 +
2.634 + fun write_goal (cplexMaximize term) =
2.635 + (writeln "MAXIMIZE"; write_term term; writeln "")
2.636 + | write_goal (cplexMinimize term) =
2.637 + (writeln "MINIMIZE"; write_term term; writeln "")
2.638 +
2.639 + fun write_cmp cplexLe = write "<"
2.640 + | write_cmp cplexLeq = write "<="
2.641 + | write_cmp cplexEq = write "="
2.642 + | write_cmp cplexGe = write ">"
2.643 + | write_cmp cplexGeq = write ">="
2.644 +
2.645 + fun write_constr (cplexConstr (cmp, (a,b))) =
2.646 + (write_term a;
2.647 + write " ";
2.648 + write_cmp cmp;
2.649 + write " ";
2.650 + write_term b)
2.651 +
2.652 + fun write_constraints [] = ()
2.653 + | write_constraints (c::cs) =
2.654 + (if (fst c <> None)
2.655 + then
2.656 + (write (the (fst c)); write ": ")
2.657 + else
2.658 + ();
2.659 + write_constr (snd c);
2.660 + writeln "";
2.661 + write_constraints cs)
2.662 +
2.663 + fun write_bounds [] = ()
2.664 + | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) =
2.665 + ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf
2.666 + andalso (c1 = cplexLeq orelse c1 = cplexLe)
2.667 + andalso (c2 = cplexLeq orelse c2 = cplexLe)
2.668 + then
2.669 + (write_term t2; write " free")
2.670 + else
2.671 + (write_term t1; write " "; write_cmp c1; write " ";
2.672 + write_term t2; write " "; write_cmp c2; write " ";
2.673 + write_term t3)
2.674 + ); writeln ""; write_bounds bs)
2.675 + | write_bounds ((cplexBound (t1, c, t2)) :: bs) =
2.676 + (write_term t1; write " ";
2.677 + write_cmp c; write " ";
2.678 + write_term t2; writeln ""; write_bounds bs)
2.679 +
2.680 + val _ = write_goal goal
2.681 + val _ = (writeln ""; writeln "ST")
2.682 + val _ = write_constraints constraints
2.683 + val _ = (writeln ""; writeln "BOUNDS")
2.684 + val _ = write_bounds bounds
2.685 + val _ = (writeln ""; writeln "END")
2.686 + val _ = TextIO.closeOut f
2.687 + in
2.688 + ()
2.689 + end
2.690 +
2.691 +fun norm_Constr (constr as cplexConstr (c, (t1, t2))) =
2.692 + if not (modulo_signed is_Num t2) andalso
2.693 + modulo_signed is_Num t1
2.694 + then
2.695 + [cplexConstr (rev_cmp c, (t2, t1))]
2.696 + else if (c = cplexLe orelse c = cplexLeq) andalso
2.697 + (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf)
2.698 + then
2.699 + []
2.700 + else if (c = cplexGe orelse c = cplexGeq) andalso
2.701 + (t1 = cplexInf orelse t2 = cplexNeg cplexInf)
2.702 + then
2.703 + []
2.704 + else
2.705 + [constr]
2.706 +
2.707 +fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) =
2.708 + (norm_Constr(cplexConstr (c1, (t1, t2))))
2.709 + @ (norm_Constr(cplexConstr (c2, (t2, t3))))
2.710 + | bound2constr (cplexBound (t1, cplexEq, t2)) =
2.711 + (norm_Constr(cplexConstr (cplexLeq, (t1, t2))))
2.712 + @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1))))
2.713 + | bound2constr (cplexBound (t1, c1, t2)) =
2.714 + norm_Constr(cplexConstr (c1, (t1,t2)))
2.715 +
2.716 +val emptyset = Symtab.empty
2.717 +
2.718 +fun singleton v = Symtab.update ((v, ()), emptyset)
2.719 +
2.720 +fun merge a b = Symtab.merge (op =) (a, b)
2.721 +
2.722 +fun mergemap f ts = foldl
2.723 + (fn (table, x) => merge table (f x))
2.724 + (Symtab.empty, ts)
2.725 +
2.726 +fun diff a b = Symtab.foldl (fn (a, (k,v)) =>
2.727 + (Symtab.delete k a) handle UNDEF => a)
2.728 + (a, b)
2.729 +
2.730 +fun collect_vars (cplexVar v) = singleton v
2.731 + | collect_vars (cplexNeg t) = collect_vars t
2.732 + | collect_vars (cplexProd (t1, t2)) =
2.733 + merge (collect_vars t1) (collect_vars t2)
2.734 + | collect_vars (cplexSum ts) = mergemap collect_vars ts
2.735 + | collect_vars _ = emptyset
2.736 +
2.737 +(* Eliminates all nonfree bounds from the linear program and produces an
2.738 + equivalent program with only free bounds
2.739 + IF for the input program P holds: is_normed_cplexProg P *)
2.740 +fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) =
2.741 + let
2.742 + fun collect_constr_vars (_, cplexConstr (c, (t1,_))) =
2.743 + (collect_vars t1)
2.744 +
2.745 + val cvars = merge (collect_vars (term_of_goal goal))
2.746 + (mergemap collect_constr_vars constraints)
2.747 +
2.748 + fun collect_lower_bounded_vars
2.749 + (cplexBounds (t1, c1, cplexVar v, c2, t3)) =
2.750 + singleton v
2.751 + | collect_lower_bounded_vars
2.752 + (cplexBound (_, cplexLe, cplexVar v)) =
2.753 + singleton v
2.754 + | collect_lower_bounded_vars
2.755 + (cplexBound (_, cplexLeq, cplexVar v)) =
2.756 + singleton v
2.757 + | collect_lower_bounded_vars
2.758 + (cplexBound (cplexVar v, cplexGe,_)) =
2.759 + singleton v
2.760 + | collect_lower_bounded_vars
2.761 + (cplexBound (cplexVar v, cplexGeq, _)) =
2.762 + singleton v
2.763 + | collect_lower_bounded_vars
2.764 + (cplexBound (cplexVar v, cplexEq, _)) =
2.765 + singleton v
2.766 + | collect_lower_bounded_vars _ = emptyset
2.767 +
2.768 + val lvars = mergemap collect_lower_bounded_vars bounds
2.769 + val positive_vars = diff cvars lvars
2.770 + val zero = cplexNum "0"
2.771 +
2.772 + fun make_pos_constr v =
2.773 + (None, cplexConstr (cplexGeq, ((cplexVar v), zero)))
2.774 +
2.775 + fun make_free_bound v =
2.776 + cplexBounds (cplexNeg cplexInf, cplexLeq,
2.777 + cplexVar v, cplexLeq,
2.778 + cplexInf)
2.779 +
2.780 + val pos_constrs = rev(Symtab.foldl
2.781 + (fn (l, (k,v)) => (make_pos_constr k) :: l)
2.782 + ([], positive_vars))
2.783 + val bound_constrs = map (fn x => (None, x))
2.784 + (flat (map bound2constr bounds))
2.785 + val constraints' = constraints @ pos_constrs @ bound_constrs
2.786 + val bounds' = rev(Symtab.foldl (fn (l, (v,_)) =>
2.787 + (make_free_bound v)::l)
2.788 + ([], cvars))
2.789 + in
2.790 + cplexProg (name, goal, constraints', bounds')
2.791 + end
2.792 +
2.793 +fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) =
2.794 + let
2.795 + fun relax cplexLe = cplexLeq
2.796 + | relax cplexGe = cplexGeq
2.797 + | relax x = x
2.798 +
2.799 + fun relax_constr (n, cplexConstr(c, (t1, t2))) =
2.800 + (n, cplexConstr(relax c, (t1, t2)))
2.801 +
2.802 + fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) =
2.803 + cplexBounds (t1, relax c1, t2, relax c2, t3)
2.804 + | relax_bounds (cplexBound (t1, c, t2)) =
2.805 + cplexBound (t1, relax c, t2)
2.806 + in
2.807 + cplexProg (name,
2.808 + goals,
2.809 + map relax_constr constrs,
2.810 + map relax_bounds bounds)
2.811 + end
2.812 +
2.813 +datatype cplexResult = Unbounded
2.814 + | Infeasible
2.815 + | Undefined
2.816 + | Optimal of string * ((string * string) list)
2.817 +
2.818 +fun is_separator x = forall (fn c => c = #"-") (String.explode x)
2.819 +
2.820 +fun is_sign x = (x = "+" orelse x = "-")
2.821 +
2.822 +fun is_colon x = (x = ":")
2.823 +
2.824 +fun is_resultsymbol a =
2.825 + let
2.826 + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-"
2.827 + fun is_symbol_char c = Char.isAlphaNum c orelse
2.828 + exists (fn d => d=c) symbol_char
2.829 + fun is_symbol_start c = is_symbol_char c andalso
2.830 + not (Char.isDigit c) andalso
2.831 + not (c= #".") andalso
2.832 + not (c= #"-")
2.833 + val b = String.explode a
2.834 + in
2.835 + b <> [] andalso is_symbol_start (hd b) andalso
2.836 + forall is_symbol_char b
2.837 + end
2.838 +
2.839 +val TOKEN_SIGN = 100
2.840 +val TOKEN_COLON = 101
2.841 +val TOKEN_SEPARATOR = 102
2.842 +
2.843 +fun load_cplexResult name =
2.844 + let
2.845 + val flist = [(is_NL, TOKEN_NL),
2.846 + (is_blank, TOKEN_BLANK),
2.847 + (is_num, TOKEN_NUM),
2.848 + (is_sign, TOKEN_SIGN),
2.849 + (is_colon, TOKEN_COLON),
2.850 + (is_cmp, TOKEN_CMP),
2.851 + (is_resultsymbol, TOKEN_SYMBOL),
2.852 + (is_separator, TOKEN_SEPARATOR)]
2.853 +
2.854 + val tokenize = tokenize_general flist
2.855 +
2.856 + val f = TextIO.openIn name
2.857 +
2.858 + val rest = ref []
2.859 +
2.860 + fun readToken_helper () =
2.861 + if length (!rest) > 0 then
2.862 + let val u = hd (!rest) in
2.863 + (
2.864 + rest := tl (!rest);
2.865 + Some u
2.866 + )
2.867 + end
2.868 + else
2.869 + let val s = TextIO.inputLine f in
2.870 + if s = "" then None else (rest := tokenize s; readToken_helper())
2.871 + end
2.872 +
2.873 + fun is_tt tok ty = (tok <> None andalso (fst (the tok)) = ty)
2.874 +
2.875 + fun pushToken a = if a = None then () else (rest := ((the a)::(!rest)))
2.876 +
2.877 + fun readToken () =
2.878 + let val t = readToken_helper () in
2.879 + if is_tt t TOKEN_BLANK then
2.880 + readToken ()
2.881 + else if is_tt t TOKEN_NL then
2.882 + let val t2 = readToken_helper () in
2.883 + if is_tt t2 TOKEN_SIGN then
2.884 + (pushToken (Some (TOKEN_SEPARATOR, "-")); t)
2.885 + else
2.886 + (pushToken t2; t)
2.887 + end
2.888 + else if is_tt t TOKEN_SIGN then
2.889 + let val t2 = readToken_helper () in
2.890 + if is_tt t2 TOKEN_NUM then
2.891 + (Some (TOKEN_NUM, (snd (the t))^(snd (the t2))))
2.892 + else
2.893 + (pushToken t2; t)
2.894 + end
2.895 + else
2.896 + t
2.897 + end
2.898 +
2.899 + fun readRestOfLine P =
2.900 + let
2.901 + val t = readToken ()
2.902 + in
2.903 + if is_tt t TOKEN_NL orelse t = None
2.904 + then P
2.905 + else readRestOfLine P
2.906 + end
2.907 +
2.908 + fun readHeader () =
2.909 + let
2.910 + fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ())))
2.911 + fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ())))
2.912 + val t1 = readToken ()
2.913 + val t2 = readToken ()
2.914 + in
2.915 + if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON
2.916 + then
2.917 + case to_upper (snd (the t1)) of
2.918 + "STATUS" => (readStatus ())::(readHeader ())
2.919 + | "OBJECTIVE" => (readObjective())::(readHeader ())
2.920 + | _ => (readRestOfLine (); readHeader ())
2.921 + else
2.922 + (pushToken t2; pushToken t1; [])
2.923 + end
2.924 +
2.925 + fun skip_until_sep () =
2.926 + let val x = readToken () in
2.927 + if is_tt x TOKEN_SEPARATOR then
2.928 + readRestOfLine ()
2.929 + else
2.930 + skip_until_sep ()
2.931 + end
2.932 +
2.933 + fun load_value () =
2.934 + let
2.935 + val t1 = readToken ()
2.936 + val t2 = readToken ()
2.937 + in
2.938 + if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then
2.939 + let
2.940 + val t = readToken ()
2.941 + val state = if is_tt t TOKEN_NL then readToken () else t
2.942 + val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected")
2.943 + val k = readToken ()
2.944 + in
2.945 + if is_tt k TOKEN_NUM then
2.946 + readRestOfLine (Some (snd (the t2), snd (the k)))
2.947 + else
2.948 + raise (Load_cplexResult "number expected")
2.949 + end
2.950 + else
2.951 + (pushToken t2; pushToken t1; None)
2.952 + end
2.953 +
2.954 + fun load_values () =
2.955 + let val v = load_value () in
2.956 + if v = None then [] else (the v)::(load_values ())
2.957 + end
2.958 +
2.959 + val header = readHeader ()
2.960 +
2.961 + val result =
2.962 + case assoc (header, "STATUS") of
2.963 + Some "INFEASIBLE" => Infeasible
2.964 + | Some "UNBOUNDED" => Unbounded
2.965 + | Some "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")),
2.966 + (skip_until_sep ();
2.967 + skip_until_sep ();
2.968 + load_values ()))
2.969 + | _ => Undefined
2.970 +
2.971 + val _ = TextIO.closeIn f
2.972 + in
2.973 + result
2.974 + end
2.975 + handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s))
2.976 + | OPTION => raise (Load_cplexResult "OPTION")
2.977 + | x => raise x
2.978 +
2.979 +exception Execute of string;
2.980 +
2.981 +fun execute_cplex lpname resultname =
2.982 +let
2.983 + fun wrap s = "\""^s^"\""
2.984 + val cplex_path = getenv "CPLEX"
2.985 + val cplex = if cplex_path = "" then "glpsol" else cplex_path
2.986 + val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname)
2.987 +in
2.988 + execute command
2.989 +end
2.990 +
2.991 +fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s])))
2.992 +
2.993 +fun solve prog =
2.994 + let
2.995 + val name = LargeInt.toString (Time.toMicroseconds (Time.now ()))
2.996 + val lpname = tmp_file (name^".lp")
2.997 + val resultname = tmp_file (name^".result")
2.998 + val _ = save_cplexFile lpname prog
2.999 + val answer = execute_cplex lpname resultname
2.1000 + in
2.1001 + let
2.1002 + val result = load_cplexResult resultname
2.1003 + val _ = OS.FileSys.remove lpname
2.1004 + val _ = OS.FileSys.remove resultname
2.1005 + in
2.1006 + result
2.1007 + end
2.1008 + handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer))
2.1009 + | _ => raise (Execute answer)
2.1010 + end
2.1011 +end;
2.1012 +
2.1013 +val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45"
2.1014 +val demoout = "/home/obua/flyspeck/kepler/LP/test.out"
2.1015 +val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol"
2.1016 +
2.1017 +fun loadcplex () = Cplex.relax_strict_ineqs
2.1018 + (Cplex.load_cplexFile demofile)
2.1019 +
2.1020 +fun writecplex lp = Cplex.save_cplexFile demoout lp
2.1021 +
2.1022 +fun test () =
2.1023 + let
2.1024 + val lp = loadcplex ()
2.1025 + val lp2 = Cplex.elim_nonfree_bounds lp
2.1026 + in
2.1027 + writecplex lp2
2.1028 + end
2.1029 +
2.1030 +fun loadresult () = Cplex.load_cplexResult demoresult;
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/src/HOL/Matrix/CplexMatrixConverter.ML Fri Sep 03 17:10:36 2004 +0200
3.3 @@ -0,0 +1,132 @@
3.4 +signature MATRIX_BUILDER =
3.5 +sig
3.6 + type vector
3.7 + type matrix
3.8 +
3.9 + val empty_vector : vector
3.10 + val empty_matrix : matrix
3.11 +
3.12 + exception Nat_expected of int
3.13 + val set_elem : vector -> int -> string -> vector
3.14 + val set_vector : matrix -> int -> vector -> matrix
3.15 +end;
3.16 +
3.17 +signature CPLEX_MATRIX_CONVERTER =
3.18 +sig
3.19 + structure cplex : CPLEX
3.20 + structure matrix_builder : MATRIX_BUILDER
3.21 + type vector = matrix_builder.vector
3.22 + type matrix = matrix_builder.matrix
3.23 + type naming = int * (int -> string) * (string -> int)
3.24 +
3.25 + exception Converter of string
3.26 +
3.27 + (* program must fulfill is_normed_cplexProg and must be an element of the image of elim_nonfree_bounds *)
3.28 + (* convert_prog maximize c A b naming *)
3.29 + val convert_prog : cplex.cplexProg -> bool * vector * matrix * vector * naming
3.30 +
3.31 + (* results must be optimal, converts_results returns the optimal value as string and the solution as vector *)
3.32 + (* convert_results results name2index *)
3.33 + val convert_results : cplex.cplexResult -> (string -> int) -> string * vector
3.34 +end;
3.35 +
3.36 +functor MAKE_CPLEX_MATRIX_CONVERTER (structure cplex: CPLEX and matrix_builder: MATRIX_BUILDER) : CPLEX_MATRIX_CONVERTER =
3.37 +struct
3.38 +
3.39 +structure cplex = cplex
3.40 +structure matrix_builder = matrix_builder
3.41 +type matrix = matrix_builder.matrix
3.42 +type vector = matrix_builder.vector
3.43 +type naming = int * (int -> string) * (string -> int)
3.44 +
3.45 +open matrix_builder
3.46 +open cplex
3.47 +
3.48 +exception Converter of string;
3.49 +
3.50 +structure Inttab = TableFun(type key = int val ord = int_ord);
3.51 +
3.52 +fun neg_term (cplexNeg t) = t
3.53 + | neg_term (cplexSum ts) = cplexSum (map neg_term ts)
3.54 + | neg_term t = cplexNeg t
3.55 +
3.56 +fun convert_prog (cplexProg (s, goal, constrs, bounds)) =
3.57 + let
3.58 + fun build_naming index i2s s2i [] = (index, i2s, s2i)
3.59 + | build_naming index i2s s2i (cplexBounds (cplexNeg cplexInf, cplexLeq, cplexVar v, cplexLeq, cplexInf)::bounds)
3.60 + = build_naming (index+1) (Inttab.update ((index, v), i2s)) (Symtab.update_new ((v, index), s2i)) bounds
3.61 + | build_naming _ _ _ _ = raise (Converter "nonfree bound")
3.62 +
3.63 + val (varcount, i2s_tab, s2i_tab) = build_naming 0 Inttab.empty Symtab.empty bounds
3.64 +
3.65 + fun i2s i = case Inttab.lookup (i2s_tab, i) of None => raise (Converter "index not found")
3.66 + | Some n => n
3.67 + fun s2i s = case Symtab.lookup (s2i_tab, s) of None => raise (Converter ("name not found: "^s))
3.68 + | Some i => i
3.69 + fun num2str positive (cplexNeg t) = num2str (not positive) t
3.70 + | num2str positive (cplexNum num) = if positive then num else "-"^num
3.71 + | num2str _ _ = raise (Converter "term is not a (possibly signed) number")
3.72 +
3.73 + fun setprod vec positive (cplexNeg t) = setprod vec (not positive) t
3.74 + | setprod vec positive (cplexVar v) = set_elem vec (s2i v) (if positive then "1" else "-1")
3.75 + | setprod vec positive (cplexProd (cplexNum num, cplexVar v)) =
3.76 + set_elem vec (s2i v) (if positive then num else "-"^num)
3.77 + | setprod _ _ _ = raise (Converter "term is not a normed product")
3.78 +
3.79 + fun sum2vec (cplexSum ts) = foldl (fn (vec, t) => setprod vec true t) (empty_vector, ts)
3.80 + | sum2vec t = setprod empty_vector true t
3.81 +
3.82 + fun constrs2Ab j A b [] = (A, b)
3.83 + | constrs2Ab j A b ((_, cplexConstr (cplexLeq, (t1,t2)))::cs) =
3.84 + constrs2Ab (j+1) (set_vector A j (sum2vec t1)) (set_elem b j (num2str true t2)) cs
3.85 + | constrs2Ab j A b ((_, cplexConstr (cplexGeq, (t1,t2)))::cs) =
3.86 + constrs2Ab (j+1) (set_vector A j (sum2vec (neg_term t1))) (set_elem b j (num2str true (neg_term t2))) cs
3.87 + | constrs2Ab j A b ((_, cplexConstr (cplexEq, (t1,t2)))::cs) =
3.88 + constrs2Ab j A b ((None, cplexConstr (cplexLeq, (t1,t2)))::
3.89 + (None, cplexConstr (cplexGeq, (t1, t2)))::cs)
3.90 + | constrs2Ab _ _ _ _ = raise (Converter "no strict constraints allowed")
3.91 +
3.92 + val (A, b) = constrs2Ab 0 empty_matrix empty_vector constrs
3.93 +
3.94 + val (goal_maximize, goal_term) =
3.95 + case goal of
3.96 + (cplexMaximize t) => (true, t)
3.97 + | (cplexMinimize t) => (false, t)
3.98 + in
3.99 + (goal_maximize, sum2vec goal_term, A, b, (varcount, i2s, s2i))
3.100 + end
3.101 +
3.102 +fun convert_results (cplex.Optimal (opt, entries)) name2index =
3.103 + let
3.104 + fun setv (v, (name, value)) = (matrix_builder.set_elem v (name2index name) value)
3.105 + in
3.106 + (opt, foldl setv (matrix_builder.empty_vector, entries))
3.107 + end
3.108 + | convert_results _ _ = raise (Converter "No optimal result")
3.109 +
3.110 +
3.111 +end;
3.112 +
3.113 +structure SimpleMatrixBuilder : MATRIX_BUILDER =
3.114 +struct
3.115 +type vector = (int * string) list
3.116 +type matrix = (int * vector) list
3.117 +
3.118 +val empty_matrix = []
3.119 +val empty_vector = []
3.120 +
3.121 +exception Nat_expected of int;
3.122 +
3.123 +fun set_elem v i s = v @ [(i, s)]
3.124 +
3.125 +fun set_vector m i v = m @ [(i, v)]
3.126 +
3.127 +end;
3.128 +
3.129 +
3.130 +
3.131 +structure SimpleCplexMatrixConverter = MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = SimpleMatrixBuilder);
3.132 +
3.133 +
3.134 +
3.135 +
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/src/HOL/Matrix/ExactFloatingPoint.ML Fri Sep 03 17:10:36 2004 +0200
4.3 @@ -0,0 +1,215 @@
4.4 +structure ExactFloatingPoint :
4.5 +sig
4.6 + exception Destruct_floatstr of string
4.7 +
4.8 + val destruct_floatstr : (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string
4.9 +
4.10 + exception Floating_point of string
4.11 +
4.12 + type floatrep = IntInf.int * IntInf.int
4.13 + val approx_dec_by_bin : IntInf.int -> floatrep -> floatrep * floatrep
4.14 + val approx_decstr_by_bin : int -> string -> floatrep * floatrep
4.15 +end
4.16 +=
4.17 +struct
4.18 +
4.19 +fun fst (a,b) = a
4.20 +fun snd (a,b) = b
4.21 +
4.22 +val filter = List.filter;
4.23 +
4.24 +exception Destruct_floatstr of string;
4.25 +
4.26 +fun destruct_floatstr isDigit isExp number =
4.27 + let
4.28 + val numlist = filter (not o Char.isSpace) (String.explode number)
4.29 +
4.30 + fun countsigns ((#"+")::cs) = countsigns cs
4.31 + | countsigns ((#"-")::cs) =
4.32 + let
4.33 + val (positive, rest) = countsigns cs
4.34 + in
4.35 + (not positive, rest)
4.36 + end
4.37 + | countsigns cs = (true, cs)
4.38 +
4.39 + fun readdigits [] = ([], [])
4.40 + | readdigits (q as c::cs) =
4.41 + if (isDigit c) then
4.42 + let
4.43 + val (digits, rest) = readdigits cs
4.44 + in
4.45 + (c::digits, rest)
4.46 + end
4.47 + else
4.48 + ([], q)
4.49 +
4.50 + fun readfromexp_helper cs =
4.51 + let
4.52 + val (positive, rest) = countsigns cs
4.53 + val (digits, rest') = readdigits rest
4.54 + in
4.55 + case rest' of
4.56 + [] => (positive, digits)
4.57 + | _ => raise (Destruct_floatstr number)
4.58 + end
4.59 +
4.60 + fun readfromexp [] = (true, [])
4.61 + | readfromexp (c::cs) =
4.62 + if isExp c then
4.63 + readfromexp_helper cs
4.64 + else
4.65 + raise (Destruct_floatstr number)
4.66 +
4.67 + fun readfromdot [] = ([], readfromexp [])
4.68 + | readfromdot ((#".")::cs) =
4.69 + let
4.70 + val (digits, rest) = readdigits cs
4.71 + val exp = readfromexp rest
4.72 + in
4.73 + (digits, exp)
4.74 + end
4.75 + | readfromdot cs = readfromdot ((#".")::cs)
4.76 +
4.77 + val (positive, numlist) = countsigns numlist
4.78 + val (digits1, numlist) = readdigits numlist
4.79 + val (digits2, exp) = readfromdot numlist
4.80 + in
4.81 + (positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp))
4.82 + end
4.83 +
4.84 +type floatrep = IntInf.int * IntInf.int
4.85 +
4.86 +exception Floating_point of string;
4.87 +
4.88 +val ln2_10 = (Math.ln 10.0)/(Math.ln 2.0)
4.89 +
4.90 +fun intmul a b = IntInf.* (a,b)
4.91 +fun intsub a b = IntInf.- (a,b)
4.92 +fun intadd a b = IntInf.+ (a,b)
4.93 +fun intpow a b = IntInf.pow (a, IntInf.toInt b);
4.94 +fun intle a b = IntInf.<= (a, b);
4.95 +fun intless a b = IntInf.< (a, b);
4.96 +fun intneg a = IntInf.~ a;
4.97 +val zero = IntInf.fromInt 0;
4.98 +val one = IntInf.fromInt 1;
4.99 +val two = IntInf.fromInt 2;
4.100 +val ten = IntInf.fromInt 10;
4.101 +val five = IntInf.fromInt 5;
4.102 +
4.103 +fun find_most_significant q r =
4.104 + let
4.105 + fun int2real i =
4.106 + case Real.fromString (IntInf.toString i) of
4.107 + SOME r => r
4.108 + | NONE => raise (Floating_point "int2real")
4.109 + fun subtract (q, r) (q', r') =
4.110 + if intle r r' then
4.111 + (intsub q (intmul q' (intpow ten (intsub r' r))), r)
4.112 + else
4.113 + (intsub (intmul q (intpow ten (intsub r r'))) q', r')
4.114 + fun bin2dec d =
4.115 + if intle zero d then
4.116 + (intpow two d, zero)
4.117 + else
4.118 + (intpow five (intneg d), d)
4.119 +
4.120 + val L = IntInf.fromInt (Real.floor (int2real (IntInf.fromInt (IntInf.log2 q)) + (int2real r) * ln2_10))
4.121 + val L1 = intadd L one
4.122 +
4.123 + val (q1, r1) = subtract (q, r) (bin2dec L1)
4.124 + in
4.125 + if intle zero q1 then
4.126 + let
4.127 + val (q2, r2) = subtract (q, r) (bin2dec (intadd L1 one))
4.128 + in
4.129 + if intle zero q2 then
4.130 + raise (Floating_point "find_most_significant")
4.131 + else
4.132 + (L1, (q1, r1))
4.133 + end
4.134 + else
4.135 + let
4.136 + val (q0, r0) = subtract (q, r) (bin2dec L)
4.137 + in
4.138 + if intle zero q0 then
4.139 + (L, (q0, r0))
4.140 + else
4.141 + raise (Floating_point "find_most_significant")
4.142 + end
4.143 + end
4.144 +
4.145 +fun approx_dec_by_bin n (q,r) =
4.146 + let
4.147 + fun addseq acc d' [] = acc
4.148 + | addseq acc d' (d::ds) = addseq (intadd acc (intpow two (intsub d d'))) d' ds
4.149 +
4.150 + fun seq2bin [] = (zero, zero)
4.151 + | seq2bin (d::ds) = (intadd (addseq zero d ds) one, d)
4.152 +
4.153 + fun approx d_seq d0 precision (q,r) =
4.154 + if q = zero then
4.155 + let val x = seq2bin d_seq in
4.156 + (x, x)
4.157 + end
4.158 + else
4.159 + let
4.160 + val (d, (q', r')) = find_most_significant q r
4.161 + in
4.162 + if intless precision (intsub d0 d) then
4.163 + let
4.164 + val d' = intsub d0 precision
4.165 + val x1 = seq2bin (d_seq)
4.166 + val x2 = (intadd (intmul (fst x1) (intpow two (intsub (snd x1) d'))) one, d') (* = seq2bin (d'::d_seq) *)
4.167 + in
4.168 + (x1, x2)
4.169 + end
4.170 + else
4.171 + approx (d::d_seq) d0 precision (q', r')
4.172 + end
4.173 +
4.174 + fun approx_start precision (q, r) =
4.175 + if q = zero then
4.176 + ((zero, zero), (zero, zero))
4.177 + else
4.178 + let
4.179 + val (d, (q', r')) = find_most_significant q r
4.180 + in
4.181 + if intle precision zero then
4.182 + let
4.183 + val x1 = seq2bin [d]
4.184 + in
4.185 + if q' = zero then
4.186 + (x1, x1)
4.187 + else
4.188 + (x1, seq2bin [intadd d one])
4.189 + end
4.190 + else
4.191 + approx [d] d precision (q', r')
4.192 + end
4.193 + in
4.194 + if intle zero q then
4.195 + approx_start n (q,r)
4.196 + else
4.197 + let
4.198 + val ((a1,b1), (a2, b2)) = approx_start n (intneg q, r)
4.199 + in
4.200 + ((intneg a2, b2), (intneg a1, b1))
4.201 + end
4.202 + end
4.203 +
4.204 +fun approx_decstr_by_bin n decstr =
4.205 + let
4.206 + fun str2int s = case IntInf.fromString s of SOME x => x | NONE => zero
4.207 + fun signint p x = if p then x else intneg x
4.208 +
4.209 + val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr
4.210 + val s = IntInf.fromInt (size d2)
4.211 +
4.212 + val q = signint p (intadd (intmul (str2int d1) (intpow ten s)) (str2int d2))
4.213 + val r = intsub (signint ep (str2int e)) s
4.214 + in
4.215 + approx_dec_by_bin (IntInf.fromInt n) (q,r)
4.216 + end
4.217 +
4.218 +end;
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/src/HOL/Matrix/Float.thy Fri Sep 03 17:10:36 2004 +0200
5.3 @@ -0,0 +1,508 @@
5.4 +theory Float = Hyperreal:
5.5 +
5.6 +constdefs
5.7 + pow2 :: "int \<Rightarrow> real"
5.8 + "pow2 a == if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a))))"
5.9 + float :: "int * int \<Rightarrow> real"
5.10 + "float x == (real (fst x)) * (pow2 (snd x))"
5.11 +
5.12 +lemma pow2_0[simp]: "pow2 0 = 1"
5.13 +by (simp add: pow2_def)
5.14 +
5.15 +lemma pow2_1[simp]: "pow2 1 = 2"
5.16 +by (simp add: pow2_def)
5.17 +
5.18 +lemma pow2_neg: "pow2 x = inverse (pow2 (-x))"
5.19 +by (simp add: pow2_def)
5.20 +
5.21 +lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)"
5.22 +proof -
5.23 + have h: "! n. nat (2 + int n) - Suc 0 = nat (1 + int n)" by arith
5.24 + have g: "! a b. a - -1 = a + (1::int)" by arith
5.25 + have pos: "! n. pow2 (int n + 1) = 2 * pow2 (int n)"
5.26 + apply (auto, induct_tac n)
5.27 + apply (simp_all add: pow2_def)
5.28 + apply (rule_tac m1="2" and n1="nat (2 + int na)" in ssubst[OF realpow_num_eq_if])
5.29 + apply (auto simp add: h)
5.30 + apply arith
5.31 + done
5.32 + show ?thesis
5.33 + proof (induct a)
5.34 + case (1 n)
5.35 + from pos show ?case by (simp add: ring_eq_simps)
5.36 + next
5.37 + case (2 n)
5.38 + show ?case
5.39 + apply (auto)
5.40 + apply (subst pow2_neg[of "- int n"])
5.41 + apply (subst pow2_neg[of "-1 - int n"])
5.42 + apply (auto simp add: g pos)
5.43 + done
5.44 + qed
5.45 +qed
5.46 +
5.47 +lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)"
5.48 +proof (induct b)
5.49 + case (1 n)
5.50 + show ?case
5.51 + proof (induct n)
5.52 + case 0
5.53 + show ?case by simp
5.54 + next
5.55 + case (Suc m)
5.56 + show ?case by (auto simp add: ring_eq_simps pow2_add1 prems)
5.57 + qed
5.58 +next
5.59 + case (2 n)
5.60 + show ?case
5.61 + proof (induct n)
5.62 + case 0
5.63 + show ?case
5.64 + apply (auto)
5.65 + apply (subst pow2_neg[of "a + -1"])
5.66 + apply (subst pow2_neg[of "-1"])
5.67 + apply (simp)
5.68 + apply (insert pow2_add1[of "-a"])
5.69 + apply (simp add: ring_eq_simps)
5.70 + apply (subst pow2_neg[of "-a"])
5.71 + apply (simp)
5.72 + done
5.73 + case (Suc m)
5.74 + have a: "int m - (a + -2) = 1 + (int m - a + 1)" by arith
5.75 + have b: "int m - -2 = 1 + (int m + 1)" by arith
5.76 + show ?case
5.77 + apply (auto)
5.78 + apply (subst pow2_neg[of "a + (-2 - int m)"])
5.79 + apply (subst pow2_neg[of "-2 - int m"])
5.80 + apply (auto simp add: ring_eq_simps)
5.81 + apply (subst a)
5.82 + apply (subst b)
5.83 + apply (simp only: pow2_add1)
5.84 + apply (subst pow2_neg[of "int m - a + 1"])
5.85 + apply (subst pow2_neg[of "int m + 1"])
5.86 + apply auto
5.87 + apply (insert prems)
5.88 + apply (auto simp add: ring_eq_simps)
5.89 + done
5.90 + qed
5.91 +qed
5.92 +
5.93 +lemma "float (a, e) + float (b, e) = float (a + b, e)"
5.94 +by (simp add: float_def ring_eq_simps)
5.95 +
5.96 +constdefs
5.97 + int_of_real :: "real \<Rightarrow> int"
5.98 + "int_of_real x == SOME y. real y = x"
5.99 + real_is_int :: "real \<Rightarrow> bool"
5.100 + "real_is_int x == ? (u::int). x = real u"
5.101 +
5.102 +lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))"
5.103 +by (auto simp add: real_is_int_def int_of_real_def)
5.104 +
5.105 +lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \<Longrightarrow> float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)"
5.106 +by (simp add: float_def real_is_int_def2 pow2_add[symmetric])
5.107 +
5.108 +lemma pow2_int: "pow2 (int c) = (2::real)^c"
5.109 +by (simp add: pow2_def)
5.110 +
5.111 +lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)"
5.112 +by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric])
5.113 +
5.114 +lemma real_is_int_real[simp]: "real_is_int (real (x::int))"
5.115 +by (auto simp add: real_is_int_def int_of_real_def)
5.116 +
5.117 +lemma int_of_real_real[simp]: "int_of_real (real x) = x"
5.118 +by (simp add: int_of_real_def)
5.119 +
5.120 +lemma real_int_of_real[simp]: "real_is_int x \<Longrightarrow> real (int_of_real x) = x"
5.121 +by (auto simp add: int_of_real_def real_is_int_def)
5.122 +
5.123 +lemma real_is_int_add_int_of_real: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a+b)) = (int_of_real a) + (int_of_real b)"
5.124 +by (auto simp add: int_of_real_def real_is_int_def)
5.125 +
5.126 +lemma real_is_int_add[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a+b)"
5.127 +apply (subst real_is_int_def2)
5.128 +apply (simp add: real_is_int_add_int_of_real real_int_of_real)
5.129 +done
5.130 +
5.131 +lemma int_of_real_sub: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> (int_of_real (a-b)) = (int_of_real a) - (int_of_real b)"
5.132 +by (auto simp add: int_of_real_def real_is_int_def)
5.133 +
5.134 +lemma real_is_int_sub[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a-b)"
5.135 +apply (subst real_is_int_def2)
5.136 +apply (simp add: int_of_real_sub real_int_of_real)
5.137 +done
5.138 +
5.139 +lemma real_is_int_rep: "real_is_int x \<Longrightarrow> ?! (a::int). real a = x"
5.140 +by (auto simp add: real_is_int_def)
5.141 +
5.142 +lemma int_of_real_mult:
5.143 + assumes "real_is_int a" "real_is_int b"
5.144 + shows "(int_of_real (a*b)) = (int_of_real a) * (int_of_real b)"
5.145 +proof -
5.146 + from prems have a: "?! (a'::int). real a' = a" by (rule_tac real_is_int_rep, auto)
5.147 + from prems have b: "?! (b'::int). real b' = b" by (rule_tac real_is_int_rep, auto)
5.148 + from a obtain a'::int where a':"a = real a'" by auto
5.149 + from b obtain b'::int where b':"b = real b'" by auto
5.150 + have r: "real a' * real b' = real (a' * b')" by auto
5.151 + show ?thesis
5.152 + apply (simp add: a' b')
5.153 + apply (subst r)
5.154 + apply (simp only: int_of_real_real)
5.155 + done
5.156 +qed
5.157 +
5.158 +lemma real_is_int_mult[simp]: "real_is_int a \<Longrightarrow> real_is_int b \<Longrightarrow> real_is_int (a*b)"
5.159 +apply (subst real_is_int_def2)
5.160 +apply (simp add: int_of_real_mult)
5.161 +done
5.162 +
5.163 +lemma real_is_int_0[simp]: "real_is_int (0::real)"
5.164 +by (simp add: real_is_int_def int_of_real_def)
5.165 +
5.166 +lemma real_is_int_1[simp]: "real_is_int (1::real)"
5.167 +proof -
5.168 + have "real_is_int (1::real) = real_is_int(real (1::int))" by auto
5.169 + also have "\<dots> = True" by (simp only: real_is_int_real)
5.170 + ultimately show ?thesis by auto
5.171 +qed
5.172 +
5.173 +lemma real_is_int_n1: "real_is_int (-1::real)"
5.174 +proof -
5.175 + have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
5.176 + also have "\<dots> = True" by (simp only: real_is_int_real)
5.177 + ultimately show ?thesis by auto
5.178 +qed
5.179 +
5.180 +lemma real_is_int_number_of[simp]: "real_is_int ((number_of::bin\<Rightarrow>real) x)"
5.181 +proof -
5.182 + have neg1: "real_is_int (-1::real)"
5.183 + proof -
5.184 + have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto
5.185 + also have "\<dots> = True" by (simp only: real_is_int_real)
5.186 + ultimately show ?thesis by auto
5.187 + qed
5.188 +
5.189 + {
5.190 + fix x::int
5.191 + have "!! y. real_is_int ((number_of::bin\<Rightarrow>real) (Abs_Bin x))"
5.192 + apply (simp add: number_of_eq)
5.193 + apply (subst Abs_Bin_inverse)
5.194 + apply (simp add: Bin_def)
5.195 + apply (induct x)
5.196 + apply (induct_tac n)
5.197 + apply (simp)
5.198 + apply (simp)
5.199 + apply (induct_tac n)
5.200 + apply (simp add: neg1)
5.201 + proof -
5.202 + fix n :: nat
5.203 + assume rn: "(real_is_int (of_int (- (int (Suc n)))))"
5.204 + have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp
5.205 + show "real_is_int (of_int (- (int (Suc (Suc n)))))"
5.206 + apply (simp only: s of_int_add)
5.207 + apply (rule real_is_int_add)
5.208 + apply (simp add: neg1)
5.209 + apply (simp only: rn)
5.210 + done
5.211 + qed
5.212 + }
5.213 + note Abs_Bin = this
5.214 + {
5.215 + fix x :: bin
5.216 + have "? u. x = Abs_Bin u"
5.217 + apply (rule exI[where x = "Rep_Bin x"])
5.218 + apply (simp add: Rep_Bin_inverse)
5.219 + done
5.220 + }
5.221 + then obtain u::int where "x = Abs_Bin u" by auto
5.222 + with Abs_Bin show ?thesis by auto
5.223 +qed
5.224 +
5.225 +lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)"
5.226 +by (simp add: int_of_real_def)
5.227 +
5.228 +lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)"
5.229 +proof -
5.230 + have 1: "(1::real) = real (1::int)" by auto
5.231 + show ?thesis by (simp only: 1 int_of_real_real)
5.232 +qed
5.233 +
5.234 +lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b"
5.235 +proof -
5.236 + have "real_is_int (number_of b)" by simp
5.237 + then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep)
5.238 + then obtain u::int where u:"number_of b = real u" by auto
5.239 + have "number_of b = real ((number_of b)::int)"
5.240 + by (simp add: number_of_eq real_of_int_def)
5.241 + have ub: "number_of b = real ((number_of b)::int)"
5.242 + by (simp add: number_of_eq real_of_int_def)
5.243 + from uu u ub have unb: "u = number_of b"
5.244 + by blast
5.245 + have "int_of_real (number_of b) = u" by (simp add: u)
5.246 + with unb show ?thesis by simp
5.247 +qed
5.248 +
5.249 +lemma float_transfer_even: "even a \<Longrightarrow> float (a, b) = float (a div 2, b+1)"
5.250 + apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified])
5.251 + apply (simp_all add: pow2_def even_def real_is_int_def ring_eq_simps)
5.252 + apply (auto)
5.253 +proof -
5.254 + fix q::int
5.255 + have a:"b - (-1\<Colon>int) = (1\<Colon>int) + b" by arith
5.256 + show "(float (q, (b - (-1\<Colon>int)))) = (float (q, ((1\<Colon>int) + b)))"
5.257 + by (simp add: a)
5.258 +qed
5.259 +
5.260 +consts
5.261 + norm_float :: "int*int \<Rightarrow> int*int"
5.262 +
5.263 +lemma int_div_zdiv: "int (a div b) = (int a) div (int b)"
5.264 +apply (subst split_div, auto)
5.265 +apply (subst split_zdiv, auto)
5.266 +apply (rule_tac a="int (b * i) + int j" and b="int b" and r="int j" and r'=ja in IntDiv.unique_quotient)
5.267 +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
5.268 +done
5.269 +
5.270 +lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)"
5.271 +apply (subst split_mod, auto)
5.272 +apply (subst split_zmod, auto)
5.273 +apply (rule_tac a="int (b * i) + int j" and b="int b" and q="int i" and q'=ia in IntDiv.unique_remainder)
5.274 +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat)
5.275 +done
5.276 +
5.277 +lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a"
5.278 +by arith
5.279 +
5.280 +lemma terminating_norm_float: "\<forall>a. (a::int) \<noteq> 0 \<and> even a \<longrightarrow> a \<noteq> 0 \<and> \<bar>a div 2\<bar> < \<bar>a\<bar>"
5.281 +apply (auto)
5.282 +apply (rule abs_div_2_less)
5.283 +apply (auto)
5.284 +done
5.285 +
5.286 +ML {* simp_depth_limit := 2 *}
5.287 +recdef norm_float "measure (% (a,b). nat (abs a))"
5.288 + "norm_float (a,b) = (if (a \<noteq> 0) & (even a) then norm_float (a div 2, b+1) else (if a=0 then (0,0) else (a,b)))"
5.289 +(hints simp: terminating_norm_float)
5.290 +ML {* simp_depth_limit := 1000 *}
5.291 +
5.292 +
5.293 +lemma norm_float: "float x = float (norm_float x)"
5.294 +proof -
5.295 + {
5.296 + fix a b :: int
5.297 + have norm_float_pair: "float (a,b) = float (norm_float (a,b))"
5.298 + proof (induct a b rule: norm_float.induct)
5.299 + case (1 u v)
5.300 + show ?case
5.301 + proof cases
5.302 + assume u: "u \<noteq> 0 \<and> even u"
5.303 + with prems have ind: "float (u div 2, v + 1) = float (norm_float (u div 2, v + 1))" by auto
5.304 + with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even)
5.305 + then show ?thesis
5.306 + apply (subst norm_float.simps)
5.307 + apply (simp add: ind)
5.308 + done
5.309 + next
5.310 + assume "~(u \<noteq> 0 \<and> even u)"
5.311 + then show ?thesis
5.312 + by (simp add: prems float_def)
5.313 + qed
5.314 + qed
5.315 + }
5.316 + note helper = this
5.317 + have "? a b. x = (a,b)" by auto
5.318 + then obtain a b where "x = (a, b)" by blast
5.319 + then show ?thesis by (simp only: helper)
5.320 +qed
5.321 +
5.322 +lemma pow2_int: "pow2 (int n) = 2^n"
5.323 + by (simp add: pow2_def)
5.324 +
5.325 +lemma float_add:
5.326 + "float (a1, e1) + float (a2, e2) =
5.327 + (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1)
5.328 + else float (a1*2^(nat (e1-e2))+a2, e2))"
5.329 + apply (simp add: float_def ring_eq_simps)
5.330 + apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric])
5.331 + done
5.332 +
5.333 +lemma float_mult:
5.334 + "float (a1, e1) * float (a2, e2) =
5.335 + (float (a1 * a2, e1 + e2))"
5.336 + by (simp add: float_def pow2_add)
5.337 +
5.338 +lemma float_minus:
5.339 + "- (float (a,b)) = float (-a, b)"
5.340 + by (simp add: float_def)
5.341 +
5.342 +lemma zero_less_pow2:
5.343 + "0 < pow2 x"
5.344 +proof -
5.345 + {
5.346 + fix y
5.347 + have "0 <= y \<Longrightarrow> 0 < pow2 y"
5.348 + by (induct y, induct_tac n, simp_all add: pow2_add)
5.349 + }
5.350 + note helper=this
5.351 + show ?thesis
5.352 + apply (case_tac "0 <= x")
5.353 + apply (simp add: helper)
5.354 + apply (subst pow2_neg)
5.355 + apply (simp add: helper)
5.356 + done
5.357 +qed
5.358 +
5.359 +lemma zero_le_float:
5.360 + "(0 <= float (a,b)) = (0 <= a)"
5.361 + apply (auto simp add: float_def)
5.362 + apply (auto simp add: zero_le_mult_iff zero_less_pow2)
5.363 + apply (insert zero_less_pow2[of b])
5.364 + apply (simp_all)
5.365 + done
5.366 +
5.367 +lemma float_abs:
5.368 + "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))"
5.369 + apply (auto simp add: abs_if)
5.370 + apply (simp_all add: zero_le_float[symmetric, of a b] float_minus)
5.371 + done
5.372 +
5.373 +lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1"
5.374 + by auto
5.375 +
5.376 +lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)"
5.377 + by simp
5.378 +
5.379 +lemma add_right_zero: "a + 0 = (a::'a::comm_monoid_add)"
5.380 + by simp
5.381 +
5.382 +lemma mult_left_one: "1 * a = (a::'a::semiring_1)"
5.383 + by simp
5.384 +
5.385 +lemma mult_right_one: "a * 1 = (a::'a::semiring_1)"
5.386 + by simp
5.387 +
5.388 +lemma int_pow_0: "(a::int)^(Numeral0) = 1"
5.389 + by simp
5.390 +
5.391 +lemma int_pow_1: "(a::int)^(Numeral1) = a"
5.392 + by simp
5.393 +
5.394 +lemma zero_eq_Numeral0_nring: "(0::'a::number_ring) = Numeral0"
5.395 + by simp
5.396 +
5.397 +lemma one_eq_Numeral1_nring: "(1::'a::number_ring) = Numeral1"
5.398 + by simp
5.399 +
5.400 +lemma zero_eq_Numeral0_nat: "(0::nat) = Numeral0"
5.401 + by simp
5.402 +
5.403 +lemma one_eq_Numeral1_nat: "(1::nat) = Numeral1"
5.404 + by simp
5.405 +
5.406 +lemma zpower_Pls: "(z::int)^Numeral0 = Numeral1"
5.407 + by simp
5.408 +
5.409 +lemma zpower_Min: "(z::int)^((-1)::nat) = Numeral1"
5.410 +proof -
5.411 + have 1:"((-1)::nat) = 0"
5.412 + by simp
5.413 + show ?thesis by (simp add: 1)
5.414 +qed
5.415 +
5.416 +lemma fst_cong: "a=a' \<Longrightarrow> fst (a,b) = fst (a',b)"
5.417 + by simp
5.418 +
5.419 +lemma snd_cong: "b=b' \<Longrightarrow> snd (a,b) = snd (a,b')"
5.420 + by simp
5.421 +
5.422 +lemma lift_bool: "x \<Longrightarrow> x=True"
5.423 + by simp
5.424 +
5.425 +lemma nlift_bool: "~x \<Longrightarrow> x=False"
5.426 + by simp
5.427 +
5.428 +lemma not_false_eq_true: "(~ False) = True" by simp
5.429 +
5.430 +lemma not_true_eq_false: "(~ True) = False" by simp
5.431 +
5.432 +
5.433 +lemmas binarith =
5.434 + Pls_0_eq Min_1_eq
5.435 + bin_pred_Pls bin_pred_Min bin_pred_1 bin_pred_0
5.436 + bin_succ_Pls bin_succ_Min bin_succ_1 bin_succ_0
5.437 + bin_add_Pls bin_add_Min bin_add_BIT_0 bin_add_BIT_10
5.438 + bin_add_BIT_11 bin_minus_Pls bin_minus_Min bin_minus_1
5.439 + bin_minus_0 bin_mult_Pls bin_mult_Min bin_mult_1 bin_mult_0
5.440 + bin_add_Pls_right bin_add_Min_right
5.441 +
5.442 +thm eq_number_of_eq
5.443 +
5.444 +lemma int_eq_number_of_eq: "(((number_of v)::int)=(number_of w)) = iszero ((number_of (bin_add v (bin_minus w)))::int)"
5.445 + by simp
5.446 +
5.447 +lemma int_iszero_number_of_Pls: "iszero (Numeral0::int)"
5.448 + by (simp only: iszero_number_of_Pls)
5.449 +
5.450 +lemma int_nonzero_number_of_Min: "~(iszero ((-1)::int))"
5.451 + by simp
5.452 +thm iszero_number_of_1
5.453 +
5.454 +lemma int_iszero_number_of_0: "iszero ((number_of (w BIT False))::int) = iszero ((number_of w)::int)"
5.455 + by simp
5.456 +
5.457 +lemma int_iszero_number_of_1: "\<not> iszero ((number_of (w BIT True))::int)"
5.458 + by simp
5.459 +
5.460 +lemma int_less_number_of_eq_neg: "(((number_of x)::int) < number_of y) = neg ((number_of (bin_add x (bin_minus y)))::int)"
5.461 + by simp
5.462 +
5.463 +lemma int_not_neg_number_of_Pls: "\<not> (neg (Numeral0::int))"
5.464 + by simp
5.465 +
5.466 +lemma int_neg_number_of_Min: "neg (-1::int)"
5.467 + by simp
5.468 +
5.469 +lemma int_neg_number_of_BIT: "neg ((number_of (w BIT x))::int) = neg ((number_of w)::int)"
5.470 + by simp
5.471 +
5.472 +lemma int_le_number_of_eq: "(((number_of x)::int) \<le> number_of y) = (\<not> neg ((number_of (bin_add y (bin_minus x)))::int))"
5.473 + by simp
5.474 +
5.475 +lemmas intarithrel =
5.476 + (*abs_zero abs_one*)
5.477 + int_eq_number_of_eq
5.478 + lift_bool[OF int_iszero_number_of_Pls] nlift_bool[OF int_nonzero_number_of_Min] int_iszero_number_of_0
5.479 + lift_bool[OF int_iszero_number_of_1] int_less_number_of_eq_neg nlift_bool[OF int_not_neg_number_of_Pls] lift_bool[OF int_neg_number_of_Min]
5.480 + int_neg_number_of_BIT int_le_number_of_eq
5.481 +
5.482 +lemma int_number_of_add_sym: "((number_of v)::int) + number_of w = number_of (bin_add v w)"
5.483 + by simp
5.484 +
5.485 +lemma int_number_of_diff_sym: "((number_of v)::int) - number_of w = number_of (bin_add v (bin_minus w))"
5.486 + by simp
5.487 +
5.488 +lemma int_number_of_mult_sym: "((number_of v)::int) * number_of w = number_of (bin_mult v w)"
5.489 + by simp
5.490 +
5.491 +lemma int_number_of_minus_sym: "- ((number_of v)::int) = number_of (bin_minus v)"
5.492 + by simp
5.493 +
5.494 +lemmas intarith = int_number_of_add_sym int_number_of_minus_sym int_number_of_diff_sym int_number_of_mult_sym
5.495 +
5.496 +lemmas natarith = (*zero_eq_Numeral0_nat one_eq_Numeral1_nat*) add_nat_number_of
5.497 + diff_nat_number_of mult_nat_number_of eq_nat_number_of less_nat_number_of
5.498 +
5.499 +lemmas powerarith = nat_number_of zpower_number_of_even
5.500 + zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring]
5.501 + zpower_Pls zpower_Min
5.502 +
5.503 +lemmas floatarith[simplified norm_0_1] = float_add float_mult float_minus float_abs zero_le_float
5.504 +
5.505 +lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false
5.506 +
5.507 +(* needed for the verifying code generator *)
5.508 +lemmas arith_no_let = arith[simplified Let_def]
5.509 +
5.510 +end
5.511 +
6.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
6.2 +++ b/src/HOL/Matrix/FloatArith.ML Fri Sep 03 17:10:36 2004 +0200
6.3 @@ -0,0 +1,56 @@
6.4 +structure FloatArith =
6.5 +struct
6.6 +
6.7 +type float = IntInf.int * IntInf.int
6.8 +
6.9 +val izero = IntInf.fromInt 0
6.10 +fun imul a b = IntInf.* (a,b)
6.11 +fun isub a b = IntInf.- (a,b)
6.12 +fun iadd a b = IntInf.+ (a,b)
6.13 +
6.14 +val floatzero = (izero, izero)
6.15 +
6.16 +fun positive_part (a,b) =
6.17 + (if IntInf.< (a,izero) then izero else a, b)
6.18 +
6.19 +fun negative_part (a,b) =
6.20 + (if IntInf.< (a,izero) then a else izero, b)
6.21 +
6.22 +fun is_negative (a,b) =
6.23 + if IntInf.< (a, izero) then true else false
6.24 +
6.25 +fun is_positive (a,b) =
6.26 + if IntInf.< (izero, a) then true else false
6.27 +
6.28 +fun is_zero (a,b) =
6.29 + if a = izero then true else false
6.30 +
6.31 +fun ipow2 a = IntInf.pow ((IntInf.fromInt 2), IntInf.toInt a)
6.32 +
6.33 +fun add (a1, b1) (a2, b2) =
6.34 + if b1 < b2 then
6.35 + (iadd a1 (imul a2 (ipow2 (isub b2 b1))), b1)
6.36 + else
6.37 + (iadd (imul a1 (ipow2 (isub b1 b2))) a2, b2)
6.38 +
6.39 +fun sub (a1, b1) (a2, b2) =
6.40 + if b1 < b2 then
6.41 + (isub a1 (imul a2 (ipow2 (isub b2 b1))), b1)
6.42 + else
6.43 + (isub (imul a1 (ipow2 (isub b1 b2))) a2, b2)
6.44 +
6.45 +fun neg (a, b) = (IntInf.~ a, b)
6.46 +
6.47 +fun is_equal a b = is_zero (sub a b)
6.48 +
6.49 +fun is_less a b = is_negative (sub a b)
6.50 +
6.51 +fun max a b = if is_less a b then b else a
6.52 +
6.53 +fun min a b = if is_less a b then a else b
6.54 +
6.55 +fun abs a = if is_negative a then neg a else a
6.56 +
6.57 +fun mul (a1, b1) (a2, b2) = (imul a1 a2, iadd b1 b2)
6.58 +
6.59 +end;
7.1 --- a/src/HOL/Matrix/LinProg.thy Fri Sep 03 10:27:05 2004 +0200
7.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
7.3 @@ -1,79 +0,0 @@
7.4 -(* Title: HOL/Matrix/LinProg.thy
7.5 - ID: $Id$
7.6 - Author: Steven Obua
7.7 -*)
7.8 -
7.9 -theory LinProg = Matrix:
7.10 -
7.11 -lemma linprog_by_duality_approx_general:
7.12 - assumes
7.13 - fmuladdprops:
7.14 - "! a b c d. a <= b & c <= d \<longrightarrow> fadd a c <= fadd b d"
7.15 - "! c a b. 0 <= c & a <= b \<longrightarrow> fmul c a <= fmul c b"
7.16 - "! a. fmul 0 a = 0"
7.17 - "! a. fmul a 0 = 0"
7.18 - "fadd 0 0 = 0"
7.19 - "associative fadd"
7.20 - "commutative fadd"
7.21 - "associative fmul"
7.22 - "distributive fmul fadd"
7.23 - and specificprops:
7.24 - "mult_matrix fmul fadd (combine_matrix fadd A dA) x <= (b::('a::{order, zero}) matrix)"
7.25 - "mult_matrix fmul fadd y A = c"
7.26 - "0 <= y"
7.27 - shows
7.28 - "combine_matrix fadd (mult_matrix fmul fadd c x) (mult_matrix fmul fadd (mult_matrix fmul fadd y dA) x)
7.29 - <= mult_matrix fmul fadd y b"
7.30 -proof -
7.31 - let ?mul = "mult_matrix fmul fadd"
7.32 - let ?add = "combine_matrix fadd"
7.33 - let ?t1 = "?mul y (?mul (?add A dA) x)"
7.34 - have a: "?t1 <= ?mul y b" by (rule le_left_mult, simp_all!)
7.35 - have b: "?t1 = ?mul (?mul y (?add A dA)) x" by (simp! add: mult_matrix_assoc_simple[THEN sym])
7.36 - have assoc: "associative ?add" by (simp! add: combine_matrix_assoc)
7.37 - have r_distr: "r_distributive ?mul ?add"
7.38 - apply (rule r_distributive_matrix)
7.39 - by (simp! add: distributive_def)+
7.40 - have l_distr: "l_distributive ?mul ?add"
7.41 - apply (rule l_distributive_matrix)
7.42 - by (simp! add: distributive_def)+
7.43 - have c:"?mul y (?add A dA) = ?add (?mul y A) (?mul y dA)"
7.44 - by (insert r_distr, simp add: r_distributive_def)
7.45 - have d:"?mul (?add (?mul y A) (?mul y dA)) x = ?add (?mul (?mul y A) x) (?mul (?mul y dA) x)"
7.46 - by (insert l_distr, simp add: l_distributive_def)
7.47 - have e:"\<dots> = ?add (?mul c x) (?mul (?mul y dA) x)" by (simp!)
7.48 - from a b c d e show "?add (?mul c x) (?mul (?mul y dA) x) <= ?mul y b" by simp
7.49 -qed
7.50 -
7.51 -lemma linprog_by_duality_approx:
7.52 - assumes
7.53 - "(A + dA) * x <= (b::('a::lordered_ring) matrix)"
7.54 - "y * A = c"
7.55 - "0 <= y"
7.56 - shows
7.57 - "c * x + (y * dA) * x <= y * b"
7.58 -apply (simp add: times_matrix_def plus_matrix_def)
7.59 -apply (rule linprog_by_duality_approx_general)
7.60 -apply (simp_all)
7.61 -apply (simp_all add: associative_def add_assoc mult_assoc)
7.62 -apply (simp_all add: commutative_def add_commute)
7.63 -apply (simp_all add: distributive_def l_distributive_def r_distributive_def left_distrib right_distrib)
7.64 -apply (simp_all! add: plus_matrix_def times_matrix_def)
7.65 -apply (simp add: add_mono)
7.66 -apply (simp add: mult_left_mono)
7.67 -done
7.68 -
7.69 -lemma linprog_by_duality:
7.70 - assumes
7.71 - "A * x <= (b::('a::lordered_ring) matrix)"
7.72 - "y * A = c"
7.73 - "0 <= y"
7.74 - shows
7.75 - "c * x <= y * b"
7.76 -proof -
7.77 - have a:"(A + 0) * x <= b" by (simp!)
7.78 - have b:"c * x + (y*0)*x <= y * b" by (rule linprog_by_duality_approx, simp_all!)
7.79 - show "c * x <= y*b" by (insert b, simp)
7.80 -qed
7.81 -
7.82 -end
8.1 --- a/src/HOL/Matrix/Matrix.thy Fri Sep 03 10:27:05 2004 +0200
8.2 +++ b/src/HOL/Matrix/Matrix.thy Fri Sep 03 17:10:36 2004 +0200
8.3 @@ -21,7 +21,10 @@
8.4 times_matrix_def: "A * B == mult_matrix (op *) (op +) A B"
8.5
8.6 lemma is_meet_combine_matrix_meet: "is_meet (combine_matrix meet)"
8.7 -by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
8.8 + by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le)
8.9 +
8.10 +lemma is_join_combine_matrix_join: "is_join (combine_matrix join)"
8.11 + by (simp_all add: is_join_def le_matrix_def join_left_le join_right_le join_imp_le)
8.12
8.13 instance matrix :: (lordered_ab_group) lordered_ab_group_meet
8.14 proof
8.15 @@ -284,7 +287,20 @@
8.16 apply (auto)
8.17 done
8.18
8.19 +lemma Rep_minus[simp]: "Rep_matrix (-(A::_::lordered_ab_group)) x y = - (Rep_matrix A x y)"
8.20 + by (simp add: minus_matrix_def)
8.21
8.22 +lemma join_matrix: "join (A::('a::lordered_ring) matrix) B = combine_matrix join A B"
8.23 + apply (insert join_unique[of "(combine_matrix join)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_join_combine_matrix_join])
8.24 + apply (simp)
8.25 + done
8.26
8.27 +lemma meet_matrix: "meet (A::('a::lordered_ring) matrix) B = combine_matrix meet A B"
8.28 + apply (insert meet_unique[of "(combine_matrix meet)::('a matrix \<Rightarrow> 'a matrix \<Rightarrow> 'a matrix)", simplified is_meet_combine_matrix_meet])
8.29 + apply (simp)
8.30 + done
8.31 +
8.32 +lemma Rep_abs[simp]: "Rep_matrix (abs (A::_::lordered_ring)) x y = abs (Rep_matrix A x y)"
8.33 + by (simp add: abs_lattice join_matrix)
8.34
8.35 end
9.1 --- a/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 10:27:05 2004 +0200
9.2 +++ b/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 17:10:36 2004 +0200
9.3 @@ -1392,7 +1392,44 @@
9.4 apply (rule le_foldseq)
9.5 by (auto)
9.6
9.7 +lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
9.8 +lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
9.9 +
9.10 lemma singleton_matrix_le[simp]: "(singleton_matrix j i a <= singleton_matrix j i b) = (a <= (b::_::order))"
9.11 by (auto simp add: le_matrix_def)
9.12
9.13 +lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
9.14 + apply (auto)
9.15 + apply (simp add: le_matrix_def)
9.16 + apply (drule_tac j=j and i=i in spec2)
9.17 + apply (simp)
9.18 + apply (simp add: le_matrix_def)
9.19 + done
9.20 +
9.21 +lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
9.22 + apply (auto)
9.23 + apply (simp add: le_matrix_def)
9.24 + apply (drule_tac j=j and i=i in spec2)
9.25 + apply (simp)
9.26 + apply (simp add: le_matrix_def)
9.27 + done
9.28 +
9.29 +lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
9.30 + apply (auto simp add: le_matrix_def neg_def)
9.31 + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
9.32 + apply (auto)
9.33 + done
9.34 +
9.35 +lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
9.36 + apply (auto simp add: le_matrix_def neg_def)
9.37 + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
9.38 + apply (auto)
9.39 + done
9.40 +
9.41 +lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))"
9.42 + apply (auto simp add: le_matrix_def neg_def)
9.43 + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
9.44 + apply (auto)
9.45 + done
9.46 +
9.47 end
10.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
10.2 +++ b/src/HOL/Matrix/MatrixLP.ML Fri Sep 03 17:10:36 2004 +0200
10.3 @@ -0,0 +1,83 @@
10.4 +use "MatrixLP_gensimp.ML";
10.5 +
10.6 +signature MATRIX_LP =
10.7 +sig
10.8 + val lp_dual_estimate : string -> int -> thm
10.9 + val simplify : thm -> thm
10.10 +end
10.11 +
10.12 +structure MatrixLP : MATRIX_LP =
10.13 +struct
10.14 +
10.15 +val _ = SimprocsCodegen.use_simprocs_code sg geninput
10.16 +
10.17 +val simp_le_spmat = simp_SparseMatrix_le_spmat_'_mult__'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_bool'
10.18 +val simp_add_spmat = simp_SparseMatrix_add_spmat_'_mult__'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
10.19 +val simp_diff_spmat = simp_SparseMatrix_diff_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
10.20 +val simp_abs_spmat = simp_SparseMatrix_abs_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
10.21 +val simp_mult_spmat = simp_SparseMatrix_mult_spmat_'List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real''''_List_list_'_mult__'nat'_'List_list_'_mult__'nat'_'RealDef_real'''''
10.22 +val simp_sorted_sparse_matrix = simp_SparseMatrix_sorted_sparse_matrix
10.23 +val simp_sorted_spvec = simp_SparseMatrix_sorted_spvec
10.24 +
10.25 +fun single_arg simp ct = snd (simp (snd (Thm.dest_comb ct)))
10.26 +
10.27 +fun double_arg simp ct =
10.28 + let
10.29 + val (a, y) = Thm.dest_comb ct
10.30 + val (_, x) = Thm.dest_comb a
10.31 + in
10.32 + snd (simp x y)
10.33 + end
10.34 +
10.35 +fun spmc ct =
10.36 + (case term_of ct of
10.37 + ((Const ("SparseMatrix.le_spmat", _)) $ _) =>
10.38 + single_arg simp_le_spmat ct
10.39 + | ((Const ("SparseMatrix.add_spmat", _)) $ _) =>
10.40 + single_arg simp_add_spmat ct
10.41 + | ((Const ("SparseMatrix.diff_spmat", _)) $ _ $ _) =>
10.42 + double_arg simp_diff_spmat ct
10.43 + | ((Const ("SparseMatrix.abs_spmat", _)) $ _) =>
10.44 + single_arg simp_abs_spmat ct
10.45 + | ((Const ("SparseMatrix.mult_spmat", _)) $ _ $ _) =>
10.46 + double_arg simp_mult_spmat ct
10.47 + | ((Const ("SparseMatrix.sorted_sparse_matrix", _)) $ _) =>
10.48 + single_arg simp_sorted_sparse_matrix ct
10.49 + | ((Const ("SparseMatrix.sorted_spvec", ty)) $ _) =>
10.50 + single_arg simp_sorted_spvec ct
10.51 + | _ => raise Fail_conv)
10.52 +
10.53 +fun lp_dual_estimate lptfile prec =
10.54 + let
10.55 + val th = inst_real (thm "SparseMatrix.spm_linprog_dual_estimate_1")
10.56 + val sg = sign_of (theory "MatrixLP")
10.57 + val (y, (A1, A2), (c1, c2), b, r) =
10.58 + let
10.59 + open fspmlp
10.60 + val l = load lptfile prec false
10.61 + in
10.62 + (y l, A l, c l, b l, r l)
10.63 + end
10.64 + fun var s x = (cterm_of sg (Var ((s,0), FloatSparseMatrixBuilder.real_spmatT)), x)
10.65 + val th = Thm.instantiate ([], [var "A1" A1, var "A2" A2, var "y" y, var "c1" c1, var "c2" c2, var "r" r, var "b" b]) th
10.66 + in
10.67 + th
10.68 + end
10.69 +
10.70 +fun simplify th =
10.71 + let
10.72 + val simp_th = botc spmc (cprop_of th)
10.73 + val th = strip_shyps (equal_elim simp_th th)
10.74 + fun removeTrue th = removeTrue (implies_elim th TrueI) handle _ => th
10.75 + in
10.76 + removeTrue th
10.77 + end
10.78 +
10.79 +fun float2real (x,y) = (Real.fromInt x) * (Math.pow (2.0, Real.fromInt y))
10.80 +
10.81 +end
10.82 +
10.83 +val result = ref ([]:thm list)
10.84 +
10.85 +fun run c = timeit (fn () => (result := [MatrixLP.simplify c]; ()))
10.86 +
11.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
11.2 +++ b/src/HOL/Matrix/MatrixLP.thy Fri Sep 03 17:10:36 2004 +0200
11.3 @@ -0,0 +1,4 @@
11.4 +theory MatrixLP = Float + SparseMatrix:
11.5 +
11.6 +end
11.7 +
12.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
12.2 +++ b/src/HOL/Matrix/MatrixLP_gensimp.ML Fri Sep 03 17:10:36 2004 +0200
12.3 @@ -0,0 +1,35 @@
12.4 +open BasisLibrary;
12.5 +
12.6 +use "Cplex.ML";
12.7 +use "CplexMatrixConverter.ML";
12.8 +use "ExactFloatingPoint.ML";
12.9 +use "FloatArith.ML";
12.10 +use "FloatSparseMatrixBuilder.ML";
12.11 +use "fspmlp.ML";
12.12 +use "codegen_prep.ML";
12.13 +use "eq_codegen.ML";
12.14 +use "conv.ML";
12.15 +
12.16 +
12.17 +val th = theory "MatrixLP"
12.18 +val sg = sign_of th
12.19 +
12.20 +fun prep x = standard (mk_meta_eq x)
12.21 +
12.22 +fun inst_real thm = standard (Thm.instantiate ([(fst (hd (term_tvars (prop_of thm))),
12.23 + ctyp_of sg HOLogic.realT)], []) thm);
12.24 +
12.25 +val spm_lp_dual_estimate_simp =
12.26 + (thms "Float.arith_no_let") @
12.27 + (map inst_real (thms "SparseMatrix.sparse_row_matrix_arith_simps")) @
12.28 + (thms "SparseMatrix.sorted_sp_simps") @
12.29 + [thm "Product_Type.fst_conv", thm "Product_Type.snd_conv"] @
12.30 + (thms "SparseMatrix.boolarith")
12.31 +
12.32 +val simp_with = Simplifier.simplify (HOL_basic_ss addsimps [thm "SparseMatrix.if_case_eq", thm "Float.norm_0_1"])
12.33 +
12.34 +val spm_lp_dual_estimate_simp = map simp_with spm_lp_dual_estimate_simp
12.35 +
12.36 +val geninput = codegen_prep.prepare_thms spm_lp_dual_estimate_simp
12.37 +
12.38 +val _ = SimprocsCodegen.use_simprocs_code sg geninput
13.1 --- a/src/HOL/Matrix/ROOT.ML Fri Sep 03 10:27:05 2004 +0200
13.2 +++ b/src/HOL/Matrix/ROOT.ML Fri Sep 03 17:10:36 2004 +0200
13.3 @@ -1,10 +1,1 @@
13.4 -(* Title: HOL/Matrix/ROOT.ML
13.5 - ID: $Id$
13.6 - Author: Steven Obua
13.7 - License: 2004 Technische Universität München
13.8 -
13.9 -Theory of matrices with an application of matrix theory to linear
13.10 -programming.
13.11 -*)
13.12 -
13.13 -use_thy "LinProg";
13.14 +use_thy "MatrixLP"
14.1 --- a/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 10:27:05 2004 +0200
14.2 +++ b/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 17:10:36 2004 +0200
14.3 @@ -94,9 +94,62 @@
14.4 done
14.5
14.6 consts
14.7 + abs_spvec :: "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
14.8 + minus_spvec :: "('a::lordered_ring) spvec \<Rightarrow> 'a spvec"
14.9 smult_spvec :: "('a::lordered_ring) \<Rightarrow> 'a spvec \<Rightarrow> 'a spvec"
14.10 addmult_spvec :: "('a::lordered_ring) * 'a spvec * 'a spvec \<Rightarrow> 'a spvec"
14.11
14.12 +primrec
14.13 + "minus_spvec [] = []"
14.14 + "minus_spvec (a#as) = (fst a, -(snd a))#(minus_spvec as)"
14.15 +
14.16 +primrec
14.17 + "abs_spvec [] = []"
14.18 + "abs_spvec (a#as) = (fst a, abs (snd a))#(abs_spvec as)"
14.19 +
14.20 +lemma sparse_row_vector_minus:
14.21 + "sparse_row_vector (minus_spvec v) = - (sparse_row_vector v)"
14.22 + apply (induct v)
14.23 + apply (simp_all add: sparse_row_vector_cons)
14.24 + apply (simp add: Rep_matrix_inject[symmetric])
14.25 + apply (rule ext)+
14.26 + apply simp
14.27 + done
14.28 +
14.29 +lemma sparse_row_vector_abs:
14.30 + "sorted_spvec v \<Longrightarrow> sparse_row_vector (abs_spvec v) = abs (sparse_row_vector v)"
14.31 + apply (induct v)
14.32 + apply (simp_all add: sparse_row_vector_cons)
14.33 + apply (frule_tac sorted_spvec_cons1, simp)
14.34 + apply (simp only: Rep_matrix_inject[symmetric])
14.35 + apply (rule ext)+
14.36 + apply auto
14.37 + apply (subgoal_tac "Rep_matrix (sparse_row_vector list) 0 a = 0")
14.38 + apply (simp)
14.39 + apply (rule sorted_sparse_row_vector_zero)
14.40 + apply auto
14.41 + done
14.42 +
14.43 +lemma sorted_spvec_minus_spvec:
14.44 + "sorted_spvec v \<Longrightarrow> sorted_spvec (minus_spvec v)"
14.45 + apply (induct v)
14.46 + apply (simp)
14.47 + apply (frule sorted_spvec_cons1, simp)
14.48 + apply (simp add: sorted_spvec.simps)
14.49 + apply (case_tac list)
14.50 + apply (simp_all)
14.51 + done
14.52 +
14.53 +lemma sorted_spvec_abs_spvec:
14.54 + "sorted_spvec v \<Longrightarrow> sorted_spvec (abs_spvec v)"
14.55 + apply (induct v)
14.56 + apply (simp)
14.57 + apply (frule sorted_spvec_cons1, simp)
14.58 + apply (simp add: sorted_spvec.simps)
14.59 + apply (case_tac list)
14.60 + apply (simp_all)
14.61 + done
14.62 +
14.63 defs
14.64 smult_spvec_def: "smult_spvec y arr == map (% a. (fst a, y * snd a)) arr"
14.65
14.66 @@ -542,9 +595,6 @@
14.67 else
14.68 (le_spvec(snd a, snd b) & le_spmat (as, bs))))"
14.69
14.70 -lemma spec2: "! j i. P j i \<Longrightarrow> P j i" by blast
14.71 -lemma neg_imp: "(\<not> Q \<longrightarrow> \<not> P) \<Longrightarrow> P \<longrightarrow> Q" by blast
14.72 -
14.73 constdefs
14.74 disj_matrices :: "('a::zero) matrix \<Rightarrow> 'a matrix \<Rightarrow> bool"
14.75 "disj_matrices A B == (! j i. (Rep_matrix A j i \<noteq> 0) \<longrightarrow> (Rep_matrix B j i = 0)) & (! j i. (Rep_matrix B j i \<noteq> 0) \<longrightarrow> (Rep_matrix A j i = 0))"
14.76 @@ -597,22 +647,6 @@
14.77 (B + A <= C) = (A <= C & (B::('a::lordered_ab_group) matrix) <= 0)"
14.78 by (auto simp add: disj_matrices_add[of B A 0 C,simplified] disj_matrices_commute)
14.79
14.80 -lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))"
14.81 - apply (auto)
14.82 - apply (simp add: le_matrix_def)
14.83 - apply (drule_tac j=j and i=i in spec2)
14.84 - apply (simp)
14.85 - apply (simp add: le_matrix_def)
14.86 - done
14.87 -
14.88 -lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)"
14.89 - apply (auto)
14.90 - apply (simp add: le_matrix_def)
14.91 - apply (drule_tac j=j and i=i in spec2)
14.92 - apply (simp)
14.93 - apply (simp add: le_matrix_def)
14.94 - done
14.95 -
14.96 lemma disj_sparse_row_singleton: "i <= j \<Longrightarrow> sorted_spvec((j,y)#v) \<Longrightarrow> disj_matrices (sparse_row_vector v) (singleton_matrix 0 i x)"
14.97 apply (simp add: disj_matrices_def)
14.98 apply (rule conjI)
14.99 @@ -641,27 +675,6 @@
14.100 lemma disj_singleton_matrices[simp]: "disj_matrices (singleton_matrix j i x) (singleton_matrix u v y) = (j \<noteq> u | i \<noteq> v | x = 0 | y = 0)"
14.101 by (auto simp add: disj_matrices_def)
14.102
14.103 -lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \<longrightarrow> (sorted_spvec b) \<longrightarrow> (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)"
14.104 - apply (rule le_spvec.induct)
14.105 - apply (simp_all add: sorted_spvec_cons1 disj_matrices_add_le_zero disj_matrices_add_zero_le disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.106 - apply (rule conjI, intro strip)
14.107 - apply (simp add: sorted_spvec_cons1)
14.108 - apply (subst disj_matrices_add_x_le)
14.109 - apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
14.110 - apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.111 - apply (simp, blast)
14.112 - apply (intro strip, rule conjI, intro strip)
14.113 - apply (simp add: sorted_spvec_cons1)
14.114 - apply (subst disj_matrices_add_le_x)
14.115 - apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_sparse_row_singleton[OF less_imp_le] disj_matrices_commute disj_matrices_x_add)
14.116 - apply (blast)
14.117 - apply (intro strip)
14.118 - apply (simp add: sorted_spvec_cons1)
14.119 - apply (case_tac "a=aa", simp_all)
14.120 - apply (subst disj_matrices_add)
14.121 - apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.122 - done
14.123 -
14.124 lemma disj_move_sparse_vec_mat[simplified disj_matrices_commute]:
14.125 "j <= a \<Longrightarrow> sorted_spvec((a,c)#as) \<Longrightarrow> disj_matrices (move_matrix (sparse_row_vector b) (int j) i) (sparse_row_matrix as)"
14.126 apply (auto simp add: neg_def disj_matrices_def)
14.127 @@ -685,24 +698,28 @@
14.128 apply (rule nrows, rule order_trans[of _ 1], simp, drule nrows_notzero, drule less_le_trans[OF _ nrows_spvec], arith)+
14.129 done
14.130
14.131 -lemma move_matrix_le_zero[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))"
14.132 - apply (auto simp add: le_matrix_def neg_def)
14.133 - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
14.134 - apply (auto)
14.135 +lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \<longrightarrow> (sorted_spvec b) \<longrightarrow> (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)"
14.136 + apply (rule le_spvec.induct)
14.137 + apply (simp_all add: sorted_spvec_cons1 disj_matrices_add_le_zero disj_matrices_add_zero_le
14.138 + disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.139 + apply (rule conjI, intro strip)
14.140 + apply (simp add: sorted_spvec_cons1)
14.141 + apply (subst disj_matrices_add_x_le)
14.142 + apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute)
14.143 + apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.144 + apply (simp, blast)
14.145 + apply (intro strip, rule conjI, intro strip)
14.146 + apply (simp add: sorted_spvec_cons1)
14.147 + apply (subst disj_matrices_add_le_x)
14.148 + apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_sparse_row_singleton[OF less_imp_le] disj_matrices_commute disj_matrices_x_add)
14.149 + apply (blast)
14.150 + apply (intro strip)
14.151 + apply (simp add: sorted_spvec_cons1)
14.152 + apply (case_tac "a=aa", simp_all)
14.153 + apply (subst disj_matrices_add)
14.154 + apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute)
14.155 done
14.156
14.157 -lemma move_matrix_zero_le[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)"
14.158 - apply (auto simp add: le_matrix_def neg_def)
14.159 - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
14.160 - apply (auto)
14.161 - done
14.162 -
14.163 -lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \<Longrightarrow> 0 <= i \<Longrightarrow> (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))"
14.164 - apply (auto simp add: le_matrix_def neg_def)
14.165 - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2)
14.166 - apply (auto)
14.167 - done
14.168 -
14.169 lemma le_spvec_empty2_sparse_row[rule_format]: "(sorted_spvec b) \<longrightarrow> (le_spvec (b,[]) = (sparse_row_vector b <= 0))"
14.170 apply (induct b)
14.171 apply (simp_all add: sorted_spvec_cons1)
14.172 @@ -752,39 +769,169 @@
14.173 apply (simp add: sorted_spvec_cons1 le_spvec_iff_sparse_row_le)
14.174 done
14.175
14.176 -term smult_spvec
14.177 -term addmult_spvec
14.178 -term add_spvec
14.179 -term mult_spvec_spmat
14.180 -term mult_spmat
14.181 -term add_spmat
14.182 -term le_spvec
14.183 -term le_spmat
14.184 -term sorted_spvec
14.185 -term sorted_spmat
14.186 +ML {* simp_depth_limit := 999 *}
14.187
14.188 -thm sparse_row_mult_spmat
14.189 -thm sparse_row_add_spmat
14.190 -thm le_spmat_iff_sparse_row_le
14.191 +consts
14.192 + abs_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
14.193 + minus_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat"
14.194
14.195 -thm sorted_spvec_mult_spmat
14.196 -thm sorted_spmat_mult_spmat
14.197 -thm sorted_spvec_add_spmat
14.198 -thm sorted_spmat_add_spmat
14.199 +primrec
14.200 + "abs_spmat [] = []"
14.201 + "abs_spmat (a#as) = (fst a, abs_spvec (snd a))#(abs_spmat as)"
14.202
14.203 -thm smult_spvec_empty
14.204 -thm smult_spvec_cons
14.205 -thm addmult_spvec.simps
14.206 -thm add_spvec.simps
14.207 -thm add_spmat.simps
14.208 -thm mult_spvec_spmat.simps
14.209 -thm mult_spmat.simps
14.210 -thm le_spvec.simps
14.211 -thm le_spmat.simps
14.212 -thm sorted_spvec.simps
14.213 -thm sorted_spmat.simps
14.214 +primrec
14.215 + "minus_spmat [] = []"
14.216 + "minus_spmat (a#as) = (fst a, minus_spvec (snd a))#(minus_spmat as)"
14.217 +
14.218 +lemma sparse_row_matrix_minus:
14.219 + "sparse_row_matrix (minus_spmat A) = - (sparse_row_matrix A)"
14.220 + apply (induct A)
14.221 + apply (simp_all add: sparse_row_vector_minus sparse_row_matrix_cons)
14.222 + apply (subst Rep_matrix_inject[symmetric])
14.223 + apply (rule ext)+
14.224 + apply simp
14.225 + done
14.226 +
14.227 +lemma Rep_sparse_row_vector_zero: "x \<noteq> 0 \<Longrightarrow> Rep_matrix (sparse_row_vector v) x y = 0"
14.228 +proof -
14.229 + assume x:"x \<noteq> 0"
14.230 + have r:"nrows (sparse_row_vector v) <= Suc 0" by (rule nrows_spvec)
14.231 + show ?thesis
14.232 + apply (rule nrows)
14.233 + apply (subgoal_tac "Suc 0 <= x")
14.234 + apply (insert r)
14.235 + apply (simp only:)
14.236 + apply (insert x)
14.237 + apply arith
14.238 + done
14.239 +qed
14.240 +
14.241 +lemma sparse_row_matrix_abs:
14.242 + "sorted_spvec A \<Longrightarrow> sorted_spmat A \<Longrightarrow> sparse_row_matrix (abs_spmat A) = abs (sparse_row_matrix A)"
14.243 + apply (induct A)
14.244 + apply (simp_all add: sparse_row_vector_abs sparse_row_matrix_cons)
14.245 + apply (frule_tac sorted_spvec_cons1, simp)
14.246 + apply (subst Rep_matrix_inject[symmetric])
14.247 + apply (rule ext)+
14.248 + apply auto
14.249 + apply (case_tac "x=a")
14.250 + apply (simp)
14.251 + apply (subst sorted_sparse_row_matrix_zero)
14.252 + apply auto
14.253 + apply (subst Rep_sparse_row_vector_zero)
14.254 + apply (simp_all add: neg_def)
14.255 + done
14.256 +
14.257 +lemma sorted_spvec_minus_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (minus_spmat A)"
14.258 + apply (induct A)
14.259 + apply (simp)
14.260 + apply (frule sorted_spvec_cons1, simp)
14.261 + apply (simp add: sorted_spvec.simps)
14.262 + apply (case_tac list)
14.263 + apply (simp_all)
14.264 + done
14.265 +
14.266 +lemma sorted_spvec_abs_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec (abs_spmat A)"
14.267 + apply (induct A)
14.268 + apply (simp)
14.269 + apply (frule sorted_spvec_cons1, simp)
14.270 + apply (simp add: sorted_spvec.simps)
14.271 + apply (case_tac list)
14.272 + apply (simp_all)
14.273 + done
14.274 +
14.275 +lemma sorted_spmat_minus_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (minus_spmat A)"
14.276 + apply (induct A)
14.277 + apply (simp_all add: sorted_spvec_minus_spvec)
14.278 + done
14.279 +
14.280 +lemma sorted_spmat_abs_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat (abs_spmat A)"
14.281 + apply (induct A)
14.282 + apply (simp_all add: sorted_spvec_abs_spvec)
14.283 + done
14.284 +
14.285 +constdefs
14.286 + diff_spmat :: "('a::lordered_ring) spmat \<Rightarrow> 'a spmat \<Rightarrow> 'a spmat"
14.287 + "diff_spmat A B == add_spmat (A, minus_spmat B)"
14.288 +
14.289 +lemma sorted_spmat_diff_spmat: "sorted_spmat A \<Longrightarrow> sorted_spmat B \<Longrightarrow> sorted_spmat (diff_spmat A B)"
14.290 + by (simp add: diff_spmat_def sorted_spmat_minus_spmat sorted_spmat_add_spmat)
14.291 +
14.292 +lemma sorted_spvec_diff_spmat: "sorted_spvec A \<Longrightarrow> sorted_spvec B \<Longrightarrow> sorted_spvec (diff_spmat A B)"
14.293 + by (simp add: diff_spmat_def sorted_spvec_minus_spmat sorted_spvec_add_spmat)
14.294 +
14.295 +lemma sparse_row_diff_spmat: "sparse_row_matrix (diff_spmat A B ) = (sparse_row_matrix A) - (sparse_row_matrix B)"
14.296 + by (simp add: diff_spmat_def sparse_row_add_spmat sparse_row_matrix_minus)
14.297 +
14.298 +constdefs
14.299 + sorted_sparse_matrix :: "'a spmat \<Rightarrow> bool"
14.300 + "sorted_sparse_matrix A == (sorted_spvec A) & (sorted_spmat A)"
14.301 +
14.302 +lemma sorted_sparse_matrix_imp_spvec: "sorted_sparse_matrix A \<Longrightarrow> sorted_spvec A"
14.303 + by (simp add: sorted_sparse_matrix_def)
14.304 +
14.305 +lemma sorted_sparse_matrix_imp_spmat: "sorted_sparse_matrix A \<Longrightarrow> sorted_spmat A"
14.306 + by (simp add: sorted_sparse_matrix_def)
14.307 +
14.308 +lemmas sparse_row_matrix_op_simps =
14.309 + sorted_sparse_matrix_imp_spmat sorted_sparse_matrix_imp_spvec
14.310 + sparse_row_add_spmat sorted_spvec_add_spmat sorted_spmat_add_spmat
14.311 + sparse_row_diff_spmat sorted_spvec_diff_spmat sorted_spmat_diff_spmat
14.312 + sparse_row_matrix_minus sorted_spvec_minus_spmat sorted_spmat_minus_spmat
14.313 + sparse_row_mult_spmat sorted_spvec_mult_spmat sorted_spmat_mult_spmat
14.314 + sparse_row_matrix_abs sorted_spvec_abs_spmat sorted_spmat_abs_spmat
14.315 + le_spmat_iff_sparse_row_le
14.316 +
14.317 +lemma zero_eq_Numeral0: "(0::_::number_ring) = Numeral0" by simp
14.318 +
14.319 +lemmas sparse_row_matrix_arith_simps[simplified zero_eq_Numeral0] =
14.320 + mult_spmat.simps mult_spvec_spmat.simps
14.321 + addmult_spvec.simps
14.322 + smult_spvec_empty smult_spvec_cons
14.323 + add_spmat.simps add_spvec.simps
14.324 + minus_spmat.simps minus_spvec.simps
14.325 + abs_spmat.simps abs_spvec.simps
14.326 + diff_spmat_def
14.327 + le_spmat.simps le_spvec.simps
14.328 +
14.329 +lemmas sorted_sp_simps =
14.330 + sorted_spvec.simps
14.331 + sorted_spmat.simps
14.332 + sorted_sparse_matrix_def
14.333 +
14.334 +lemma bool1: "(\<not> True) = False" by blast
14.335 +lemma bool2: "(\<not> False) = True" by blast
14.336 +lemma bool3: "((P\<Colon>bool) \<and> True) = P" by blast
14.337 +lemma bool4: "(True \<and> (P\<Colon>bool)) = P" by blast
14.338 +lemma bool5: "((P\<Colon>bool) \<and> False) = False" by blast
14.339 +lemma bool6: "(False \<and> (P\<Colon>bool)) = False" by blast
14.340 +lemma bool7: "((P\<Colon>bool) \<or> True) = True" by blast
14.341 +lemma bool8: "(True \<or> (P\<Colon>bool)) = True" by blast
14.342 +lemma bool9: "((P\<Colon>bool) \<or> False) = P" by blast
14.343 +lemma bool10: "(False \<or> (P\<Colon>bool)) = P" by blast
14.344 +lemmas boolarith = bool1 bool2 bool3 bool4 bool5 bool6 bool7 bool8 bool9 bool10
14.345 +
14.346 +lemma if_case_eq: "(if b then x else y) = (case b of True => x | False => y)" by simp
14.347 +
14.348 +lemma spm_linprog_dual_estimate_1:
14.349 + assumes
14.350 + "sorted_sparse_matrix A1"
14.351 + "sorted_sparse_matrix A2"
14.352 + "sorted_sparse_matrix c1"
14.353 + "sorted_sparse_matrix c2"
14.354 + "sorted_sparse_matrix y"
14.355 + "sorted_spvec b"
14.356 + "sorted_spvec r"
14.357 + "le_spmat ([], y)"
14.358 + "A * x \<le> sparse_row_matrix (b::('a::lordered_ring) spmat)"
14.359 + "sparse_row_matrix A1 <= A"
14.360 + "A <= sparse_row_matrix A2"
14.361 + "sparse_row_matrix c1 <= c"
14.362 + "c <= sparse_row_matrix c2"
14.363 + "abs x \<le> sparse_row_matrix r"
14.364 + shows
14.365 + "c * x \<le> sparse_row_matrix (add_spmat (mult_spmat y b, mult_spmat (add_spmat (add_spmat (mult_spmat y (diff_spmat A2 A1),
14.366 + abs_spmat (diff_spmat (mult_spmat y A1) c1)), diff_spmat c2 c1)) r))"
14.367 + by (insert prems, simp add: sparse_row_matrix_op_simps linprog_dual_estimate_1[where A=A])
14.368
14.369 end
14.370 -
14.371 -
14.372 -
15.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
15.2 +++ b/src/HOL/Matrix/codegen_prep.ML Fri Sep 03 17:10:36 2004 +0200
15.3 @@ -0,0 +1,110 @@
15.4 +structure codegen_prep =
15.5 +struct
15.6 +
15.7 +exception Prepare_thms of string;
15.8 +
15.9 +fun is_meta_eq th =
15.10 + let
15.11 + fun check ((Const ("==", _)) $ _ $ _) = true
15.12 + | check _ = false
15.13 + in
15.14 + check (concl_of th)
15.15 + end
15.16 +
15.17 +fun printty ty =
15.18 + let
15.19 + fun immerse s [] = []
15.20 + | immerse s [x] = [x]
15.21 + | immerse s (x::xs) = x::s::(immerse s xs)
15.22 +
15.23 + fun py t =
15.24 + let
15.25 + val (name, args) = dest_Type t
15.26 + val args' = map printty args
15.27 + in
15.28 + concat (immerse "_" (name::args'))
15.29 + end
15.30 +
15.31 + val (args, res) = strip_type ty
15.32 + val tystrs = map py (args@[res])
15.33 +
15.34 + val s = "'"^(concat (immerse "_" tystrs))^"'"
15.35 + in
15.36 + s
15.37 + end
15.38 +
15.39 +fun head_name th =
15.40 + let val s =
15.41 + let
15.42 + val h = fst (strip_comb (hd (snd (strip_comb (concl_of th)))))
15.43 + val n = fst (dest_Const h)
15.44 + val ty = snd (dest_Const h)
15.45 + val hn = fst (dest_Const h)
15.46 + in
15.47 + hn^"_"^(printty ty) handle _ => (writeln ("warning: polymorphic "^hn); hn)
15.48 + end
15.49 + in
15.50 + s
15.51 + end
15.52 +
15.53 +val mangle_id =
15.54 + let
15.55 + fun tr #"=" = "_eq_"
15.56 + | tr #"." = "_"
15.57 + | tr #" " = "_"
15.58 + | tr #"<" = "_le_"
15.59 + | tr #">" = "_ge_"
15.60 + | tr #"(" = "_bro_"
15.61 + | tr #")" = "_brc_"
15.62 + | tr #"+" = "_plus_"
15.63 + | tr #"*" = "_mult_"
15.64 + | tr #"-" = "_minus_"
15.65 + | tr #"|" = "_or_"
15.66 + | tr #"&" = "_and_"
15.67 + | tr x = Char.toString x
15.68 + fun m s = "simp_"^(String.translate tr s)
15.69 + in
15.70 + m
15.71 + end
15.72 +
15.73 +fun group f [] = []
15.74 + | group f (x::xs) =
15.75 + let
15.76 + val g = group f xs
15.77 + val key = f x
15.78 + in
15.79 + case assoc (g, key) of
15.80 + None => (key, [x])::g
15.81 + | Some l => overwrite (g, (key, x::l))
15.82 + end
15.83 +
15.84 +fun prepare_thms ths =
15.85 + let
15.86 + val ths = (filter is_meta_eq ths) @ (map (standard o mk_meta_eq) (filter (not o is_meta_eq) ths))
15.87 + val _ = if forall Thm.no_prems ths then () else raise (Prepare_thms "premisse found")
15.88 + val thmgroups = group head_name ths
15.89 + val test_clashgroups = group (fn (a,b) => mangle_id a) thmgroups
15.90 + val _ = if (length thmgroups <> length test_clashgroups) then raise (Prepare_thms "clash after name mangling") else ()
15.91 +
15.92 + fun prep (name, ths) =
15.93 + let
15.94 + val m = mangle_id name
15.95 +
15.96 + in
15.97 + (m, true, ths)
15.98 + end
15.99 +
15.100 + val thmgroups = map prep thmgroups
15.101 + in
15.102 + (thmgroups)
15.103 + end
15.104 +
15.105 +fun writestr filename s =
15.106 + let
15.107 + val f = TextIO.openOut filename
15.108 + val _ = TextIO.output(f, s)
15.109 + val _ = TextIO.closeOut f
15.110 + in
15.111 + ()
15.112 + end
15.113 +end
16.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
16.2 +++ b/src/HOL/Matrix/conv.ML Fri Sep 03 17:10:36 2004 +0200
16.3 @@ -0,0 +1,26 @@
16.4 +exception Fail_conv;
16.5 +
16.6 +fun orelsec conv1 conv2 ct = conv1 ct handle Fail_conv => conv2 ct
16.7 +
16.8 +val allc = Thm.reflexive
16.9 +
16.10 +fun thenc conv1 conv2 ct =
16.11 + let
16.12 + fun rhs_of t = snd (Thm.dest_comb (strip_imp_concl (cprop_of t)))
16.13 +
16.14 + val eq = conv1 ct
16.15 + in
16.16 + Thm.transitive eq (conv2 (rhs_of eq))
16.17 + end
16.18 +
16.19 +fun subc conv ct =
16.20 + case term_of ct of
16.21 + _ $ _ =>
16.22 + let
16.23 + val (ct1, ct2) = Thm.dest_comb ct
16.24 + in
16.25 + Thm.combination (conv ct1) (conv ct2)
16.26 + end
16.27 + | _ => allc ct
16.28 +
16.29 +fun botc conv ct = thenc (subc (botc conv)) (orelsec conv allc) ct
16.30 \ No newline at end of file
17.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
17.2 +++ b/src/HOL/Matrix/eq_codegen.ML Fri Sep 03 17:10:36 2004 +0200
17.3 @@ -0,0 +1,493 @@
17.4 +fun inst_cterm inst ct = fst (Drule.dest_equals
17.5 + (Thm.cprop_of (Thm.instantiate inst (reflexive ct))));
17.6 +fun tyinst_cterm tyinst = inst_cterm (tyinst, []);
17.7 +
17.8 +val bla = ref ([] : term list);
17.9 +
17.10 +(******************************************************)
17.11 +(* Code generator for equational proofs *)
17.12 +(******************************************************)
17.13 +fun my_mk_meta_eq thm =
17.14 + let
17.15 + val (_, eq) = Thm.dest_comb (cprop_of thm);
17.16 + val (ct, rhs) = Thm.dest_comb eq;
17.17 + val (_, lhs) = Thm.dest_comb ct
17.18 + in Thm.implies_elim (Drule.instantiate' [Some (ctyp_of_term lhs)]
17.19 + [Some lhs, Some rhs] eq_reflection) thm
17.20 + end;
17.21 +
17.22 +structure SimprocsCodegen =
17.23 +struct
17.24 +
17.25 +val simp_thms = ref ([] : thm list);
17.26 +
17.27 +fun parens b = if b then Pretty.enclose "(" ")" else Pretty.block;
17.28 +
17.29 +fun gen_mk_val f xs ps = Pretty.block ([Pretty.str "val ",
17.30 + f (length xs > 1) (flat
17.31 + (separate [Pretty.str ",", Pretty.brk 1] (map (single o Pretty.str) xs))),
17.32 + Pretty.str " =", Pretty.brk 1] @ ps @ [Pretty.str ";"]);
17.33 +
17.34 +val mk_val = gen_mk_val parens;
17.35 +val mk_vall = gen_mk_val (K (Pretty.enclose "[" "]"));
17.36 +
17.37 +fun rename s = if s mem ThmDatabase.ml_reserved then s ^ "'" else s;
17.38 +
17.39 +fun mk_decomp_name (Var ((s, i), _)) = rename (if i=0 then s else s ^ string_of_int i)
17.40 + | mk_decomp_name (Const (s, _)) = rename (Codegen.mk_id (Sign.base_name s))
17.41 + | mk_decomp_name _ = "ct";
17.42 +
17.43 +fun decomp_term_code cn ((vs, bs, ps), (v, t)) =
17.44 + if exists (equal t o fst) bs then (vs, bs, ps)
17.45 + else (case t of
17.46 + Var _ => (vs, bs @ [(t, v)], ps)
17.47 + | Const _ => (vs, if cn then bs @ [(t, v)] else bs, ps)
17.48 + | Bound _ => (vs, bs, ps)
17.49 + | Abs (s, T, t) =>
17.50 + let
17.51 + val v1 = variant vs s;
17.52 + val v2 = variant (v1 :: vs) (mk_decomp_name t)
17.53 + in
17.54 + decomp_term_code cn ((v1 :: v2 :: vs,
17.55 + bs @ [(Free (s, T), v1)],
17.56 + ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_abs", Pretty.brk 1,
17.57 + Pretty.str "None", Pretty.brk 1, Pretty.str v]]), (v2, t))
17.58 + end
17.59 + | t $ u =>
17.60 + let
17.61 + val v1 = variant vs (mk_decomp_name t);
17.62 + val v2 = variant (v1 :: vs) (mk_decomp_name u);
17.63 + val (vs', bs', ps') = decomp_term_code cn ((v1 :: v2 :: vs, bs,
17.64 + ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_comb", Pretty.brk 1,
17.65 + Pretty.str v]]), (v1, t));
17.66 + val (vs'', bs'', ps'') = decomp_term_code cn ((vs', bs', ps'), (v2, u))
17.67 + in
17.68 + if bs'' = bs then (vs, bs, ps) else (vs'', bs'', ps'')
17.69 + end);
17.70 +
17.71 +val strip_tv = implode o tl o explode;
17.72 +
17.73 +fun mk_decomp_tname (TVar ((s, i), _)) =
17.74 + strip_tv ((if i=0 then s else s ^ string_of_int i) ^ "T")
17.75 + | mk_decomp_tname (Type (s, _)) = Codegen.mk_id (Sign.base_name s) ^ "T"
17.76 + | mk_decomp_tname _ = "cT";
17.77 +
17.78 +fun decomp_type_code ((vs, bs, ps), (v, TVar (ixn, _))) =
17.79 + if exists (equal ixn o fst) bs then (vs, bs, ps)
17.80 + else (vs, bs @ [(ixn, v)], ps)
17.81 + | decomp_type_code ((vs, bs, ps), (v, Type (_, Ts))) =
17.82 + let
17.83 + val vs' = variantlist (map mk_decomp_tname Ts, vs);
17.84 + val (vs'', bs', ps') =
17.85 + foldl decomp_type_code ((vs @ vs', bs, ps @
17.86 + [mk_vall vs' [Pretty.str "Thm.dest_ctyp", Pretty.brk 1,
17.87 + Pretty.str v]]), vs' ~~ Ts)
17.88 + in
17.89 + if bs' = bs then (vs, bs, ps) else (vs'', bs', ps')
17.90 + end;
17.91 +
17.92 +fun gen_mk_bindings s dest decomp ((vs, bs, ps), (v, x)) =
17.93 + let
17.94 + val s' = variant vs s;
17.95 + val (vs', bs', ps') = decomp ((s' :: vs, bs, ps @
17.96 + [mk_val [s'] (dest v)]), (s', x))
17.97 + in
17.98 + if bs' = bs then (vs, bs, ps) else (vs', bs', ps')
17.99 + end;
17.100 +
17.101 +val mk_term_bindings = gen_mk_bindings "ct"
17.102 + (fn s => [Pretty.str "cprop_of", Pretty.brk 1, Pretty.str s])
17.103 + (decomp_term_code true);
17.104 +
17.105 +val mk_type_bindings = gen_mk_bindings "cT"
17.106 + (fn s => [Pretty.str "Thm.ctyp_of_term", Pretty.brk 1, Pretty.str s])
17.107 + decomp_type_code;
17.108 +
17.109 +fun pretty_pattern b (Const (s, _)) = Pretty.block [Pretty.str "Const",
17.110 + Pretty.brk 1, Pretty.str ("(\"" ^ s ^ "\", _)")]
17.111 + | pretty_pattern b (t as _ $ _) = parens b
17.112 + (flat (separate [Pretty.str " $", Pretty.brk 1]
17.113 + (map (single o pretty_pattern true) (op :: (strip_comb t)))))
17.114 + | pretty_pattern b _ = Pretty.str "_";
17.115 +
17.116 +fun term_consts' t = foldl_aterms
17.117 + (fn (cs, c as Const _) => c ins cs | (cs, _) => cs) ([], t);
17.118 +
17.119 +fun mk_apps s b p [] = p
17.120 + | mk_apps s b p (q :: qs) =
17.121 + mk_apps s b (parens (b orelse not (null qs))
17.122 + [Pretty.str s, Pretty.brk 1, p, Pretty.brk 1, q]) qs;
17.123 +
17.124 +fun mk_refleq eq ct = mk_val [eq] [Pretty.str ("Thm.reflexive " ^ ct)];
17.125 +
17.126 +fun mk_tyinst ((s, i), s') =
17.127 + Pretty.block [Pretty.str ("((" ^ quote s ^ ","), Pretty.brk 1,
17.128 + Pretty.str (string_of_int i ^ "),"), Pretty.brk 1,
17.129 + Pretty.str (s' ^ ")")];
17.130 +
17.131 +fun inst_ty b ty_bs t s = (case term_tvars t of
17.132 + [] => Pretty.str s
17.133 + | Ts => parens b [Pretty.str "tyinst_cterm", Pretty.brk 1,
17.134 + Pretty.list "[" "]" (map (fn (ixn, _) => mk_tyinst
17.135 + (ixn, the (assoc (ty_bs, ixn)))) Ts),
17.136 + Pretty.brk 1, Pretty.str s]);
17.137 +
17.138 +fun mk_cterm_code b ty_bs ts xs (vals, t $ u) =
17.139 + let
17.140 + val (vals', p1) = mk_cterm_code true ty_bs ts xs (vals, t);
17.141 + val (vals'', p2) = mk_cterm_code true ty_bs ts xs (vals', u)
17.142 + in
17.143 + (vals'', parens b [Pretty.str "Thm.capply", Pretty.brk 1,
17.144 + p1, Pretty.brk 1, p2])
17.145 + end
17.146 + | mk_cterm_code b ty_bs ts xs (vals, Abs (s, T, t)) =
17.147 + let
17.148 + val u = Free (s, T);
17.149 + val Some s' = assoc (ts, u);
17.150 + val p = Pretty.str s';
17.151 + val (vals', p') = mk_cterm_code true ty_bs ts (p :: xs)
17.152 + (if null (typ_tvars T) then vals
17.153 + else vals @ [(u, (("", s'), [mk_val [s'] [inst_ty true ty_bs u s']]))], t)
17.154 + in (vals',
17.155 + parens b [Pretty.str "Thm.cabs", Pretty.brk 1, p, Pretty.brk 1, p'])
17.156 + end
17.157 + | mk_cterm_code b ty_bs ts xs (vals, Bound i) = (vals, nth_elem (i, xs))
17.158 + | mk_cterm_code b ty_bs ts xs (vals, t) = (case assoc (vals, t) of
17.159 + None =>
17.160 + let val Some s = assoc (ts, t)
17.161 + in (if is_Const t andalso not (null (term_tvars t)) then
17.162 + vals @ [(t, (("", s), [mk_val [s] [inst_ty true ty_bs t s]]))]
17.163 + else vals, Pretty.str s)
17.164 + end
17.165 + | Some ((_, s), _) => (vals, Pretty.str s));
17.166 +
17.167 +fun get_cases sg =
17.168 + Symtab.foldl (fn (tab, (k, {case_rewrites, ...})) => Symtab.update_new
17.169 + ((fst (dest_Const (head_of (fst (HOLogic.dest_eq (HOLogic.dest_Trueprop
17.170 + (prop_of (hd case_rewrites))))))), map my_mk_meta_eq case_rewrites), tab))
17.171 + (Symtab.empty, DatatypePackage.get_datatypes_sg sg);
17.172 +
17.173 +fun decomp_case th =
17.174 + let
17.175 + val (lhs, _) = Logic.dest_equals (prop_of th);
17.176 + val (f, ts) = strip_comb lhs;
17.177 + val (us, u) = split_last ts;
17.178 + val (Const (s, _), vs) = strip_comb u
17.179 + in (us, s, vs, u) end;
17.180 +
17.181 +fun rename vs t =
17.182 + let
17.183 + fun mk_subst ((vs, subs), Var ((s, i), T)) =
17.184 + let val s' = variant vs s
17.185 + in if s = s' then (vs, subs)
17.186 + else (s' :: vs, ((s, i), Var ((s', i), T)) :: subs)
17.187 + end;
17.188 + val (vs', subs) = foldl mk_subst ((vs, []), term_vars t)
17.189 + in (vs', subst_Vars subs t) end;
17.190 +
17.191 +fun is_instance sg t u = t = subst_TVars_Vartab
17.192 + (Type.typ_match (Sign.tsig_of sg) (Vartab.empty,
17.193 + (fastype_of u, fastype_of t))) u handle Type.TYPE_MATCH => false;
17.194 +
17.195 +(*
17.196 +fun lookup sg fs t = apsome snd (Library.find_first
17.197 + (is_instance sg t o fst) fs);
17.198 +*)
17.199 +
17.200 +fun lookup sg fs t = (case Library.find_first (is_instance sg t o fst) fs of
17.201 + None => (bla := (t ins !bla); None)
17.202 + | Some (_, x) => Some x);
17.203 +
17.204 +fun unint sg fs t = forall (is_none o lookup sg fs) (term_consts' t);
17.205 +
17.206 +fun mk_let s i xs ys =
17.207 + Pretty.blk (0, [Pretty.blk (i, separate Pretty.fbrk (Pretty.str s :: xs)),
17.208 + Pretty.fbrk,
17.209 + Pretty.blk (i, ([Pretty.str "in", Pretty.fbrk] @ ys)),
17.210 + Pretty.fbrk, Pretty.str "end"]);
17.211 +
17.212 +(*****************************************************************************)
17.213 +(* Generate bindings for simplifying term t *)
17.214 +(* mkeq: whether to generate reflexivity theorem for uninterpreted terms *)
17.215 +(* fs: interpreted functions *)
17.216 +(* ts: atomic terms *)
17.217 +(* vs: used identifiers *)
17.218 +(* vals: list of bindings of the form ((eq, ct), ps) where *)
17.219 +(* eq: name of equational theorem *)
17.220 +(* ct: name of simplified cterm *)
17.221 +(* ps: ML code for creating the above two items *)
17.222 +(*****************************************************************************)
17.223 +
17.224 +fun mk_simpl_code sg case_tab mkeq fs ts ty_bs thm_bs ((vs, vals), t) =
17.225 + (case assoc (vals, t) of
17.226 + Some ((eq, ct), ps) => (* binding already generated *)
17.227 + if mkeq andalso eq="" then
17.228 + let val eq' = variant vs "eq"
17.229 + in ((eq' :: vs, overwrite (vals,
17.230 + (t, ((eq', ct), ps @ [mk_refleq eq' ct])))), (eq', ct))
17.231 + end
17.232 + else ((vs, vals), (eq, ct))
17.233 + | None => (case assoc (ts, t) of
17.234 + Some v => (* atomic term *)
17.235 + let val xs = if not (null (term_tvars t)) andalso is_Const t then
17.236 + [mk_val [v] [inst_ty false ty_bs t v]] else []
17.237 + in
17.238 + if mkeq then
17.239 + let val eq = variant vs "eq"
17.240 + in ((eq :: vs, vals @
17.241 + [(t, ((eq, v), xs @ [mk_refleq eq v]))]), (eq, v))
17.242 + end
17.243 + else ((vs, if null xs then vals else vals @
17.244 + [(t, (("", v), xs))]), ("", v))
17.245 + end
17.246 + | None => (* complex term *)
17.247 + let val (f as Const (cname, _), us) = strip_comb t
17.248 + in case Symtab.lookup (case_tab, cname) of
17.249 + Some cases => (* case expression *)
17.250 + let
17.251 + val (us', u) = split_last us;
17.252 + val b = unint sg fs u;
17.253 + val ((vs1, vals1), (eq, ct)) =
17.254 + mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs ((vs, vals), u);
17.255 + val xs = variantlist (replicate (length us') "f", vs1);
17.256 + val (vals2, ps) = foldl_map
17.257 + (mk_cterm_code false ty_bs ts []) (vals1, us');
17.258 + val fvals = map (fn (x, p) => mk_val [x] [p]) (xs ~~ ps);
17.259 + val uT = fastype_of u;
17.260 + val (us'', _, _, u') = decomp_case (hd cases);
17.261 + val (vs2, ty_bs', ty_vals) = mk_type_bindings
17.262 + (mk_type_bindings ((vs1 @ xs, [], []),
17.263 + (hd xs, fastype_of (hd us''))), (ct, fastype_of u'));
17.264 + val insts1 = map mk_tyinst ty_bs';
17.265 + val i = length vals2;
17.266 +
17.267 + fun mk_case_code ((vs, vals), (f, (name, eqn))) =
17.268 + let
17.269 + val (fvs, cname, cvs, _) = decomp_case eqn;
17.270 + val Ts = binder_types (fastype_of f);
17.271 + val ys = variantlist (map (fst o fst o dest_Var) cvs, vs);
17.272 + val cvs' = map Var (map (rpair 0) ys ~~ Ts);
17.273 + val rs = cvs' ~~ cvs;
17.274 + val lhs = list_comb (Const (cname, Ts ---> uT), cvs');
17.275 + val rhs = foldl betapply (f, cvs');
17.276 + val (vs', tm_bs, tm_vals) = decomp_term_code false
17.277 + ((vs @ ys, [], []), (ct, lhs));
17.278 + val ((vs'', all_vals), (eq', ct')) = mk_simpl_code sg case_tab
17.279 + false fs (tm_bs @ ts) ty_bs thm_bs ((vs', vals), rhs);
17.280 + val (old_vals, eq_vals) = splitAt (i, all_vals);
17.281 + val vs''' = vs @ filter (fn v => exists
17.282 + (fn (_, ((v', _), _)) => v = v') old_vals) (vs'' \\ vs');
17.283 + val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
17.284 + inst_ty false ty_bs' t (the (assoc (thm_bs, t))), Pretty.str ",",
17.285 + Pretty.brk 1, Pretty.str (s ^ ")")]) ((fvs ~~ xs) @
17.286 + (map (fn (v, s) => (the (assoc (rs, v)), s)) tm_bs));
17.287 + val eq'' = if null insts1 andalso null insts2 then Pretty.str name
17.288 + else parens (eq' <> "") [Pretty.str
17.289 + (if null cvs then "Thm.instantiate" else "Drule.instantiate"),
17.290 + Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
17.291 + Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
17.292 + Pretty.str ")", Pretty.brk 1, Pretty.str name];
17.293 + val eq''' = if eq' = "" then eq'' else
17.294 + Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
17.295 + eq'', Pretty.brk 1, Pretty.str eq']
17.296 + in
17.297 + ((vs''', old_vals), Pretty.block [pretty_pattern false lhs,
17.298 + Pretty.str " =>",
17.299 + Pretty.brk 1, mk_let "let" 2 (tm_vals @ flat (map (snd o snd) eq_vals))
17.300 + [Pretty.str ("(" ^ ct' ^ ","), Pretty.brk 1, eq''', Pretty.str ")"]])
17.301 + end;
17.302 +
17.303 + val case_names = map (fn i => Sign.base_name cname ^ "_" ^
17.304 + string_of_int i) (1 upto length cases);
17.305 + val ((vs3, vals3), case_ps) = foldl_map mk_case_code
17.306 + ((vs2, vals2), us' ~~ (case_names ~~ cases));
17.307 + val eq' = variant vs3 "eq";
17.308 + val ct' = variant (eq' :: vs3) "ct";
17.309 + val eq'' = variant (eq' :: ct' :: vs3) "eq";
17.310 + val case_vals =
17.311 + fvals @ ty_vals @
17.312 + [mk_val [ct', eq'] ([Pretty.str "(case", Pretty.brk 1,
17.313 + Pretty.str ("term_of " ^ ct ^ " of"), Pretty.brk 1] @
17.314 + flat (separate [Pretty.brk 1, Pretty.str "| "]
17.315 + (map single case_ps)) @ [Pretty.str ")"])]
17.316 + in
17.317 + if b then
17.318 + ((eq' :: ct' :: vs3, vals3 @
17.319 + [(t, ((eq', ct'), case_vals))]), (eq', ct'))
17.320 + else
17.321 + let val ((vs4, vals4), (_, ctcase)) = mk_simpl_code sg case_tab false
17.322 + fs ts ty_bs thm_bs ((eq' :: eq'' :: ct' :: vs3, vals3), f)
17.323 + in
17.324 + ((vs4, vals4 @ [(t, ((eq'', ct'), case_vals @
17.325 + [mk_val [eq''] [Pretty.str "Thm.transitive", Pretty.brk 1,
17.326 + Pretty.str "(Thm.combination", Pretty.brk 1,
17.327 + Pretty.str "(Thm.reflexive", Pretty.brk 1,
17.328 + mk_apps "Thm.capply" true (Pretty.str ctcase)
17.329 + (map Pretty.str xs),
17.330 + Pretty.str ")", Pretty.brk 1, Pretty.str (eq ^ ")"),
17.331 + Pretty.brk 1, Pretty.str eq']]))]), (eq'', ct'))
17.332 + end
17.333 + end
17.334 +
17.335 + | None =>
17.336 + let
17.337 + val b = forall (unint sg fs) us;
17.338 + val (q, eqs) = foldl_map
17.339 + (mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs) ((vs, vals), us);
17.340 + val ((vs', vals'), (eqf, ctf)) = if is_some (lookup sg fs f) andalso b
17.341 + then (q, ("", ""))
17.342 + else mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs (q, f);
17.343 + val ct = variant vs' "ct";
17.344 + val eq = variant (ct :: vs') "eq";
17.345 + val ctv = mk_val [ct] [mk_apps "Thm.capply" false
17.346 + (Pretty.str ctf) (map (Pretty.str o snd) eqs)];
17.347 + fun combp b = mk_apps "Thm.combination" b
17.348 + (Pretty.str eqf) (map (Pretty.str o fst) eqs)
17.349 + in
17.350 + case (lookup sg fs f, b) of
17.351 + (None, true) => (* completely uninterpreted *)
17.352 + if mkeq then ((ct :: eq :: vs', vals' @
17.353 + [(t, ((eq, ct), [ctv, mk_refleq eq ct]))]), (eq, ct))
17.354 + else ((ct :: vs', vals' @ [(t, (("", ct), [ctv]))]), ("", ct))
17.355 + | (None, false) => (* function uninterpreted *)
17.356 + ((eq :: ct :: vs', vals' @
17.357 + [(t, ((eq, ct), [ctv, mk_val [eq] [combp false]]))]), (eq, ct))
17.358 + | (Some (s, _, _), true) => (* arguments uninterpreted *)
17.359 + ((eq :: ct :: vs', vals' @
17.360 + [(t, ((eq, ct), [mk_val [ct, eq] (separate (Pretty.brk 1)
17.361 + (Pretty.str s :: map (Pretty.str o snd) eqs))]))]), (eq, ct))
17.362 + | (Some (s, _, _), false) => (* function and arguments interpreted *)
17.363 + let val eq' = variant (eq :: ct :: vs') "eq"
17.364 + in ((eq' :: eq :: ct :: vs', vals' @ [(t, ((eq', ct),
17.365 + [mk_val [ct, eq] (separate (Pretty.brk 1)
17.366 + (Pretty.str s :: map (Pretty.str o snd) eqs)),
17.367 + mk_val [eq'] [Pretty.str "Thm.transitive", Pretty.brk 1,
17.368 + combp true, Pretty.brk 1, Pretty.str eq]]))]), (eq', ct))
17.369 + end
17.370 + end
17.371 + end));
17.372 +
17.373 +fun lhs_of thm = fst (Logic.dest_equals (prop_of thm));
17.374 +fun rhs_of thm = snd (Logic.dest_equals (prop_of thm));
17.375 +
17.376 +fun mk_funs_code sg case_tab fs fs' =
17.377 + let
17.378 + val case_thms = mapfilter (fn s => (case Symtab.lookup (case_tab, s) of
17.379 + None => None
17.380 + | Some thms => Some (unsuffix "_case" (Sign.base_name s) ^ ".cases",
17.381 + map (fn i => Sign.base_name s ^ "_" ^ string_of_int i)
17.382 + (1 upto length thms) ~~ thms)))
17.383 + (foldr add_term_consts (map (prop_of o snd)
17.384 + (flat (map (#3 o snd) fs')), []));
17.385 + val case_vals = map (fn (s, cs) => mk_vall (map fst cs)
17.386 + [Pretty.str "map my_mk_meta_eq", Pretty.brk 1,
17.387 + Pretty.str ("(thms \"" ^ s ^ "\")")]) case_thms;
17.388 + val (vs, thm_bs, thm_vals) = foldl mk_term_bindings (([], [], []),
17.389 + flat (map (map (apsnd prop_of) o #3 o snd) fs') @
17.390 + map (apsnd prop_of) (flat (map snd case_thms)));
17.391 +
17.392 + fun mk_fun_code (prfx, (fname, d, eqns)) =
17.393 + let
17.394 + val (f, ts) = strip_comb (lhs_of (snd (hd eqns)));
17.395 + val args = variantlist (replicate (length ts) "ct", vs);
17.396 + val (vs', ty_bs, ty_vals) = foldl mk_type_bindings
17.397 + ((vs @ args, [], []), args ~~ map fastype_of ts);
17.398 + val insts1 = map mk_tyinst ty_bs;
17.399 +
17.400 + fun mk_eqn_code (name, eqn) =
17.401 + let
17.402 + val (_, argts) = strip_comb (lhs_of eqn);
17.403 + val (vs'', tm_bs, tm_vals) = foldl (decomp_term_code false)
17.404 + ((vs', [], []), args ~~ argts);
17.405 + val ((vs''', eq_vals), (eq, ct)) = mk_simpl_code sg case_tab false fs
17.406 + (tm_bs @ filter_out (is_Var o fst) thm_bs) ty_bs thm_bs
17.407 + ((vs'', []), rhs_of eqn);
17.408 + val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(",
17.409 + inst_ty false ty_bs t (the (assoc (thm_bs, t))), Pretty.str ",", Pretty.brk 1,
17.410 + Pretty.str (s ^ ")")]) tm_bs
17.411 + val eq' = if null insts1 andalso null insts2 then Pretty.str name
17.412 + else parens (eq <> "") [Pretty.str "Thm.instantiate",
17.413 + Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1,
17.414 + Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2,
17.415 + Pretty.str ")", Pretty.brk 1, Pretty.str name];
17.416 + val eq'' = if eq = "" then eq' else
17.417 + Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1,
17.418 + eq', Pretty.brk 1, Pretty.str eq]
17.419 + in
17.420 + Pretty.block [parens (length argts > 1)
17.421 + (Pretty.commas (map (pretty_pattern false) argts)),
17.422 + Pretty.str " =>",
17.423 + Pretty.brk 1, mk_let "let" 2 (ty_vals @ tm_vals @ flat (map (snd o snd) eq_vals))
17.424 + [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1, eq'', Pretty.str ")"]]
17.425 + end;
17.426 +
17.427 + val default = if d then
17.428 + let
17.429 + val Some s = assoc (thm_bs, f);
17.430 + val ct = variant vs' "ct"
17.431 + in [Pretty.brk 1, Pretty.str "handle", Pretty.brk 1,
17.432 + Pretty.str "Match =>", Pretty.brk 1, mk_let "let" 2
17.433 + (ty_vals @ (if null (term_tvars f) then [] else
17.434 + [mk_val [s] [inst_ty false ty_bs f s]]) @
17.435 + [mk_val [ct] [mk_apps "Thm.capply" false (Pretty.str s)
17.436 + (map Pretty.str args)]])
17.437 + [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1,
17.438 + Pretty.str "Thm.reflexive", Pretty.brk 1, Pretty.str (ct ^ ")")]]
17.439 + end
17.440 + else []
17.441 + in
17.442 + ("and ", Pretty.block (separate (Pretty.brk 1)
17.443 + (Pretty.str (prfx ^ fname) :: map Pretty.str args) @
17.444 + [Pretty.str " =", Pretty.brk 1, Pretty.str "(case", Pretty.brk 1,
17.445 + Pretty.list "(" ")" (map (fn s => Pretty.str ("term_of " ^ s)) args),
17.446 + Pretty.str " of", Pretty.brk 1] @
17.447 + flat (separate [Pretty.brk 1, Pretty.str "| "]
17.448 + (map (single o mk_eqn_code) eqns)) @ [Pretty.str ")"] @ default))
17.449 + end;
17.450 +
17.451 + val (_, decls) = foldl_map mk_fun_code ("fun ", map snd fs')
17.452 + in
17.453 + mk_let "local" 2 (case_vals @ thm_vals) (separate Pretty.fbrk decls)
17.454 + end;
17.455 +
17.456 +fun mk_simprocs_code sg eqns =
17.457 + let
17.458 + val case_tab = get_cases sg;
17.459 + fun get_head th = head_of (fst (Logic.dest_equals (prop_of th)));
17.460 + fun attach_term (x as (_, _, (_, th) :: _)) = (get_head th, x);
17.461 + val eqns' = map attach_term eqns;
17.462 + fun mk_node (s, _, (_, th) :: _) = (s, get_head th);
17.463 + fun mk_edges (s, _, ths) = map (pair s) (distinct
17.464 + (mapfilter (fn t => apsome #1 (lookup sg eqns' t))
17.465 + (flat (map (term_consts' o prop_of o snd) ths))));
17.466 + val gr = foldr (uncurry Graph.add_edge)
17.467 + (map (pair "" o #1) eqns @ flat (map mk_edges eqns),
17.468 + foldr (uncurry Graph.new_node)
17.469 + (("", Bound 0) :: map mk_node eqns, Graph.empty));
17.470 + val keys = rev (Graph.all_succs gr [""] \ "");
17.471 + fun gr_ord (x :: _, y :: _) =
17.472 + int_ord (find_index (equal x) keys, find_index (equal y) keys);
17.473 + val scc = map (fn xs => filter (fn (_, (s, _, _)) => s mem xs) eqns')
17.474 + (sort gr_ord (Graph.strong_conn gr \ [""]));
17.475 + in
17.476 + flat (separate [Pretty.str ";", Pretty.fbrk, Pretty.str " ", Pretty.fbrk]
17.477 + (map (fn eqns'' => [mk_funs_code sg case_tab eqns' eqns'']) scc)) @
17.478 + [Pretty.str ";", Pretty.fbrk]
17.479 + end;
17.480 +
17.481 +fun use_simprocs_code sg eqns =
17.482 + let
17.483 + fun attach_name (i, x) = (i+1, ("simp_thm_" ^ string_of_int i, x));
17.484 + fun attach_names (i, (s, b, eqs)) =
17.485 + let val (i', eqs') = foldl_map attach_name (i, eqs)
17.486 + in (i', (s, b, eqs')) end;
17.487 + val (_, eqns') = foldl_map attach_names (1, eqns);
17.488 + val (names, thms) = split_list (flat (map #3 eqns'));
17.489 + val s = setmp print_mode [] Pretty.string_of
17.490 + (mk_let "local" 2 [mk_vall names [Pretty.str "!SimprocsCodegen.simp_thms"]]
17.491 + (mk_simprocs_code sg eqns'))
17.492 + in
17.493 + (simp_thms := thms; use_text Context.ml_output false s)
17.494 + end;
17.495 +
17.496 +end;
18.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
18.2 +++ b/src/HOL/Matrix/fspmlp.ML Fri Sep 03 17:10:36 2004 +0200
18.3 @@ -0,0 +1,303 @@
18.4 +signature FSPMLP =
18.5 +sig
18.6 + type linprog
18.7 +
18.8 + val y : linprog -> cterm
18.9 + val A : linprog -> cterm * cterm
18.10 + val b : linprog -> cterm
18.11 + val c : linprog -> cterm * cterm
18.12 + val r : linprog -> cterm
18.13 +
18.14 + exception Load of string
18.15 +
18.16 + val load : string -> int -> bool -> linprog
18.17 +end
18.18 +
18.19 +structure fspmlp : FSPMLP =
18.20 +struct
18.21 +
18.22 +type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm
18.23 +
18.24 +fun y (c1, c2, c3, c4, c5) = c1
18.25 +fun A (c1, c2, c3, c4, c5) = c2
18.26 +fun b (c1, c2, c3, c4, c5) = c3
18.27 +fun c (c1, c2, c3, c4, c5) = c4
18.28 +fun r (c1, c2, c3, c4, c5) = c5
18.29 +
18.30 +structure CplexFloatSparseMatrixConverter =
18.31 +MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder);
18.32 +
18.33 +datatype bound_type = LOWER | UPPER
18.34 +
18.35 +fun intbound_ord ((i1, b1),(i2,b2)) =
18.36 + if i1 < i2 then LESS
18.37 + else if i1 = i2 then
18.38 + (if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER)
18.39 + else GREATER
18.40 +
18.41 +structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord));
18.42 +
18.43 +structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord);
18.44 +(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *)
18.45 +(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *)
18.46 +
18.47 +exception Internal of string;
18.48 +
18.49 +fun add_row_bound g dest_key row_index row_bound =
18.50 + let
18.51 + val x =
18.52 + case VarGraph.lookup (g, dest_key) of
18.53 + None => (None, Inttab.update ((row_index, (row_bound, [])), Inttab.empty))
18.54 + | Some (sure_bound, f) =>
18.55 + (sure_bound,
18.56 + case Inttab.lookup (f, row_index) of
18.57 + None => Inttab.update ((row_index, (row_bound, [])), f)
18.58 + | Some _ => raise (Internal "add_row_bound"))
18.59 + in
18.60 + VarGraph.update ((dest_key, x), g)
18.61 + end
18.62 +
18.63 +fun update_sure_bound g (key as (_, btype)) bound =
18.64 + let
18.65 + val x =
18.66 + case VarGraph.lookup (g, key) of
18.67 + None => (Some bound, Inttab.empty)
18.68 + | Some (None, f) => (Some bound, f)
18.69 + | Some (Some old_bound, f) =>
18.70 + (Some ((case btype of
18.71 + UPPER => FloatArith.min
18.72 + | LOWER => FloatArith.max)
18.73 + old_bound bound), f)
18.74 + in
18.75 + VarGraph.update ((key, x), g)
18.76 + end
18.77 +
18.78 +fun get_sure_bound g key =
18.79 + case VarGraph.lookup (g, key) of
18.80 + None => None
18.81 + | Some (sure_bound, _) => sure_bound
18.82 +
18.83 +(*fun get_row_bound g key row_index =
18.84 + case VarGraph.lookup (g, key) of
18.85 + None => None
18.86 + | Some (sure_bound, f) =>
18.87 + (case Inttab.lookup (f, row_index) of
18.88 + None => None
18.89 + | Some (row_bound, _) => (sure_bound, row_bound))*)
18.90 +
18.91 +fun add_edge g src_key dest_key row_index coeff =
18.92 + case VarGraph.lookup (g, dest_key) of
18.93 + None => raise (Internal "add_edge: dest_key not found")
18.94 + | Some (sure_bound, f) =>
18.95 + (case Inttab.lookup (f, row_index) of
18.96 + None => raise (Internal "add_edge: row_index not found")
18.97 + | Some (row_bound, sources) =>
18.98 + VarGraph.update ((dest_key, (sure_bound, Inttab.update ((row_index, (row_bound, (coeff, src_key) :: sources)), f))), g))
18.99 +
18.100 +fun split_graph g =
18.101 + let
18.102 + fun split ((r1, r2), (key, (sure_bound, _))) =
18.103 + case sure_bound of
18.104 + None => (r1, r2)
18.105 + | Some bound =>
18.106 + (case key of
18.107 + (u, UPPER) => (r1, Inttab.update ((u, bound), r2))
18.108 + | (u, LOWER) => (Inttab.update ((u, bound), r1), r2))
18.109 + in
18.110 + VarGraph.foldl split ((Inttab.empty, Inttab.empty), g)
18.111 + end
18.112 +
18.113 +fun it2list t =
18.114 + let
18.115 + fun tolist (l, a) = a::l
18.116 + in
18.117 + Inttab.foldl tolist ([], t)
18.118 + end
18.119 +
18.120 +(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm).
18.121 + If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *)
18.122 +fun propagate_sure_bounds safe names g =
18.123 + let
18.124 + (* returns None if no new sure bound could be calculated, otherwise the new sure bound is returned *)
18.125 + fun calc_sure_bound_from_sources g (key as (_, btype)) =
18.126 + let
18.127 + fun mult_upper x (lower, upper) =
18.128 + if FloatArith.is_negative x then
18.129 + FloatArith.mul x lower
18.130 + else
18.131 + FloatArith.mul x upper
18.132 +
18.133 + fun mult_lower x (lower, upper) =
18.134 + if FloatArith.is_negative x then
18.135 + FloatArith.mul x upper
18.136 + else
18.137 + FloatArith.mul x lower
18.138 +
18.139 + val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower
18.140 +
18.141 + fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) =
18.142 + let
18.143 + fun add_src_bound (sum, (coeff, src_key)) =
18.144 + case sum of
18.145 + None => None
18.146 + | Some x =>
18.147 + (case get_sure_bound g src_key of
18.148 + None => None
18.149 + | Some src_sure_bound => Some (FloatArith.add x (mult_btype src_sure_bound coeff)))
18.150 + in
18.151 + case foldl add_src_bound (Some row_bound, sources) of
18.152 + None => sure_bound
18.153 + | new_sure_bound as (Some new_bound) =>
18.154 + (case sure_bound of
18.155 + None => new_sure_bound
18.156 + | Some old_bound =>
18.157 + Some (case btype of
18.158 + UPPER => FloatArith.min old_bound new_bound
18.159 + | LOWER => FloatArith.max old_bound new_bound))
18.160 + end
18.161 + in
18.162 + case VarGraph.lookup (g, key) of
18.163 + None => None
18.164 + | Some (sure_bound, f) =>
18.165 + let
18.166 + val x = Inttab.foldl calc_sure_bound (sure_bound, f)
18.167 + in
18.168 + if x = sure_bound then None else x
18.169 + end
18.170 + end
18.171 +
18.172 + fun propagate ((g, b), (key, _)) =
18.173 + case calc_sure_bound_from_sources g key of
18.174 + None => (g,b)
18.175 + | Some bound => (update_sure_bound g key bound,
18.176 + if safe then
18.177 + case get_sure_bound g key of
18.178 + None => true
18.179 + | _ => b
18.180 + else
18.181 + true)
18.182 +
18.183 + val (g, b) = VarGraph.foldl propagate ((g, false), g)
18.184 + in
18.185 + if b then propagate_sure_bounds safe names g else g
18.186 + end
18.187 +
18.188 +exception Load of string;
18.189 +
18.190 +fun calcr safe_propagation xlen names prec A b =
18.191 + let
18.192 + val empty = Inttab.empty
18.193 +
18.194 + fun instab t i x = Inttab.update ((i,x), t)
18.195 +
18.196 + fun isnegstr x = String.isPrefix "-" x
18.197 + fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x
18.198 +
18.199 + fun test_1 (lower, upper) =
18.200 + if lower = upper then
18.201 + (if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1
18.202 + else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1
18.203 + else 0)
18.204 + else 0
18.205 +
18.206 + fun calcr (g, (row_index, a)) =
18.207 + let
18.208 + val b = FloatSparseMatrixBuilder.v_elem_at b row_index
18.209 + val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of None => "0" | Some b => b)
18.210 + val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) =>
18.211 + (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a
18.212 +
18.213 + fun fold_dest_nodes (g, (dest_index, dest_value)) =
18.214 + let
18.215 + val dest_test = test_1 dest_value
18.216 + in
18.217 + if dest_test = 0 then
18.218 + g
18.219 + else let
18.220 + val (dest_key as (_, dest_btype), row_bound) =
18.221 + if dest_test = ~1 then
18.222 + ((dest_index, LOWER), FloatArith.neg b2)
18.223 + else
18.224 + ((dest_index, UPPER), b2)
18.225 +
18.226 + fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) =
18.227 + if src_index = dest_index then g
18.228 + else
18.229 + let
18.230 + val coeff = case dest_btype of
18.231 + UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower)
18.232 + | LOWER => src_value
18.233 + in
18.234 + if FloatArith.is_negative src_lower then
18.235 + add_edge g (src_index, UPPER) dest_key row_index coeff
18.236 + else
18.237 + add_edge g (src_index, LOWER) dest_key row_index coeff
18.238 + end
18.239 + in
18.240 + foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a)
18.241 + end
18.242 + end
18.243 + in
18.244 + case approx_a of
18.245 + [] => g
18.246 + | [(u, a)] =>
18.247 + let
18.248 + val atest = test_1 a
18.249 + in
18.250 + if atest = ~1 then
18.251 + update_sure_bound g (u, LOWER) (FloatArith.neg b2)
18.252 + else if atest = 1 then
18.253 + update_sure_bound g (u, UPPER) b2
18.254 + else
18.255 + g
18.256 + end
18.257 + | _ => foldl fold_dest_nodes (g, approx_a)
18.258 + end
18.259 +
18.260 + val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A
18.261 + val g = propagate_sure_bounds safe_propagation names g
18.262 +
18.263 + val (r1, r2) = split_graph g
18.264 +
18.265 + fun abs_estimate i r1 r2 =
18.266 + if i = 0 then FloatSparseMatrixBuilder.empty_spmat
18.267 + else
18.268 + let
18.269 + val index = xlen-i
18.270 + val r = abs_estimate (i-1) r1 r2
18.271 + val b1 = case Inttab.lookup (r1, index) of None => raise (Load ("x-value not bounded from below: "^(names index))) | Some x => x
18.272 + val b2 = case Inttab.lookup (r2, index) of None => raise (Load ("x-value not bounded from above: "^(names index))) | Some x => x
18.273 + val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2)
18.274 + val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec
18.275 + in
18.276 + FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r
18.277 + end
18.278 + in
18.279 + FloatSparseMatrixBuilder.sign_term (abs_estimate xlen r1 r2)
18.280 + end
18.281 +
18.282 +fun load filename prec safe_propagation =
18.283 + let
18.284 + val prog = Cplex.load_cplexFile filename
18.285 + val prog = Cplex.elim_nonfree_bounds prog
18.286 + val prog = Cplex.relax_strict_ineqs prog
18.287 + val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog
18.288 + val r = calcr safe_propagation xlen names prec A b
18.289 + val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems"
18.290 + val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b
18.291 + val results = Cplex.solve dualprog
18.292 + val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof
18.293 + val A = FloatSparseMatrixBuilder.cut_matrix v None A
18.294 + fun id x = x
18.295 + val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v
18.296 + val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b)
18.297 + val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c
18.298 + val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v
18.299 + val A = FloatSparseMatrixBuilder.approx_matrix prec id A
18.300 + val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b
18.301 + val c = FloatSparseMatrixBuilder.approx_matrix prec id c
18.302 + in
18.303 + (y1, A, b2, c, r)
18.304 + end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s)))
18.305 +
18.306 +end
18.307 \ No newline at end of file
19.1 --- a/src/HOL/OrderedGroup.thy Fri Sep 03 10:27:05 2004 +0200
19.2 +++ b/src/HOL/OrderedGroup.thy Fri Sep 03 17:10:36 2004 +0200
19.3 @@ -842,6 +842,42 @@
19.4 lemma minus_add_cancel: "-(a::'a::ab_group_add) + (a + b) = b"
19.5 by (simp add: add_assoc[symmetric])
19.6
19.7 +lemma le_add_right_mono:
19.8 + assumes
19.9 + "a <= b + (c::'a::pordered_ab_group_add)"
19.10 + "c <= d"
19.11 + shows "a <= b + d"
19.12 + apply (rule_tac order_trans[where y = "b+c"])
19.13 + apply (simp_all add: prems)
19.14 + done
19.15 +
19.16 +lemmas group_eq_simps =
19.17 + mult_ac
19.18 + add_ac
19.19 + add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
19.20 + diff_eq_eq eq_diff_eq
19.21 +
19.22 +lemma estimate_by_abs:
19.23 +"a + b <= (c::'a::lordered_ab_group_abs) \<Longrightarrow> a <= c + abs b"
19.24 +proof -
19.25 + assume 1: "a+b <= c"
19.26 + have 2: "a <= c+(-b)"
19.27 + apply (insert 1)
19.28 + apply (drule_tac add_right_mono[where c="-b"])
19.29 + apply (simp add: group_eq_simps)
19.30 + done
19.31 + have 3: "(-b) <= abs b" by (rule abs_ge_minus_self)
19.32 + show ?thesis by (rule le_add_right_mono[OF 2 3])
19.33 +qed
19.34 +
19.35 +lemma abs_of_ge_0: "0 <= (y::'a::lordered_ab_group_abs) \<Longrightarrow> abs y = y"
19.36 +proof -
19.37 + assume 1:"0 <= y"
19.38 + have 2:"-y <= 0" by (simp add: 1)
19.39 + from 1 2 have 3:"-y <= y" by (simp only:)
19.40 + show ?thesis by (simp add: abs_lattice ge_imp_join_eq[OF 3])
19.41 +qed
19.42 +
19.43 ML {*
19.44 val add_zero_left = thm"add_0";
19.45 val add_zero_right = thm"add_0_right";
20.1 --- a/src/HOL/Ring_and_Field.thy Fri Sep 03 10:27:05 2004 +0200
20.2 +++ b/src/HOL/Ring_and_Field.thy Fri Sep 03 17:10:36 2004 +0200
20.3 @@ -546,14 +546,14 @@
20.4 done
20.5
20.6 text{*This list of rewrites decides ring equalities by ordered rewriting.*}
20.7 -lemmas ring_eq_simps =
20.8 - mult_ac
20.9 +lemmas ring_eq_simps =
20.10 +(* mult_ac*)
20.11 left_distrib right_distrib left_diff_distrib right_diff_distrib
20.12 - add_ac
20.13 + group_eq_simps
20.14 +(* add_ac
20.15 add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2
20.16 - diff_eq_eq eq_diff_eq
20.17 + diff_eq_eq eq_diff_eq *)
20.18
20.19 -
20.20 subsection {* Fields *}
20.21
20.22 lemma right_inverse [simp]:
20.23 @@ -1503,6 +1503,90 @@
20.24 text{*Moving this up spoils many proofs using @{text mult_le_cancel_right}*}
20.25 declare times_divide_eq_left [simp]
20.26
20.27 +lemma linprog_dual_estimate:
20.28 + assumes
20.29 + "A * x \<le> (b::'a::lordered_ring)"
20.30 + "0 \<le> y"
20.31 + "abs (A - A') \<le> \<delta>A"
20.32 + "b \<le> b'"
20.33 + "abs (c - c') \<le> \<delta>c"
20.34 + "abs x \<le> r"
20.35 + shows
20.36 + "c * x \<le> y * b' + (y * \<delta>A + abs (y * A' - c') + \<delta>c) * r"
20.37 +proof -
20.38 + from prems have 1: "y * b <= y * b'" by (simp add: mult_left_mono)
20.39 + from prems have 2: "y * (A * x) <= y * b" by (simp add: mult_left_mono)
20.40 + have 3: "y * (A * x) = c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x" by (simp add: ring_eq_simps)
20.41 + from 1 2 3 have 4: "c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x <= y * b'" by simp
20.42 + have 5: "c * x <= y * b' + abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"
20.43 + by (simp only: 4 estimate_by_abs)
20.44 + have 6: "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <= abs (y * (A - A') + (y * A' - c') + (c'-c)) * abs x"
20.45 + by (simp add: abs_le_mult)
20.46 + have 7: "(abs (y * (A - A') + (y * A' - c') + (c'-c))) * abs x <= (abs (y * (A-A') + (y*A'-c')) + abs(c'-c)) * abs x"
20.47 + by (simp add: abs_triangle_ineq mult_right_mono)
20.48 + have 8: " (abs (y * (A-A') + (y*A'-c')) + abs(c'-c)) * abs x <= (abs (y * (A-A')) + abs (y*A'-c') + abs(c'-c)) * abs x"
20.49 + by (simp add: abs_triangle_ineq mult_right_mono)
20.50 + have 9: "(abs (y * (A-A')) + abs (y*A'-c') + abs(c'-c)) * abs x <= (abs y * abs (A-A') + abs (y*A'-c') + abs (c'-c)) * abs x"
20.51 + by (simp add: abs_le_mult mult_right_mono)
20.52 + have 10: "c'-c = -(c-c')" by (simp add: ring_eq_simps)
20.53 + have 11: "abs (c'-c) = abs (c-c')"
20.54 + by (subst 10, subst abs_minus_cancel, simp)
20.55 + have 12: "(abs y * abs (A-A') + abs (y*A'-c') + abs (c'-c)) * abs x <= (abs y * abs (A-A') + abs (y*A'-c') + \<delta>c) * abs x"
20.56 + by (simp add: 11 prems mult_right_mono)
20.57 + have 13: "(abs y * abs (A-A') + abs (y*A'-c') + \<delta>c) * abs x <= (abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * abs x"
20.58 + by (simp add: prems mult_right_mono mult_left_mono)
20.59 + have r: "(abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * abs x <= (abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * r"
20.60 + apply (rule mult_left_mono)
20.61 + apply (simp add: prems)
20.62 + apply (rule_tac add_mono[of "0::'a" _ "0", simplified])+
20.63 + apply (rule mult_left_mono[of "0" "\<delta>A", simplified])
20.64 + apply (simp_all)
20.65 + apply (rule order_trans[where y="abs (A-A')"], simp_all add: prems)
20.66 + apply (rule order_trans[where y="abs (c-c')"], simp_all add: prems)
20.67 + done
20.68 + from 6 7 8 9 12 13 r have 14:" abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <=(abs y * \<delta>A + abs (y*A'-c') + \<delta>c) * r"
20.69 + by (simp)
20.70 + show ?thesis
20.71 + apply (rule_tac le_add_right_mono[of _ _ "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"])
20.72 + apply (simp_all add: 5 14[simplified abs_of_ge_0[of y, simplified prems]])
20.73 + done
20.74 +qed
20.75 +
20.76 +lemma le_ge_imp_abs_diff_1:
20.77 + assumes
20.78 + "A1 <= (A::'a::lordered_ring)"
20.79 + "A <= A2"
20.80 + shows "abs (A-A1) <= A2-A1"
20.81 +proof -
20.82 + have "0 <= A - A1"
20.83 + proof -
20.84 + have 1: "A - A1 = A + (- A1)" by simp
20.85 + show ?thesis by (simp only: 1 add_right_mono[of A1 A "-A1", simplified, simplified prems])
20.86 + qed
20.87 + then have "abs (A-A1) = A-A1" by (rule abs_of_ge_0)
20.88 + with prems show "abs (A-A1) <= (A2-A1)" by simp
20.89 +qed
20.90 +
20.91 +lemma linprog_dual_estimate_1:
20.92 + assumes
20.93 + "A * x \<le> (b::'a::lordered_ring)"
20.94 + "0 \<le> y"
20.95 + "A1 <= A"
20.96 + "A <= A2"
20.97 + "c1 <= c"
20.98 + "c <= c2"
20.99 + "abs x \<le> r"
20.100 + shows
20.101 + "c * x \<le> y * b + (y * (A2 - A1) + abs (y * A1 - c1) + (c2 - c1)) * r"
20.102 +proof -
20.103 + from prems have delta_A: "abs (A-A1) <= (A2-A1)" by (simp add: le_ge_imp_abs_diff_1)
20.104 + from prems have delta_c: "abs (c-c1) <= (c2-c1)" by (simp add: le_ge_imp_abs_diff_1)
20.105 + show ?thesis
20.106 + apply (rule_tac linprog_dual_estimate)
20.107 + apply (auto intro: delta_A delta_c simp add: prems)
20.108 + done
20.109 +qed
20.110 +
20.111 ML {*
20.112 val left_distrib = thm "left_distrib";
20.113 val right_distrib = thm "right_distrib";