# HG changeset patch # User obua # Date 1094224236 -7200 # Node ID 5f621aa35c259af3876d0f3a1448f7303b25c131 # Parent e7616269fdca84592c31a8e81b12b6a916b8abc4 Matrix theory, linear programming diff -r e7616269fdca -r 5f621aa35c25 src/HOL/IsaMakefile --- a/src/HOL/IsaMakefile Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/IsaMakefile Fri Sep 03 17:10:36 2004 +0200 @@ -1,5 +1,3 @@ -# -# $Id$ # # IsaMakefile for HOL # @@ -80,8 +78,7 @@ $(SRC)/Provers/quasi.ML \ $(SRC)/Provers/simplifier.ML $(SRC)/Provers/splitter.ML \ $(SRC)/Provers/trancl.ML \ - $(SRC)/TFL/isand.ML $(SRC)/TFL/casesplit.ML $(SRC)/TFL/dcterm.ML\ - $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \ + $(SRC)/TFL/dcterm.ML $(SRC)/TFL/post.ML $(SRC)/TFL/rules.ML \ $(SRC)/TFL/tfl.ML $(SRC)/TFL/thms.ML $(SRC)/TFL/thry.ML \ $(SRC)/TFL/usyntax.ML $(SRC)/TFL/utils.ML \ Datatype.thy Datatype_Universe.ML Datatype_Universe.thy \ @@ -97,8 +94,7 @@ Nat.thy NatArith.ML NatArith.thy \ Power.thy PreList.thy Product_Type.ML Product_Type.thy \ Refute.thy ROOT.ML \ - Recdef.thy Reconstruction.thy Record.thy\ - Relation.ML Relation.thy Relation_Power.ML \ + Recdef.thy Record.thy Relation.ML Relation.thy Relation_Power.ML \ Relation_Power.thy LOrder.thy OrderedGroup.thy OrderedGroup.ML Ring_and_Field.thy\ Set.ML Set.thy SetInterval.ML SetInterval.thy \ Sum_Type.ML Sum_Type.thy Tools/datatype_abs_proofs.ML Tools/datatype_aux.ML \ @@ -109,7 +105,7 @@ Tools/primrec_package.ML \ Tools/prop_logic.ML \ Tools/recdef_package.ML Tools/recfun_codegen.ML \ - Tools/record_package.ML Tools/reconstruction.ML\ + Tools/record_package.ML \ Tools/refute.ML Tools/refute_isar.ML \ Tools/rewrite_hol_proof.ML \ Tools/sat_solver.ML \ @@ -635,14 +631,18 @@ ## HOL-Matrix -HOL-Matrix: HOL $(LOG)/HOL-Matrix.gz +HOL-Matrix: HOL HOL-Complex $(LOG)/HOL-Matrix.gz -$(LOG)/HOL-Matrix.gz: $(OUT)/HOL Matrix/ROOT.ML \ - Matrix/Matrix.thy Matrix/LinProg.thy Matrix/MatrixGeneral.thy \ - Matrix/SparseMatrix.thy Matrix/document/root.tex - @$(ISATOOL) usedir $(OUT)/HOL Matrix - - +$(LOG)/HOL-Matrix.gz: $(OUT)/HOL-Complex \ + Matrix/MatrixGeneral.thy Matrix/Matrix.thy Matrix/SparseMatrix.thy \ + Matrix/Float.thy Matrix/FloatArith.ML Matrix/ExactFloatingPoint.ML \ + Matrix/Cplex.ML Matrix/CplexMatrixConverter.ML \ + Matrix/FloatSparseMatrixBuilder.ML \ + Matrix/conv.ML Matrix/eq_codegen.ML Matrix/codegen_prep.ML \ + Matrix/fspmlp.ML \ + Matrix/MatrixLP_gensimp.ML Matrix/MatrixLP.ML Matrix/MatrixLP.thy \ + Matrix/document/root.tex Matrix/ROOT.ML + @cd Matrix; $(ISATOOL) usedir -b $(OUT)/HOL-Complex HOL-Matrix ## TLA diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/Cplex.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/Cplex.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,1027 @@ +signature CPLEX = +sig + + datatype cplexTerm = cplexVar of string | cplexNum of string | cplexInf + | cplexNeg of cplexTerm + | cplexProd of cplexTerm * cplexTerm + | cplexSum of (cplexTerm list) + + datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq + + datatype cplexGoal = cplexMinimize of cplexTerm + | cplexMaximize of cplexTerm + + datatype cplexConstr = cplexConstr of cplexComp * + (cplexTerm * cplexTerm) + + datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm + * cplexComp * cplexTerm + | cplexBound of cplexTerm * cplexComp * cplexTerm + + datatype cplexProg = cplexProg of string + * cplexGoal + * ((string option * cplexConstr) + list) + * cplexBounds list + + datatype cplexResult = Unbounded + | Infeasible + | Undefined + | Optimal of string * + (((* name *) string * + (* value *) string) list) + + exception Load_cplexFile of string + exception Load_cplexResult of string + exception Save_cplexFile of string + exception Execute of string + + val load_cplexFile : string -> cplexProg + + val save_cplexFile : string -> cplexProg -> unit + + val elim_nonfree_bounds : cplexProg -> cplexProg + + val relax_strict_ineqs : cplexProg -> cplexProg + + val is_normed_cplexProg : cplexProg -> bool + + val load_cplexResult : string -> cplexResult + + val solve : cplexProg -> cplexResult +end; + +structure Cplex : CPLEX = +struct + +exception Load_cplexFile of string; +exception Load_cplexResult of string; +exception Save_cplexFile of string; + +datatype cplexTerm = cplexVar of string + | cplexNum of string + | cplexInf + | cplexNeg of cplexTerm + | cplexProd of cplexTerm * cplexTerm + | cplexSum of (cplexTerm list) +datatype cplexComp = cplexLe | cplexLeq | cplexEq | cplexGe | cplexGeq +datatype cplexGoal = cplexMinimize of cplexTerm | cplexMaximize of cplexTerm +datatype cplexConstr = cplexConstr of cplexComp * (cplexTerm * cplexTerm) +datatype cplexBounds = cplexBounds of cplexTerm * cplexComp * cplexTerm + * cplexComp * cplexTerm + | cplexBound of cplexTerm * cplexComp * cplexTerm +datatype cplexProg = cplexProg of string + * cplexGoal + * ((string option * cplexConstr) list) + * cplexBounds list + +fun rev_cmp cplexLe = cplexGe + | rev_cmp cplexLeq = cplexGeq + | rev_cmp cplexGe = cplexLe + | rev_cmp cplexGeq = cplexLeq + | rev_cmp cplexEq = cplexEq + +fun the None = raise (Load_cplexFile "Some expected") + | the (Some x) = x; + +fun modulo_signed is_something (cplexNeg u) = is_something u + | modulo_signed is_something u = is_something u + +fun is_Num (cplexNum _) = true + | is_Num _ = false + +fun is_Inf cplexInf = true + | is_Inf _ = false + +fun is_Var (cplexVar _) = true + | is_Var _ = false + +fun is_Neg (cplexNeg x ) = true + | is_Neg _ = false + +fun is_normed_Prod (cplexProd (t1, t2)) = + (is_Num t1) andalso (is_Var t2) + | is_normed_Prod x = is_Var x + +fun is_normed_Sum (cplexSum ts) = + (ts <> []) andalso forall (modulo_signed is_normed_Prod) ts + | is_normed_Sum x = modulo_signed is_normed_Prod x + +fun is_normed_Constr (cplexConstr (c, (t1, t2))) = + (is_normed_Sum t1) andalso (modulo_signed is_Num t2) + +fun is_Num_or_Inf x = is_Inf x orelse is_Num x + +fun is_normed_Bounds (cplexBounds (t1, c1, t2, c2, t3)) = + (c1 = cplexLe orelse c1 = cplexLeq) andalso + (c2 = cplexLe orelse c2 = cplexLeq) andalso + is_Var t2 andalso + modulo_signed is_Num_or_Inf t1 andalso + modulo_signed is_Num_or_Inf t3 + | is_normed_Bounds (cplexBound (t1, c, t2)) = + (is_Var t1 andalso (modulo_signed is_Num_or_Inf t2)) + orelse + (c <> cplexEq andalso + is_Var t2 andalso (modulo_signed is_Num_or_Inf t1)) + +fun term_of_goal (cplexMinimize x) = x + | term_of_goal (cplexMaximize x) = x + +fun is_normed_cplexProg (cplexProg (name, goal, constraints, bounds)) = + is_normed_Sum (term_of_goal goal) andalso + forall (fn (_,x) => is_normed_Constr x) constraints andalso + forall is_normed_Bounds bounds + +fun is_NL s = s = "\n" + +fun is_blank s = forall (fn c => c <> #"\n" andalso Char.isSpace c) (String.explode s) + +fun is_num a = + let + val b = String.explode a + fun num4 cs = forall Char.isDigit cs + fun num3 [] = true + | num3 (ds as (c::cs)) = + if c = #"+" orelse c = #"-" then + num4 cs + else + num4 ds + fun num2 [] = true + | num2 (c::cs) = + if c = #"e" orelse c = #"E" then num3 cs + else (Char.isDigit c) andalso num2 cs + fun num1 [] = true + | num1 (c::cs) = + if c = #"." then num2 cs + else if c = #"e" orelse c = #"E" then num3 cs + else (Char.isDigit c) andalso num1 cs + fun num [] = true + | num (c::cs) = + if c = #"." then num2 cs + else (Char.isDigit c) andalso num1 cs + in + num b + end + +fun is_delimiter s = s = "+" orelse s = "-" orelse s = ":" + +fun is_cmp s = s = "<" orelse s = ">" orelse s = "<=" + orelse s = ">=" orelse s = "=" + +fun is_symbol a = + let + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~" + fun is_symbol_char c = Char.isAlphaNum c orelse + exists (fn d => d=c) symbol_char + fun is_symbol_start c = is_symbol_char c andalso + not (Char.isDigit c) andalso + not (c= #".") + val b = String.explode a + in + b <> [] andalso is_symbol_start (hd b) andalso + forall is_symbol_char b + end + +fun to_upper s = String.implode (map Char.toUpper (String.explode s)) + +fun keyword x = + let + val a = to_upper x + in + if a = "BOUNDS" orelse a = "BOUND" then + Some "BOUNDS" + else if a = "MINIMIZE" orelse a = "MINIMUM" orelse a = "MIN" then + Some "MINIMIZE" + else if a = "MAXIMIZE" orelse a = "MAXIMUM" orelse a = "MAX" then + Some "MAXIMIZE" + else if a = "ST" orelse a = "S.T." orelse a = "ST." then + Some "ST" + else if a = "FREE" orelse a = "END" then + Some a + else if a = "GENERAL" orelse a = "GENERALS" orelse a = "GEN" then + Some "GENERAL" + else if a = "INTEGER" orelse a = "INTEGERS" orelse a = "INT" then + Some "INTEGER" + else if a = "BINARY" orelse a = "BINARIES" orelse a = "BIN" then + Some "BINARY" + else if a = "INF" orelse a = "INFINITY" then + Some "INF" + else + None + end + +val TOKEN_ERROR = ~1 +val TOKEN_BLANK = 0 +val TOKEN_NUM = 1 +val TOKEN_DELIMITER = 2 +val TOKEN_SYMBOL = 3 +val TOKEN_LABEL = 4 +val TOKEN_CMP = 5 +val TOKEN_KEYWORD = 6 +val TOKEN_NL = 7 + +(* tokenize takes a list of chars as argument and returns a list of + int * string pairs, each string representing a "cplex token", + and each int being one of TOKEN_NUM, TOKEN_DELIMITER, TOKEN_CMP + or TOKEN_SYMBOL *) +fun tokenize s = + let + val flist = [(is_NL, TOKEN_NL), + (is_blank, TOKEN_BLANK), + (is_num, TOKEN_NUM), + (is_delimiter, TOKEN_DELIMITER), + (is_cmp, TOKEN_CMP), + (is_symbol, TOKEN_SYMBOL)] + fun match_helper [] s = (fn x => false, TOKEN_ERROR) + | match_helper (f::fs) s = + if ((fst f) s) then f else match_helper fs s + fun match s = match_helper flist s + fun tok s = + if s = "" then [] else + let + val h = String.substring (s,0,1) + val (f, j) = match h + fun len i = + if size s = i then i + else if f (String.substring (s,0,i+1)) then + len (i+1) + else i + in + if j < 0 then + (if h = "\\" then [] + else raise (Load_cplexFile ("token expected, found: " + ^s))) + else + let + val l = len 1 + val u = String.substring (s,0,l) + val v = String.extract (s,l,NONE) + in + if j = 0 then tok v else (j, u) :: tok v + end + end + in + tok s + end + +exception Tokenize of string; + +fun tokenize_general flist s = + let + fun match_helper [] s = raise (Tokenize s) + | match_helper (f::fs) s = + if ((fst f) s) then f else match_helper fs s + fun match s = match_helper flist s + fun tok s = + if s = "" then [] else + let + val h = String.substring (s,0,1) + val (f, j) = match h + fun len i = + if size s = i then i + else if f (String.substring (s,0,i+1)) then + len (i+1) + else i + val l = len 1 + in + (j, String.substring (s,0,l)) :: tok (String.extract (s,l,NONE)) + end + in + tok s + end + +fun load_cplexFile name = + let + val f = TextIO.openIn name + val ignore_NL = ref true + val rest = ref [] + + fun is_symbol s c = (fst c) = TOKEN_SYMBOL andalso (to_upper (snd c)) = s + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + Some u + ) + end + else + let val s = TextIO.inputLine f in + if s = "" then None else + let val t = tokenize s in + if (length t >= 2 andalso + snd(hd (tl t)) = ":") + then + rest := (TOKEN_LABEL, snd (hd t)) :: (tl (tl t)) + else if (length t >= 2) andalso is_symbol "SUBJECT" (hd (t)) + andalso is_symbol "TO" (hd (tl t)) + then + rest := (TOKEN_SYMBOL, "ST") :: (tl (tl t)) + else + rest := t; + readToken_helper () + end + end + + fun readToken_helper2 () = + let val c = readToken_helper () in + if c = None then None + else if !ignore_NL andalso fst (the c) = TOKEN_NL then + readToken_helper2 () + else if fst (the c) = TOKEN_SYMBOL + andalso keyword (snd (the c)) <> None + then Some (TOKEN_KEYWORD, the (keyword (snd (the c)))) + else c + end + + fun readToken () = readToken_helper2 () + + fun pushToken a = rest := (a::(!rest)) + + fun is_value token = + fst token = TOKEN_NUM orelse (fst token = TOKEN_KEYWORD + andalso snd token = "INF") + + fun get_value token = + if fst token = TOKEN_NUM then + cplexNum (snd token) + else if fst token = TOKEN_KEYWORD andalso snd token = "INF" + then + cplexInf + else + raise (Load_cplexFile "num expected") + + fun readTerm_Product only_num = + let val c = readToken () in + if c = None then None + else if fst (the c) = TOKEN_SYMBOL + then ( + if only_num then (pushToken (the c); None) + else Some (cplexVar (snd (the c))) + ) + else if only_num andalso is_value (the c) then + Some (get_value (the c)) + else if is_value (the c) then + let val t1 = get_value (the c) + val d = readToken () + in + if d = None then Some t1 + else if fst (the d) = TOKEN_SYMBOL then + Some (cplexProd (t1, cplexVar (snd (the d)))) + else + (pushToken (the d); Some t1) + end + else (pushToken (the c); None) + end + + fun readTerm_Signed only_signed only_num = + let + val c = readToken () + in + if c = None then None + else + let val d = the c in + if d = (TOKEN_DELIMITER, "+") then + readTerm_Product only_num + else if d = (TOKEN_DELIMITER, "-") then + Some (cplexNeg (the (readTerm_Product + only_num))) + else (pushToken d; + if only_signed then None + else readTerm_Product only_num) + end + end + + fun readTerm_Sum first_signed = + let val c = readTerm_Signed first_signed false in + if c = None then [] else (the c)::(readTerm_Sum true) + end + + fun readTerm () = + let val c = readTerm_Sum false in + if c = [] then None + else if tl c = [] then Some (hd c) + else Some (cplexSum c) + end + + fun readLabeledTerm () = + let val c = readToken () in + if c = None then (None, None) + else if fst (the c) = TOKEN_LABEL then + let val t = readTerm () in + if t = None then + raise (Load_cplexFile ("term after label "^ + (snd (the c))^ + " expected")) + else (Some (snd (the c)), t) + end + else (pushToken (the c); (None, readTerm ())) + end + + fun readGoal () = + let + val g = readToken () + in + if g = Some (TOKEN_KEYWORD, "MAXIMIZE") then + cplexMaximize (the (snd (readLabeledTerm ()))) + else if g = Some (TOKEN_KEYWORD, "MINIMIZE") then + cplexMinimize (the (snd (readLabeledTerm ()))) + else raise (Load_cplexFile "MAXIMIZE or MINIMIZE expected") + end + + fun str2cmp b = + (case b of + "<" => cplexLe + | "<=" => cplexLeq + | ">" => cplexGe + | ">=" => cplexGeq + | "=" => cplexEq + | _ => raise (Load_cplexFile (b^" is no TOKEN_CMP"))) + + fun readConstraint () = + let + val t = readLabeledTerm () + fun make_constraint b t1 t2 = + cplexConstr + (str2cmp b, + (t1, t2)) + in + if snd t = None then None + else + let val c = readToken () in + if c = None orelse fst (the c) <> TOKEN_CMP + then raise (Load_cplexFile "TOKEN_CMP expected") + else + let val n = readTerm_Signed false true in + if n = None then + raise (Load_cplexFile "num expected") + else + Some (fst t, + make_constraint (snd (the c)) + (the (snd t)) + (the n)) + end + end + end + + fun readST () = + let + fun readbody () = + let val t = readConstraint () in + if t = None then [] + else if (is_normed_Constr (snd (the t))) then + (the t)::(readbody ()) + else if (fst (the t) <> None) then + raise (Load_cplexFile + ("constraint '"^(the (fst (the t)))^ + "'is not normed")) + else + raise (Load_cplexFile + "constraint is not normed") + end + in + if readToken () = Some (TOKEN_KEYWORD, "ST") + then + readbody () + else + raise (Load_cplexFile "ST expected") + end + + fun readCmp () = + let val c = readToken () in + if c = None then None + else if fst (the c) = TOKEN_CMP then + Some (str2cmp (snd (the c))) + else (pushToken (the c); None) + end + + fun skip_NL () = + let val c = readToken () in + if c <> None andalso fst (the c) = TOKEN_NL then + skip_NL () + else + (pushToken (the c); ()) + end + + fun is_var (cplexVar _) = true + | is_var _ = false + + fun make_bounds c t1 t2 = + cplexBound (t1, c, t2) + + fun readBound () = + let + val _ = skip_NL () + val t1 = readTerm () + in + if t1 = None then None + else + let + val c1 = readCmp () + in + if c1 = None then + let + val c = readToken () + in + if c = Some (TOKEN_KEYWORD, "FREE") then + Some ( + cplexBounds (cplexNeg cplexInf, + cplexLeq, + the t1, + cplexLeq, + cplexInf)) + else + raise (Load_cplexFile "FREE expected") + end + else + let + val t2 = readTerm () + in + if t2 = None then + raise (Load_cplexFile "term expected") + else + let val c2 = readCmp () in + if c2 = None then + Some (make_bounds (the c1) + (the t1) + (the t2)) + else + Some ( + cplexBounds (the t1, + the c1, + the t2, + the c2, + the (readTerm()))) + end + end + end + end + + fun readBounds () = + let + fun makestring b = "?" + fun readbody () = + let + val b = readBound () + in + if b = None then [] + else if (is_normed_Bounds (the b)) then + (the b)::(readbody()) + else ( + raise (Load_cplexFile + ("bounds are not normed in: "^ + (makestring (the b))))) + end + in + if readToken () = Some (TOKEN_KEYWORD, "BOUNDS") then + readbody () + else raise (Load_cplexFile "BOUNDS expected") + end + + fun readEnd () = + if readToken () = Some (TOKEN_KEYWORD, "END") then () + else raise (Load_cplexFile "END expected") + + val result_Goal = readGoal () + val result_ST = readST () + val _ = ignore_NL := false + val result_Bounds = readBounds () + val _ = ignore_NL := true + val _ = readEnd () + val _ = TextIO.closeIn f + in + cplexProg (name, result_Goal, result_ST, result_Bounds) + end + +fun save_cplexFile filename (cplexProg (name, goal, constraints, bounds)) = + let + val f = TextIO.openOut filename + + fun basic_write s = TextIO.output(f, s) + + val linebuf = ref "" + fun buf_flushline s = + (basic_write (!linebuf); + basic_write "\n"; + linebuf := s) + fun buf_add s = linebuf := (!linebuf) ^ s + + fun write s = + if (String.size s) + (String.size (!linebuf)) >= 250 then + buf_flushline (" "^s) + else + buf_add s + + fun writeln s = (buf_add s; buf_flushline "") + + fun write_term (cplexVar x) = write x + | write_term (cplexNum x) = write x + | write_term cplexInf = write "inf" + | write_term (cplexProd (cplexNum "1", b)) = write_term b + | write_term (cplexProd (a, b)) = + (write_term a; write " "; write_term b) + | write_term (cplexNeg x) = (write " - "; write_term x) + | write_term (cplexSum ts) = write_terms ts + and write_terms [] = () + | write_terms (t::ts) = + (if (not (is_Neg t)) then write " + " else (); + write_term t; write_terms ts) + + fun write_goal (cplexMaximize term) = + (writeln "MAXIMIZE"; write_term term; writeln "") + | write_goal (cplexMinimize term) = + (writeln "MINIMIZE"; write_term term; writeln "") + + fun write_cmp cplexLe = write "<" + | write_cmp cplexLeq = write "<=" + | write_cmp cplexEq = write "=" + | write_cmp cplexGe = write ">" + | write_cmp cplexGeq = write ">=" + + fun write_constr (cplexConstr (cmp, (a,b))) = + (write_term a; + write " "; + write_cmp cmp; + write " "; + write_term b) + + fun write_constraints [] = () + | write_constraints (c::cs) = + (if (fst c <> None) + then + (write (the (fst c)); write ": ") + else + (); + write_constr (snd c); + writeln ""; + write_constraints cs) + + fun write_bounds [] = () + | write_bounds ((cplexBounds (t1,c1,t2,c2,t3))::bs) = + ((if t1 = cplexNeg cplexInf andalso t3 = cplexInf + andalso (c1 = cplexLeq orelse c1 = cplexLe) + andalso (c2 = cplexLeq orelse c2 = cplexLe) + then + (write_term t2; write " free") + else + (write_term t1; write " "; write_cmp c1; write " "; + write_term t2; write " "; write_cmp c2; write " "; + write_term t3) + ); writeln ""; write_bounds bs) + | write_bounds ((cplexBound (t1, c, t2)) :: bs) = + (write_term t1; write " "; + write_cmp c; write " "; + write_term t2; writeln ""; write_bounds bs) + + val _ = write_goal goal + val _ = (writeln ""; writeln "ST") + val _ = write_constraints constraints + val _ = (writeln ""; writeln "BOUNDS") + val _ = write_bounds bounds + val _ = (writeln ""; writeln "END") + val _ = TextIO.closeOut f + in + () + end + +fun norm_Constr (constr as cplexConstr (c, (t1, t2))) = + if not (modulo_signed is_Num t2) andalso + modulo_signed is_Num t1 + then + [cplexConstr (rev_cmp c, (t2, t1))] + else if (c = cplexLe orelse c = cplexLeq) andalso + (t1 = (cplexNeg cplexInf) orelse t2 = cplexInf) + then + [] + else if (c = cplexGe orelse c = cplexGeq) andalso + (t1 = cplexInf orelse t2 = cplexNeg cplexInf) + then + [] + else + [constr] + +fun bound2constr (cplexBounds (t1,c1,t2,c2,t3)) = + (norm_Constr(cplexConstr (c1, (t1, t2)))) + @ (norm_Constr(cplexConstr (c2, (t2, t3)))) + | bound2constr (cplexBound (t1, cplexEq, t2)) = + (norm_Constr(cplexConstr (cplexLeq, (t1, t2)))) + @ (norm_Constr(cplexConstr (cplexLeq, (t2, t1)))) + | bound2constr (cplexBound (t1, c1, t2)) = + norm_Constr(cplexConstr (c1, (t1,t2))) + +val emptyset = Symtab.empty + +fun singleton v = Symtab.update ((v, ()), emptyset) + +fun merge a b = Symtab.merge (op =) (a, b) + +fun mergemap f ts = foldl + (fn (table, x) => merge table (f x)) + (Symtab.empty, ts) + +fun diff a b = Symtab.foldl (fn (a, (k,v)) => + (Symtab.delete k a) handle UNDEF => a) + (a, b) + +fun collect_vars (cplexVar v) = singleton v + | collect_vars (cplexNeg t) = collect_vars t + | collect_vars (cplexProd (t1, t2)) = + merge (collect_vars t1) (collect_vars t2) + | collect_vars (cplexSum ts) = mergemap collect_vars ts + | collect_vars _ = emptyset + +(* Eliminates all nonfree bounds from the linear program and produces an + equivalent program with only free bounds + IF for the input program P holds: is_normed_cplexProg P *) +fun elim_nonfree_bounds (cplexProg (name, goal, constraints, bounds)) = + let + fun collect_constr_vars (_, cplexConstr (c, (t1,_))) = + (collect_vars t1) + + val cvars = merge (collect_vars (term_of_goal goal)) + (mergemap collect_constr_vars constraints) + + fun collect_lower_bounded_vars + (cplexBounds (t1, c1, cplexVar v, c2, t3)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (_, cplexLe, cplexVar v)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (_, cplexLeq, cplexVar v)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexGe,_)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexGeq, _)) = + singleton v + | collect_lower_bounded_vars + (cplexBound (cplexVar v, cplexEq, _)) = + singleton v + | collect_lower_bounded_vars _ = emptyset + + val lvars = mergemap collect_lower_bounded_vars bounds + val positive_vars = diff cvars lvars + val zero = cplexNum "0" + + fun make_pos_constr v = + (None, cplexConstr (cplexGeq, ((cplexVar v), zero))) + + fun make_free_bound v = + cplexBounds (cplexNeg cplexInf, cplexLeq, + cplexVar v, cplexLeq, + cplexInf) + + val pos_constrs = rev(Symtab.foldl + (fn (l, (k,v)) => (make_pos_constr k) :: l) + ([], positive_vars)) + val bound_constrs = map (fn x => (None, x)) + (flat (map bound2constr bounds)) + val constraints' = constraints @ pos_constrs @ bound_constrs + val bounds' = rev(Symtab.foldl (fn (l, (v,_)) => + (make_free_bound v)::l) + ([], cvars)) + in + cplexProg (name, goal, constraints', bounds') + end + +fun relax_strict_ineqs (cplexProg (name, goals, constrs, bounds)) = + let + fun relax cplexLe = cplexLeq + | relax cplexGe = cplexGeq + | relax x = x + + fun relax_constr (n, cplexConstr(c, (t1, t2))) = + (n, cplexConstr(relax c, (t1, t2))) + + fun relax_bounds (cplexBounds (t1, c1, t2, c2, t3)) = + cplexBounds (t1, relax c1, t2, relax c2, t3) + | relax_bounds (cplexBound (t1, c, t2)) = + cplexBound (t1, relax c, t2) + in + cplexProg (name, + goals, + map relax_constr constrs, + map relax_bounds bounds) + end + +datatype cplexResult = Unbounded + | Infeasible + | Undefined + | Optimal of string * ((string * string) list) + +fun is_separator x = forall (fn c => c = #"-") (String.explode x) + +fun is_sign x = (x = "+" orelse x = "-") + +fun is_colon x = (x = ":") + +fun is_resultsymbol a = + let + val symbol_char = String.explode "!\"#$%&()/,.;?@_`'{}|~-" + fun is_symbol_char c = Char.isAlphaNum c orelse + exists (fn d => d=c) symbol_char + fun is_symbol_start c = is_symbol_char c andalso + not (Char.isDigit c) andalso + not (c= #".") andalso + not (c= #"-") + val b = String.explode a + in + b <> [] andalso is_symbol_start (hd b) andalso + forall is_symbol_char b + end + +val TOKEN_SIGN = 100 +val TOKEN_COLON = 101 +val TOKEN_SEPARATOR = 102 + +fun load_cplexResult name = + let + val flist = [(is_NL, TOKEN_NL), + (is_blank, TOKEN_BLANK), + (is_num, TOKEN_NUM), + (is_sign, TOKEN_SIGN), + (is_colon, TOKEN_COLON), + (is_cmp, TOKEN_CMP), + (is_resultsymbol, TOKEN_SYMBOL), + (is_separator, TOKEN_SEPARATOR)] + + val tokenize = tokenize_general flist + + val f = TextIO.openIn name + + val rest = ref [] + + fun readToken_helper () = + if length (!rest) > 0 then + let val u = hd (!rest) in + ( + rest := tl (!rest); + Some u + ) + end + else + let val s = TextIO.inputLine f in + if s = "" then None else (rest := tokenize s; readToken_helper()) + end + + fun is_tt tok ty = (tok <> None andalso (fst (the tok)) = ty) + + fun pushToken a = if a = None then () else (rest := ((the a)::(!rest))) + + fun readToken () = + let val t = readToken_helper () in + if is_tt t TOKEN_BLANK then + readToken () + else if is_tt t TOKEN_NL then + let val t2 = readToken_helper () in + if is_tt t2 TOKEN_SIGN then + (pushToken (Some (TOKEN_SEPARATOR, "-")); t) + else + (pushToken t2; t) + end + else if is_tt t TOKEN_SIGN then + let val t2 = readToken_helper () in + if is_tt t2 TOKEN_NUM then + (Some (TOKEN_NUM, (snd (the t))^(snd (the t2)))) + else + (pushToken t2; t) + end + else + t + end + + fun readRestOfLine P = + let + val t = readToken () + in + if is_tt t TOKEN_NL orelse t = None + then P + else readRestOfLine P + end + + fun readHeader () = + let + fun readStatus () = readRestOfLine ("STATUS", snd (the (readToken ()))) + fun readObjective () = readRestOfLine ("OBJECTIVE", snd (the (readToken (); readToken (); readToken ()))) + val t1 = readToken () + val t2 = readToken () + in + if is_tt t1 TOKEN_SYMBOL andalso is_tt t2 TOKEN_COLON + then + case to_upper (snd (the t1)) of + "STATUS" => (readStatus ())::(readHeader ()) + | "OBJECTIVE" => (readObjective())::(readHeader ()) + | _ => (readRestOfLine (); readHeader ()) + else + (pushToken t2; pushToken t1; []) + end + + fun skip_until_sep () = + let val x = readToken () in + if is_tt x TOKEN_SEPARATOR then + readRestOfLine () + else + skip_until_sep () + end + + fun load_value () = + let + val t1 = readToken () + val t2 = readToken () + in + if is_tt t1 TOKEN_NUM andalso is_tt t2 TOKEN_SYMBOL then + let + val t = readToken () + val state = if is_tt t TOKEN_NL then readToken () else t + val _ = if is_tt state TOKEN_SYMBOL then () else raise (Load_cplexResult "state expected") + val k = readToken () + in + if is_tt k TOKEN_NUM then + readRestOfLine (Some (snd (the t2), snd (the k))) + else + raise (Load_cplexResult "number expected") + end + else + (pushToken t2; pushToken t1; None) + end + + fun load_values () = + let val v = load_value () in + if v = None then [] else (the v)::(load_values ()) + end + + val header = readHeader () + + val result = + case assoc (header, "STATUS") of + Some "INFEASIBLE" => Infeasible + | Some "UNBOUNDED" => Unbounded + | Some "OPTIMAL" => Optimal (the (assoc (header, "OBJECTIVE")), + (skip_until_sep (); + skip_until_sep (); + load_values ())) + | _ => Undefined + + val _ = TextIO.closeIn f + in + result + end + handle (Tokenize s) => raise (Load_cplexResult ("Tokenize: "^s)) + | OPTION => raise (Load_cplexResult "OPTION") + | x => raise x + +exception Execute of string; + +fun execute_cplex lpname resultname = +let + fun wrap s = "\""^s^"\"" + val cplex_path = getenv "CPLEX" + val cplex = if cplex_path = "" then "glpsol" else cplex_path + val command = (wrap cplex)^" --lpt "^(wrap lpname)^" --output "^(wrap resultname) +in + execute command +end + +fun tmp_file s = Path.pack (Path.expand (File.tmp_path (Path.make [s]))) + +fun solve prog = + let + val name = LargeInt.toString (Time.toMicroseconds (Time.now ())) + val lpname = tmp_file (name^".lp") + val resultname = tmp_file (name^".result") + val _ = save_cplexFile lpname prog + val answer = execute_cplex lpname resultname + in + let + val result = load_cplexResult resultname + val _ = OS.FileSys.remove lpname + val _ = OS.FileSys.remove resultname + in + result + end + handle (Load_cplexResult s) => raise (Execute ("Load_cplexResult: "^s^"\nExecute: "^answer)) + | _ => raise (Execute answer) + end +end; + +val demofile = "/home/obua/flyspeck/kepler/LP/cplexPent2.lp45" +val demoout = "/home/obua/flyspeck/kepler/LP/test.out" +val demoresult = "/home/obua/flyspeck/kepler/LP/try/test2.sol" + +fun loadcplex () = Cplex.relax_strict_ineqs + (Cplex.load_cplexFile demofile) + +fun writecplex lp = Cplex.save_cplexFile demoout lp + +fun test () = + let + val lp = loadcplex () + val lp2 = Cplex.elim_nonfree_bounds lp + in + writecplex lp2 + end + +fun loadresult () = Cplex.load_cplexResult demoresult; diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/CplexMatrixConverter.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/CplexMatrixConverter.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,132 @@ +signature MATRIX_BUILDER = +sig + type vector + type matrix + + val empty_vector : vector + val empty_matrix : matrix + + exception Nat_expected of int + val set_elem : vector -> int -> string -> vector + val set_vector : matrix -> int -> vector -> matrix +end; + +signature CPLEX_MATRIX_CONVERTER = +sig + structure cplex : CPLEX + structure matrix_builder : MATRIX_BUILDER + type vector = matrix_builder.vector + type matrix = matrix_builder.matrix + type naming = int * (int -> string) * (string -> int) + + exception Converter of string + + (* program must fulfill is_normed_cplexProg and must be an element of the image of elim_nonfree_bounds *) + (* convert_prog maximize c A b naming *) + val convert_prog : cplex.cplexProg -> bool * vector * matrix * vector * naming + + (* results must be optimal, converts_results returns the optimal value as string and the solution as vector *) + (* convert_results results name2index *) + val convert_results : cplex.cplexResult -> (string -> int) -> string * vector +end; + +functor MAKE_CPLEX_MATRIX_CONVERTER (structure cplex: CPLEX and matrix_builder: MATRIX_BUILDER) : CPLEX_MATRIX_CONVERTER = +struct + +structure cplex = cplex +structure matrix_builder = matrix_builder +type matrix = matrix_builder.matrix +type vector = matrix_builder.vector +type naming = int * (int -> string) * (string -> int) + +open matrix_builder +open cplex + +exception Converter of string; + +structure Inttab = TableFun(type key = int val ord = int_ord); + +fun neg_term (cplexNeg t) = t + | neg_term (cplexSum ts) = cplexSum (map neg_term ts) + | neg_term t = cplexNeg t + +fun convert_prog (cplexProg (s, goal, constrs, bounds)) = + let + fun build_naming index i2s s2i [] = (index, i2s, s2i) + | build_naming index i2s s2i (cplexBounds (cplexNeg cplexInf, cplexLeq, cplexVar v, cplexLeq, cplexInf)::bounds) + = build_naming (index+1) (Inttab.update ((index, v), i2s)) (Symtab.update_new ((v, index), s2i)) bounds + | build_naming _ _ _ _ = raise (Converter "nonfree bound") + + val (varcount, i2s_tab, s2i_tab) = build_naming 0 Inttab.empty Symtab.empty bounds + + fun i2s i = case Inttab.lookup (i2s_tab, i) of None => raise (Converter "index not found") + | Some n => n + fun s2i s = case Symtab.lookup (s2i_tab, s) of None => raise (Converter ("name not found: "^s)) + | Some i => i + fun num2str positive (cplexNeg t) = num2str (not positive) t + | num2str positive (cplexNum num) = if positive then num else "-"^num + | num2str _ _ = raise (Converter "term is not a (possibly signed) number") + + fun setprod vec positive (cplexNeg t) = setprod vec (not positive) t + | setprod vec positive (cplexVar v) = set_elem vec (s2i v) (if positive then "1" else "-1") + | setprod vec positive (cplexProd (cplexNum num, cplexVar v)) = + set_elem vec (s2i v) (if positive then num else "-"^num) + | setprod _ _ _ = raise (Converter "term is not a normed product") + + fun sum2vec (cplexSum ts) = foldl (fn (vec, t) => setprod vec true t) (empty_vector, ts) + | sum2vec t = setprod empty_vector true t + + fun constrs2Ab j A b [] = (A, b) + | constrs2Ab j A b ((_, cplexConstr (cplexLeq, (t1,t2)))::cs) = + constrs2Ab (j+1) (set_vector A j (sum2vec t1)) (set_elem b j (num2str true t2)) cs + | constrs2Ab j A b ((_, cplexConstr (cplexGeq, (t1,t2)))::cs) = + constrs2Ab (j+1) (set_vector A j (sum2vec (neg_term t1))) (set_elem b j (num2str true (neg_term t2))) cs + | constrs2Ab j A b ((_, cplexConstr (cplexEq, (t1,t2)))::cs) = + constrs2Ab j A b ((None, cplexConstr (cplexLeq, (t1,t2))):: + (None, cplexConstr (cplexGeq, (t1, t2)))::cs) + | constrs2Ab _ _ _ _ = raise (Converter "no strict constraints allowed") + + val (A, b) = constrs2Ab 0 empty_matrix empty_vector constrs + + val (goal_maximize, goal_term) = + case goal of + (cplexMaximize t) => (true, t) + | (cplexMinimize t) => (false, t) + in + (goal_maximize, sum2vec goal_term, A, b, (varcount, i2s, s2i)) + end + +fun convert_results (cplex.Optimal (opt, entries)) name2index = + let + fun setv (v, (name, value)) = (matrix_builder.set_elem v (name2index name) value) + in + (opt, foldl setv (matrix_builder.empty_vector, entries)) + end + | convert_results _ _ = raise (Converter "No optimal result") + + +end; + +structure SimpleMatrixBuilder : MATRIX_BUILDER = +struct +type vector = (int * string) list +type matrix = (int * vector) list + +val empty_matrix = [] +val empty_vector = [] + +exception Nat_expected of int; + +fun set_elem v i s = v @ [(i, s)] + +fun set_vector m i v = m @ [(i, v)] + +end; + + + +structure SimpleCplexMatrixConverter = MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = SimpleMatrixBuilder); + + + + diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/ExactFloatingPoint.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/ExactFloatingPoint.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,215 @@ +structure ExactFloatingPoint : +sig + exception Destruct_floatstr of string + + val destruct_floatstr : (char -> bool) -> (char -> bool) -> string -> bool * string * string * bool * string + + exception Floating_point of string + + type floatrep = IntInf.int * IntInf.int + val approx_dec_by_bin : IntInf.int -> floatrep -> floatrep * floatrep + val approx_decstr_by_bin : int -> string -> floatrep * floatrep +end += +struct + +fun fst (a,b) = a +fun snd (a,b) = b + +val filter = List.filter; + +exception Destruct_floatstr of string; + +fun destruct_floatstr isDigit isExp number = + let + val numlist = filter (not o Char.isSpace) (String.explode number) + + fun countsigns ((#"+")::cs) = countsigns cs + | countsigns ((#"-")::cs) = + let + val (positive, rest) = countsigns cs + in + (not positive, rest) + end + | countsigns cs = (true, cs) + + fun readdigits [] = ([], []) + | readdigits (q as c::cs) = + if (isDigit c) then + let + val (digits, rest) = readdigits cs + in + (c::digits, rest) + end + else + ([], q) + + fun readfromexp_helper cs = + let + val (positive, rest) = countsigns cs + val (digits, rest') = readdigits rest + in + case rest' of + [] => (positive, digits) + | _ => raise (Destruct_floatstr number) + end + + fun readfromexp [] = (true, []) + | readfromexp (c::cs) = + if isExp c then + readfromexp_helper cs + else + raise (Destruct_floatstr number) + + fun readfromdot [] = ([], readfromexp []) + | readfromdot ((#".")::cs) = + let + val (digits, rest) = readdigits cs + val exp = readfromexp rest + in + (digits, exp) + end + | readfromdot cs = readfromdot ((#".")::cs) + + val (positive, numlist) = countsigns numlist + val (digits1, numlist) = readdigits numlist + val (digits2, exp) = readfromdot numlist + in + (positive, String.implode digits1, String.implode digits2, fst exp, String.implode (snd exp)) + end + +type floatrep = IntInf.int * IntInf.int + +exception Floating_point of string; + +val ln2_10 = (Math.ln 10.0)/(Math.ln 2.0) + +fun intmul a b = IntInf.* (a,b) +fun intsub a b = IntInf.- (a,b) +fun intadd a b = IntInf.+ (a,b) +fun intpow a b = IntInf.pow (a, IntInf.toInt b); +fun intle a b = IntInf.<= (a, b); +fun intless a b = IntInf.< (a, b); +fun intneg a = IntInf.~ a; +val zero = IntInf.fromInt 0; +val one = IntInf.fromInt 1; +val two = IntInf.fromInt 2; +val ten = IntInf.fromInt 10; +val five = IntInf.fromInt 5; + +fun find_most_significant q r = + let + fun int2real i = + case Real.fromString (IntInf.toString i) of + SOME r => r + | NONE => raise (Floating_point "int2real") + fun subtract (q, r) (q', r') = + if intle r r' then + (intsub q (intmul q' (intpow ten (intsub r' r))), r) + else + (intsub (intmul q (intpow ten (intsub r r'))) q', r') + fun bin2dec d = + if intle zero d then + (intpow two d, zero) + else + (intpow five (intneg d), d) + + val L = IntInf.fromInt (Real.floor (int2real (IntInf.fromInt (IntInf.log2 q)) + (int2real r) * ln2_10)) + val L1 = intadd L one + + val (q1, r1) = subtract (q, r) (bin2dec L1) + in + if intle zero q1 then + let + val (q2, r2) = subtract (q, r) (bin2dec (intadd L1 one)) + in + if intle zero q2 then + raise (Floating_point "find_most_significant") + else + (L1, (q1, r1)) + end + else + let + val (q0, r0) = subtract (q, r) (bin2dec L) + in + if intle zero q0 then + (L, (q0, r0)) + else + raise (Floating_point "find_most_significant") + end + end + +fun approx_dec_by_bin n (q,r) = + let + fun addseq acc d' [] = acc + | addseq acc d' (d::ds) = addseq (intadd acc (intpow two (intsub d d'))) d' ds + + fun seq2bin [] = (zero, zero) + | seq2bin (d::ds) = (intadd (addseq zero d ds) one, d) + + fun approx d_seq d0 precision (q,r) = + if q = zero then + let val x = seq2bin d_seq in + (x, x) + end + else + let + val (d, (q', r')) = find_most_significant q r + in + if intless precision (intsub d0 d) then + let + val d' = intsub d0 precision + val x1 = seq2bin (d_seq) + val x2 = (intadd (intmul (fst x1) (intpow two (intsub (snd x1) d'))) one, d') (* = seq2bin (d'::d_seq) *) + in + (x1, x2) + end + else + approx (d::d_seq) d0 precision (q', r') + end + + fun approx_start precision (q, r) = + if q = zero then + ((zero, zero), (zero, zero)) + else + let + val (d, (q', r')) = find_most_significant q r + in + if intle precision zero then + let + val x1 = seq2bin [d] + in + if q' = zero then + (x1, x1) + else + (x1, seq2bin [intadd d one]) + end + else + approx [d] d precision (q', r') + end + in + if intle zero q then + approx_start n (q,r) + else + let + val ((a1,b1), (a2, b2)) = approx_start n (intneg q, r) + in + ((intneg a2, b2), (intneg a1, b1)) + end + end + +fun approx_decstr_by_bin n decstr = + let + fun str2int s = case IntInf.fromString s of SOME x => x | NONE => zero + fun signint p x = if p then x else intneg x + + val (p, d1, d2, ep, e) = destruct_floatstr Char.isDigit (fn e => e = #"e" orelse e = #"E") decstr + val s = IntInf.fromInt (size d2) + + val q = signint p (intadd (intmul (str2int d1) (intpow ten s)) (str2int d2)) + val r = intsub (signint ep (str2int e)) s + in + approx_dec_by_bin (IntInf.fromInt n) (q,r) + end + +end; diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/Float.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/Float.thy Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,508 @@ +theory Float = Hyperreal: + +constdefs + pow2 :: "int \ real" + "pow2 a == if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a))))" + float :: "int * int \ real" + "float x == (real (fst x)) * (pow2 (snd x))" + +lemma pow2_0[simp]: "pow2 0 = 1" +by (simp add: pow2_def) + +lemma pow2_1[simp]: "pow2 1 = 2" +by (simp add: pow2_def) + +lemma pow2_neg: "pow2 x = inverse (pow2 (-x))" +by (simp add: pow2_def) + +lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)" +proof - + have h: "! n. nat (2 + int n) - Suc 0 = nat (1 + int n)" by arith + have g: "! a b. a - -1 = a + (1::int)" by arith + have pos: "! n. pow2 (int n + 1) = 2 * pow2 (int n)" + apply (auto, induct_tac n) + apply (simp_all add: pow2_def) + apply (rule_tac m1="2" and n1="nat (2 + int na)" in ssubst[OF realpow_num_eq_if]) + apply (auto simp add: h) + apply arith + done + show ?thesis + proof (induct a) + case (1 n) + from pos show ?case by (simp add: ring_eq_simps) + next + case (2 n) + show ?case + apply (auto) + apply (subst pow2_neg[of "- int n"]) + apply (subst pow2_neg[of "-1 - int n"]) + apply (auto simp add: g pos) + done + qed +qed + +lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)" +proof (induct b) + case (1 n) + show ?case + proof (induct n) + case 0 + show ?case by simp + next + case (Suc m) + show ?case by (auto simp add: ring_eq_simps pow2_add1 prems) + qed +next + case (2 n) + show ?case + proof (induct n) + case 0 + show ?case + apply (auto) + apply (subst pow2_neg[of "a + -1"]) + apply (subst pow2_neg[of "-1"]) + apply (simp) + apply (insert pow2_add1[of "-a"]) + apply (simp add: ring_eq_simps) + apply (subst pow2_neg[of "-a"]) + apply (simp) + done + case (Suc m) + have a: "int m - (a + -2) = 1 + (int m - a + 1)" by arith + have b: "int m - -2 = 1 + (int m + 1)" by arith + show ?case + apply (auto) + apply (subst pow2_neg[of "a + (-2 - int m)"]) + apply (subst pow2_neg[of "-2 - int m"]) + apply (auto simp add: ring_eq_simps) + apply (subst a) + apply (subst b) + apply (simp only: pow2_add1) + apply (subst pow2_neg[of "int m - a + 1"]) + apply (subst pow2_neg[of "int m + 1"]) + apply auto + apply (insert prems) + apply (auto simp add: ring_eq_simps) + done + qed +qed + +lemma "float (a, e) + float (b, e) = float (a + b, e)" +by (simp add: float_def ring_eq_simps) + +constdefs + int_of_real :: "real \ int" + "int_of_real x == SOME y. real y = x" + real_is_int :: "real \ bool" + "real_is_int x == ? (u::int). x = real u" + +lemma real_is_int_def2: "real_is_int x = (x = real (int_of_real x))" +by (auto simp add: real_is_int_def int_of_real_def) + +lemma float_transfer: "real_is_int ((real a)*(pow2 c)) \ float (a, b) = float (int_of_real ((real a)*(pow2 c)), b - c)" +by (simp add: float_def real_is_int_def2 pow2_add[symmetric]) + +lemma pow2_int: "pow2 (int c) = (2::real)^c" +by (simp add: pow2_def) + +lemma float_transfer_nat: "float (a, b) = float (a * 2^c, b - int c)" +by (simp add: float_def pow2_int[symmetric] pow2_add[symmetric]) + +lemma real_is_int_real[simp]: "real_is_int (real (x::int))" +by (auto simp add: real_is_int_def int_of_real_def) + +lemma int_of_real_real[simp]: "int_of_real (real x) = x" +by (simp add: int_of_real_def) + +lemma real_int_of_real[simp]: "real_is_int x \ real (int_of_real x) = x" +by (auto simp add: int_of_real_def real_is_int_def) + +lemma real_is_int_add_int_of_real: "real_is_int a \ real_is_int b \ (int_of_real (a+b)) = (int_of_real a) + (int_of_real b)" +by (auto simp add: int_of_real_def real_is_int_def) + +lemma real_is_int_add[simp]: "real_is_int a \ real_is_int b \ real_is_int (a+b)" +apply (subst real_is_int_def2) +apply (simp add: real_is_int_add_int_of_real real_int_of_real) +done + +lemma int_of_real_sub: "real_is_int a \ real_is_int b \ (int_of_real (a-b)) = (int_of_real a) - (int_of_real b)" +by (auto simp add: int_of_real_def real_is_int_def) + +lemma real_is_int_sub[simp]: "real_is_int a \ real_is_int b \ real_is_int (a-b)" +apply (subst real_is_int_def2) +apply (simp add: int_of_real_sub real_int_of_real) +done + +lemma real_is_int_rep: "real_is_int x \ ?! (a::int). real a = x" +by (auto simp add: real_is_int_def) + +lemma int_of_real_mult: + assumes "real_is_int a" "real_is_int b" + shows "(int_of_real (a*b)) = (int_of_real a) * (int_of_real b)" +proof - + from prems have a: "?! (a'::int). real a' = a" by (rule_tac real_is_int_rep, auto) + from prems have b: "?! (b'::int). real b' = b" by (rule_tac real_is_int_rep, auto) + from a obtain a'::int where a':"a = real a'" by auto + from b obtain b'::int where b':"b = real b'" by auto + have r: "real a' * real b' = real (a' * b')" by auto + show ?thesis + apply (simp add: a' b') + apply (subst r) + apply (simp only: int_of_real_real) + done +qed + +lemma real_is_int_mult[simp]: "real_is_int a \ real_is_int b \ real_is_int (a*b)" +apply (subst real_is_int_def2) +apply (simp add: int_of_real_mult) +done + +lemma real_is_int_0[simp]: "real_is_int (0::real)" +by (simp add: real_is_int_def int_of_real_def) + +lemma real_is_int_1[simp]: "real_is_int (1::real)" +proof - + have "real_is_int (1::real) = real_is_int(real (1::int))" by auto + also have "\ = True" by (simp only: real_is_int_real) + ultimately show ?thesis by auto +qed + +lemma real_is_int_n1: "real_is_int (-1::real)" +proof - + have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto + also have "\ = True" by (simp only: real_is_int_real) + ultimately show ?thesis by auto +qed + +lemma real_is_int_number_of[simp]: "real_is_int ((number_of::bin\real) x)" +proof - + have neg1: "real_is_int (-1::real)" + proof - + have "real_is_int (-1::real) = real_is_int(real (-1::int))" by auto + also have "\ = True" by (simp only: real_is_int_real) + ultimately show ?thesis by auto + qed + + { + fix x::int + have "!! y. real_is_int ((number_of::bin\real) (Abs_Bin x))" + apply (simp add: number_of_eq) + apply (subst Abs_Bin_inverse) + apply (simp add: Bin_def) + apply (induct x) + apply (induct_tac n) + apply (simp) + apply (simp) + apply (induct_tac n) + apply (simp add: neg1) + proof - + fix n :: nat + assume rn: "(real_is_int (of_int (- (int (Suc n)))))" + have s: "-(int (Suc (Suc n))) = -1 + - (int (Suc n))" by simp + show "real_is_int (of_int (- (int (Suc (Suc n)))))" + apply (simp only: s of_int_add) + apply (rule real_is_int_add) + apply (simp add: neg1) + apply (simp only: rn) + done + qed + } + note Abs_Bin = this + { + fix x :: bin + have "? u. x = Abs_Bin u" + apply (rule exI[where x = "Rep_Bin x"]) + apply (simp add: Rep_Bin_inverse) + done + } + then obtain u::int where "x = Abs_Bin u" by auto + with Abs_Bin show ?thesis by auto +qed + +lemma int_of_real_0[simp]: "int_of_real (0::real) = (0::int)" +by (simp add: int_of_real_def) + +lemma int_of_real_1[simp]: "int_of_real (1::real) = (1::int)" +proof - + have 1: "(1::real) = real (1::int)" by auto + show ?thesis by (simp only: 1 int_of_real_real) +qed + +lemma int_of_real_number_of[simp]: "int_of_real (number_of b) = number_of b" +proof - + have "real_is_int (number_of b)" by simp + then have uu: "?! u::int. number_of b = real u" by (auto simp add: real_is_int_rep) + then obtain u::int where u:"number_of b = real u" by auto + have "number_of b = real ((number_of b)::int)" + by (simp add: number_of_eq real_of_int_def) + have ub: "number_of b = real ((number_of b)::int)" + by (simp add: number_of_eq real_of_int_def) + from uu u ub have unb: "u = number_of b" + by blast + have "int_of_real (number_of b) = u" by (simp add: u) + with unb show ?thesis by simp +qed + +lemma float_transfer_even: "even a \ float (a, b) = float (a div 2, b+1)" + apply (subst float_transfer[where a="a" and b="b" and c="-1", simplified]) + apply (simp_all add: pow2_def even_def real_is_int_def ring_eq_simps) + apply (auto) +proof - + fix q::int + have a:"b - (-1\int) = (1\int) + b" by arith + show "(float (q, (b - (-1\int)))) = (float (q, ((1\int) + b)))" + by (simp add: a) +qed + +consts + norm_float :: "int*int \ int*int" + +lemma int_div_zdiv: "int (a div b) = (int a) div (int b)" +apply (subst split_div, auto) +apply (subst split_zdiv, auto) +apply (rule_tac a="int (b * i) + int j" and b="int b" and r="int j" and r'=ja in IntDiv.unique_quotient) +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat) +done + +lemma int_mod_zmod: "int (a mod b) = (int a) mod (int b)" +apply (subst split_mod, auto) +apply (subst split_zmod, auto) +apply (rule_tac a="int (b * i) + int j" and b="int b" and q="int i" and q'=ia in IntDiv.unique_remainder) +apply (auto simp add: IntDiv.quorem_def int_eq_of_nat) +done + +lemma abs_div_2_less: "a \ 0 \ a \ -1 \ abs((a::int) div 2) < abs a" +by arith + +lemma terminating_norm_float: "\a. (a::int) \ 0 \ even a \ a \ 0 \ \a div 2\ < \a\" +apply (auto) +apply (rule abs_div_2_less) +apply (auto) +done + +ML {* simp_depth_limit := 2 *} +recdef norm_float "measure (% (a,b). nat (abs a))" + "norm_float (a,b) = (if (a \ 0) & (even a) then norm_float (a div 2, b+1) else (if a=0 then (0,0) else (a,b)))" +(hints simp: terminating_norm_float) +ML {* simp_depth_limit := 1000 *} + + +lemma norm_float: "float x = float (norm_float x)" +proof - + { + fix a b :: int + have norm_float_pair: "float (a,b) = float (norm_float (a,b))" + proof (induct a b rule: norm_float.induct) + case (1 u v) + show ?case + proof cases + assume u: "u \ 0 \ even u" + with prems have ind: "float (u div 2, v + 1) = float (norm_float (u div 2, v + 1))" by auto + with u have "float (u,v) = float (u div 2, v+1)" by (simp add: float_transfer_even) + then show ?thesis + apply (subst norm_float.simps) + apply (simp add: ind) + done + next + assume "~(u \ 0 \ even u)" + then show ?thesis + by (simp add: prems float_def) + qed + qed + } + note helper = this + have "? a b. x = (a,b)" by auto + then obtain a b where "x = (a, b)" by blast + then show ?thesis by (simp only: helper) +qed + +lemma pow2_int: "pow2 (int n) = 2^n" + by (simp add: pow2_def) + +lemma float_add: + "float (a1, e1) + float (a2, e2) = + (if e1<=e2 then float (a1+a2*2^(nat(e2-e1)), e1) + else float (a1*2^(nat (e1-e2))+a2, e2))" + apply (simp add: float_def ring_eq_simps) + apply (auto simp add: pow2_int[symmetric] pow2_add[symmetric]) + done + +lemma float_mult: + "float (a1, e1) * float (a2, e2) = + (float (a1 * a2, e1 + e2))" + by (simp add: float_def pow2_add) + +lemma float_minus: + "- (float (a,b)) = float (-a, b)" + by (simp add: float_def) + +lemma zero_less_pow2: + "0 < pow2 x" +proof - + { + fix y + have "0 <= y \ 0 < pow2 y" + by (induct y, induct_tac n, simp_all add: pow2_add) + } + note helper=this + show ?thesis + apply (case_tac "0 <= x") + apply (simp add: helper) + apply (subst pow2_neg) + apply (simp add: helper) + done +qed + +lemma zero_le_float: + "(0 <= float (a,b)) = (0 <= a)" + apply (auto simp add: float_def) + apply (auto simp add: zero_le_mult_iff zero_less_pow2) + apply (insert zero_less_pow2[of b]) + apply (simp_all) + done + +lemma float_abs: + "abs (float (a,b)) = (if 0 <= a then (float (a,b)) else (float (-a,b)))" + apply (auto simp add: abs_if) + apply (simp_all add: zero_le_float[symmetric, of a b] float_minus) + done + +lemma norm_0_1: "(0::_::number_ring) = Numeral0 & (1::_::number_ring) = Numeral1" + by auto + +lemma add_left_zero: "0 + a = (a::'a::comm_monoid_add)" + by simp + +lemma add_right_zero: "a + 0 = (a::'a::comm_monoid_add)" + by simp + +lemma mult_left_one: "1 * a = (a::'a::semiring_1)" + by simp + +lemma mult_right_one: "a * 1 = (a::'a::semiring_1)" + by simp + +lemma int_pow_0: "(a::int)^(Numeral0) = 1" + by simp + +lemma int_pow_1: "(a::int)^(Numeral1) = a" + by simp + +lemma zero_eq_Numeral0_nring: "(0::'a::number_ring) = Numeral0" + by simp + +lemma one_eq_Numeral1_nring: "(1::'a::number_ring) = Numeral1" + by simp + +lemma zero_eq_Numeral0_nat: "(0::nat) = Numeral0" + by simp + +lemma one_eq_Numeral1_nat: "(1::nat) = Numeral1" + by simp + +lemma zpower_Pls: "(z::int)^Numeral0 = Numeral1" + by simp + +lemma zpower_Min: "(z::int)^((-1)::nat) = Numeral1" +proof - + have 1:"((-1)::nat) = 0" + by simp + show ?thesis by (simp add: 1) +qed + +lemma fst_cong: "a=a' \ fst (a,b) = fst (a',b)" + by simp + +lemma snd_cong: "b=b' \ snd (a,b) = snd (a,b')" + by simp + +lemma lift_bool: "x \ x=True" + by simp + +lemma nlift_bool: "~x \ x=False" + by simp + +lemma not_false_eq_true: "(~ False) = True" by simp + +lemma not_true_eq_false: "(~ True) = False" by simp + + +lemmas binarith = + Pls_0_eq Min_1_eq + bin_pred_Pls bin_pred_Min bin_pred_1 bin_pred_0 + bin_succ_Pls bin_succ_Min bin_succ_1 bin_succ_0 + bin_add_Pls bin_add_Min bin_add_BIT_0 bin_add_BIT_10 + bin_add_BIT_11 bin_minus_Pls bin_minus_Min bin_minus_1 + bin_minus_0 bin_mult_Pls bin_mult_Min bin_mult_1 bin_mult_0 + bin_add_Pls_right bin_add_Min_right + +thm eq_number_of_eq + +lemma int_eq_number_of_eq: "(((number_of v)::int)=(number_of w)) = iszero ((number_of (bin_add v (bin_minus w)))::int)" + by simp + +lemma int_iszero_number_of_Pls: "iszero (Numeral0::int)" + by (simp only: iszero_number_of_Pls) + +lemma int_nonzero_number_of_Min: "~(iszero ((-1)::int))" + by simp +thm iszero_number_of_1 + +lemma int_iszero_number_of_0: "iszero ((number_of (w BIT False))::int) = iszero ((number_of w)::int)" + by simp + +lemma int_iszero_number_of_1: "\ iszero ((number_of (w BIT True))::int)" + by simp + +lemma int_less_number_of_eq_neg: "(((number_of x)::int) < number_of y) = neg ((number_of (bin_add x (bin_minus y)))::int)" + by simp + +lemma int_not_neg_number_of_Pls: "\ (neg (Numeral0::int))" + by simp + +lemma int_neg_number_of_Min: "neg (-1::int)" + by simp + +lemma int_neg_number_of_BIT: "neg ((number_of (w BIT x))::int) = neg ((number_of w)::int)" + by simp + +lemma int_le_number_of_eq: "(((number_of x)::int) \ number_of y) = (\ neg ((number_of (bin_add y (bin_minus x)))::int))" + by simp + +lemmas intarithrel = + (*abs_zero abs_one*) + int_eq_number_of_eq + lift_bool[OF int_iszero_number_of_Pls] nlift_bool[OF int_nonzero_number_of_Min] int_iszero_number_of_0 + 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] + int_neg_number_of_BIT int_le_number_of_eq + +lemma int_number_of_add_sym: "((number_of v)::int) + number_of w = number_of (bin_add v w)" + by simp + +lemma int_number_of_diff_sym: "((number_of v)::int) - number_of w = number_of (bin_add v (bin_minus w))" + by simp + +lemma int_number_of_mult_sym: "((number_of v)::int) * number_of w = number_of (bin_mult v w)" + by simp + +lemma int_number_of_minus_sym: "- ((number_of v)::int) = number_of (bin_minus v)" + by simp + +lemmas intarith = int_number_of_add_sym int_number_of_minus_sym int_number_of_diff_sym int_number_of_mult_sym + +lemmas natarith = (*zero_eq_Numeral0_nat one_eq_Numeral1_nat*) add_nat_number_of + diff_nat_number_of mult_nat_number_of eq_nat_number_of less_nat_number_of + +lemmas powerarith = nat_number_of zpower_number_of_even + zpower_number_of_odd[simplified zero_eq_Numeral0_nring one_eq_Numeral1_nring] + zpower_Pls zpower_Min + +lemmas floatarith[simplified norm_0_1] = float_add float_mult float_minus float_abs zero_le_float + +lemmas arith = binarith intarith intarithrel natarith powerarith floatarith not_false_eq_true not_true_eq_false + +(* needed for the verifying code generator *) +lemmas arith_no_let = arith[simplified Let_def] + +end + diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/FloatArith.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/FloatArith.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,56 @@ +structure FloatArith = +struct + +type float = IntInf.int * IntInf.int + +val izero = IntInf.fromInt 0 +fun imul a b = IntInf.* (a,b) +fun isub a b = IntInf.- (a,b) +fun iadd a b = IntInf.+ (a,b) + +val floatzero = (izero, izero) + +fun positive_part (a,b) = + (if IntInf.< (a,izero) then izero else a, b) + +fun negative_part (a,b) = + (if IntInf.< (a,izero) then a else izero, b) + +fun is_negative (a,b) = + if IntInf.< (a, izero) then true else false + +fun is_positive (a,b) = + if IntInf.< (izero, a) then true else false + +fun is_zero (a,b) = + if a = izero then true else false + +fun ipow2 a = IntInf.pow ((IntInf.fromInt 2), IntInf.toInt a) + +fun add (a1, b1) (a2, b2) = + if b1 < b2 then + (iadd a1 (imul a2 (ipow2 (isub b2 b1))), b1) + else + (iadd (imul a1 (ipow2 (isub b1 b2))) a2, b2) + +fun sub (a1, b1) (a2, b2) = + if b1 < b2 then + (isub a1 (imul a2 (ipow2 (isub b2 b1))), b1) + else + (isub (imul a1 (ipow2 (isub b1 b2))) a2, b2) + +fun neg (a, b) = (IntInf.~ a, b) + +fun is_equal a b = is_zero (sub a b) + +fun is_less a b = is_negative (sub a b) + +fun max a b = if is_less a b then b else a + +fun min a b = if is_less a b then a else b + +fun abs a = if is_negative a then neg a else a + +fun mul (a1, b1) (a2, b2) = (imul a1 a2, iadd b1 b2) + +end; diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/LinProg.thy --- a/src/HOL/Matrix/LinProg.thy Fri Sep 03 10:27:05 2004 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ -(* Title: HOL/Matrix/LinProg.thy - ID: $Id$ - Author: Steven Obua -*) - -theory LinProg = Matrix: - -lemma linprog_by_duality_approx_general: - assumes - fmuladdprops: - "! a b c d. a <= b & c <= d \ fadd a c <= fadd b d" - "! c a b. 0 <= c & a <= b \ fmul c a <= fmul c b" - "! a. fmul 0 a = 0" - "! a. fmul a 0 = 0" - "fadd 0 0 = 0" - "associative fadd" - "commutative fadd" - "associative fmul" - "distributive fmul fadd" - and specificprops: - "mult_matrix fmul fadd (combine_matrix fadd A dA) x <= (b::('a::{order, zero}) matrix)" - "mult_matrix fmul fadd y A = c" - "0 <= y" - shows - "combine_matrix fadd (mult_matrix fmul fadd c x) (mult_matrix fmul fadd (mult_matrix fmul fadd y dA) x) - <= mult_matrix fmul fadd y b" -proof - - let ?mul = "mult_matrix fmul fadd" - let ?add = "combine_matrix fadd" - let ?t1 = "?mul y (?mul (?add A dA) x)" - have a: "?t1 <= ?mul y b" by (rule le_left_mult, simp_all!) - have b: "?t1 = ?mul (?mul y (?add A dA)) x" by (simp! add: mult_matrix_assoc_simple[THEN sym]) - have assoc: "associative ?add" by (simp! add: combine_matrix_assoc) - have r_distr: "r_distributive ?mul ?add" - apply (rule r_distributive_matrix) - by (simp! add: distributive_def)+ - have l_distr: "l_distributive ?mul ?add" - apply (rule l_distributive_matrix) - by (simp! add: distributive_def)+ - have c:"?mul y (?add A dA) = ?add (?mul y A) (?mul y dA)" - by (insert r_distr, simp add: r_distributive_def) - have d:"?mul (?add (?mul y A) (?mul y dA)) x = ?add (?mul (?mul y A) x) (?mul (?mul y dA) x)" - by (insert l_distr, simp add: l_distributive_def) - have e:"\ = ?add (?mul c x) (?mul (?mul y dA) x)" by (simp!) - from a b c d e show "?add (?mul c x) (?mul (?mul y dA) x) <= ?mul y b" by simp -qed - -lemma linprog_by_duality_approx: - assumes - "(A + dA) * x <= (b::('a::lordered_ring) matrix)" - "y * A = c" - "0 <= y" - shows - "c * x + (y * dA) * x <= y * b" -apply (simp add: times_matrix_def plus_matrix_def) -apply (rule linprog_by_duality_approx_general) -apply (simp_all) -apply (simp_all add: associative_def add_assoc mult_assoc) -apply (simp_all add: commutative_def add_commute) -apply (simp_all add: distributive_def l_distributive_def r_distributive_def left_distrib right_distrib) -apply (simp_all! add: plus_matrix_def times_matrix_def) -apply (simp add: add_mono) -apply (simp add: mult_left_mono) -done - -lemma linprog_by_duality: - assumes - "A * x <= (b::('a::lordered_ring) matrix)" - "y * A = c" - "0 <= y" - shows - "c * x <= y * b" -proof - - have a:"(A + 0) * x <= b" by (simp!) - have b:"c * x + (y*0)*x <= y * b" by (rule linprog_by_duality_approx, simp_all!) - show "c * x <= y*b" by (insert b, simp) -qed - -end diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/Matrix.thy --- a/src/HOL/Matrix/Matrix.thy Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/Matrix/Matrix.thy Fri Sep 03 17:10:36 2004 +0200 @@ -21,7 +21,10 @@ times_matrix_def: "A * B == mult_matrix (op *) (op +) A B" lemma is_meet_combine_matrix_meet: "is_meet (combine_matrix meet)" -by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le) + by (simp_all add: is_meet_def le_matrix_def meet_left_le meet_right_le meet_imp_le) + +lemma is_join_combine_matrix_join: "is_join (combine_matrix join)" + by (simp_all add: is_join_def le_matrix_def join_left_le join_right_le join_imp_le) instance matrix :: (lordered_ab_group) lordered_ab_group_meet proof @@ -284,7 +287,20 @@ apply (auto) done +lemma Rep_minus[simp]: "Rep_matrix (-(A::_::lordered_ab_group)) x y = - (Rep_matrix A x y)" + by (simp add: minus_matrix_def) +lemma join_matrix: "join (A::('a::lordered_ring) matrix) B = combine_matrix join A B" + apply (insert join_unique[of "(combine_matrix join)::('a matrix \ 'a matrix \ 'a matrix)", simplified is_join_combine_matrix_join]) + apply (simp) + done +lemma meet_matrix: "meet (A::('a::lordered_ring) matrix) B = combine_matrix meet A B" + apply (insert meet_unique[of "(combine_matrix meet)::('a matrix \ 'a matrix \ 'a matrix)", simplified is_meet_combine_matrix_meet]) + apply (simp) + done + +lemma Rep_abs[simp]: "Rep_matrix (abs (A::_::lordered_ring)) x y = abs (Rep_matrix A x y)" + by (simp add: abs_lattice join_matrix) end diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/MatrixGeneral.thy --- a/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/Matrix/MatrixGeneral.thy Fri Sep 03 17:10:36 2004 +0200 @@ -1392,7 +1392,44 @@ apply (rule le_foldseq) by (auto) +lemma spec2: "! j i. P j i \ P j i" by blast +lemma neg_imp: "(\ Q \ \ P) \ P \ Q" by blast + lemma singleton_matrix_le[simp]: "(singleton_matrix j i a <= singleton_matrix j i b) = (a <= (b::_::order))" by (auto simp add: le_matrix_def) +lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))" + apply (auto) + apply (simp add: le_matrix_def) + apply (drule_tac j=j and i=i in spec2) + apply (simp) + apply (simp add: le_matrix_def) + done + +lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)" + apply (auto) + apply (simp add: le_matrix_def) + apply (drule_tac j=j and i=i in spec2) + apply (simp) + apply (simp add: le_matrix_def) + done + +lemma move_matrix_le_zero[simp]: "0 <= j \ 0 <= i \ (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))" + apply (auto simp add: le_matrix_def neg_def) + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) + apply (auto) + done + +lemma move_matrix_zero_le[simp]: "0 <= j \ 0 <= i \ (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)" + apply (auto simp add: le_matrix_def neg_def) + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) + apply (auto) + done + +lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \ 0 <= i \ (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))" + apply (auto simp add: le_matrix_def neg_def) + apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) + apply (auto) + done + end diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/MatrixLP.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/MatrixLP.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,83 @@ +use "MatrixLP_gensimp.ML"; + +signature MATRIX_LP = +sig + val lp_dual_estimate : string -> int -> thm + val simplify : thm -> thm +end + +structure MatrixLP : MATRIX_LP = +struct + +val _ = SimprocsCodegen.use_simprocs_code sg geninput + +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' +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''''' +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''''' +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''''' +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''''' +val simp_sorted_sparse_matrix = simp_SparseMatrix_sorted_sparse_matrix +val simp_sorted_spvec = simp_SparseMatrix_sorted_spvec + +fun single_arg simp ct = snd (simp (snd (Thm.dest_comb ct))) + +fun double_arg simp ct = + let + val (a, y) = Thm.dest_comb ct + val (_, x) = Thm.dest_comb a + in + snd (simp x y) + end + +fun spmc ct = + (case term_of ct of + ((Const ("SparseMatrix.le_spmat", _)) $ _) => + single_arg simp_le_spmat ct + | ((Const ("SparseMatrix.add_spmat", _)) $ _) => + single_arg simp_add_spmat ct + | ((Const ("SparseMatrix.diff_spmat", _)) $ _ $ _) => + double_arg simp_diff_spmat ct + | ((Const ("SparseMatrix.abs_spmat", _)) $ _) => + single_arg simp_abs_spmat ct + | ((Const ("SparseMatrix.mult_spmat", _)) $ _ $ _) => + double_arg simp_mult_spmat ct + | ((Const ("SparseMatrix.sorted_sparse_matrix", _)) $ _) => + single_arg simp_sorted_sparse_matrix ct + | ((Const ("SparseMatrix.sorted_spvec", ty)) $ _) => + single_arg simp_sorted_spvec ct + | _ => raise Fail_conv) + +fun lp_dual_estimate lptfile prec = + let + val th = inst_real (thm "SparseMatrix.spm_linprog_dual_estimate_1") + val sg = sign_of (theory "MatrixLP") + val (y, (A1, A2), (c1, c2), b, r) = + let + open fspmlp + val l = load lptfile prec false + in + (y l, A l, c l, b l, r l) + end + fun var s x = (cterm_of sg (Var ((s,0), FloatSparseMatrixBuilder.real_spmatT)), x) + 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 + in + th + end + +fun simplify th = + let + val simp_th = botc spmc (cprop_of th) + val th = strip_shyps (equal_elim simp_th th) + fun removeTrue th = removeTrue (implies_elim th TrueI) handle _ => th + in + removeTrue th + end + +fun float2real (x,y) = (Real.fromInt x) * (Math.pow (2.0, Real.fromInt y)) + +end + +val result = ref ([]:thm list) + +fun run c = timeit (fn () => (result := [MatrixLP.simplify c]; ())) + diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/MatrixLP.thy --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/MatrixLP.thy Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,4 @@ +theory MatrixLP = Float + SparseMatrix: + +end + diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/MatrixLP_gensimp.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/MatrixLP_gensimp.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,35 @@ +open BasisLibrary; + +use "Cplex.ML"; +use "CplexMatrixConverter.ML"; +use "ExactFloatingPoint.ML"; +use "FloatArith.ML"; +use "FloatSparseMatrixBuilder.ML"; +use "fspmlp.ML"; +use "codegen_prep.ML"; +use "eq_codegen.ML"; +use "conv.ML"; + + +val th = theory "MatrixLP" +val sg = sign_of th + +fun prep x = standard (mk_meta_eq x) + +fun inst_real thm = standard (Thm.instantiate ([(fst (hd (term_tvars (prop_of thm))), + ctyp_of sg HOLogic.realT)], []) thm); + +val spm_lp_dual_estimate_simp = + (thms "Float.arith_no_let") @ + (map inst_real (thms "SparseMatrix.sparse_row_matrix_arith_simps")) @ + (thms "SparseMatrix.sorted_sp_simps") @ + [thm "Product_Type.fst_conv", thm "Product_Type.snd_conv"] @ + (thms "SparseMatrix.boolarith") + +val simp_with = Simplifier.simplify (HOL_basic_ss addsimps [thm "SparseMatrix.if_case_eq", thm "Float.norm_0_1"]) + +val spm_lp_dual_estimate_simp = map simp_with spm_lp_dual_estimate_simp + +val geninput = codegen_prep.prepare_thms spm_lp_dual_estimate_simp + +val _ = SimprocsCodegen.use_simprocs_code sg geninput diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/ROOT.ML --- a/src/HOL/Matrix/ROOT.ML Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/Matrix/ROOT.ML Fri Sep 03 17:10:36 2004 +0200 @@ -1,10 +1,1 @@ -(* Title: HOL/Matrix/ROOT.ML - ID: $Id$ - Author: Steven Obua - License: 2004 Technische Universität München - -Theory of matrices with an application of matrix theory to linear -programming. -*) - -use_thy "LinProg"; +use_thy "MatrixLP" diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/SparseMatrix.thy --- a/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/Matrix/SparseMatrix.thy Fri Sep 03 17:10:36 2004 +0200 @@ -94,9 +94,62 @@ done consts + abs_spvec :: "('a::lordered_ring) spvec \ 'a spvec" + minus_spvec :: "('a::lordered_ring) spvec \ 'a spvec" smult_spvec :: "('a::lordered_ring) \ 'a spvec \ 'a spvec" addmult_spvec :: "('a::lordered_ring) * 'a spvec * 'a spvec \ 'a spvec" +primrec + "minus_spvec [] = []" + "minus_spvec (a#as) = (fst a, -(snd a))#(minus_spvec as)" + +primrec + "abs_spvec [] = []" + "abs_spvec (a#as) = (fst a, abs (snd a))#(abs_spvec as)" + +lemma sparse_row_vector_minus: + "sparse_row_vector (minus_spvec v) = - (sparse_row_vector v)" + apply (induct v) + apply (simp_all add: sparse_row_vector_cons) + apply (simp add: Rep_matrix_inject[symmetric]) + apply (rule ext)+ + apply simp + done + +lemma sparse_row_vector_abs: + "sorted_spvec v \ sparse_row_vector (abs_spvec v) = abs (sparse_row_vector v)" + apply (induct v) + apply (simp_all add: sparse_row_vector_cons) + apply (frule_tac sorted_spvec_cons1, simp) + apply (simp only: Rep_matrix_inject[symmetric]) + apply (rule ext)+ + apply auto + apply (subgoal_tac "Rep_matrix (sparse_row_vector list) 0 a = 0") + apply (simp) + apply (rule sorted_sparse_row_vector_zero) + apply auto + done + +lemma sorted_spvec_minus_spvec: + "sorted_spvec v \ sorted_spvec (minus_spvec v)" + apply (induct v) + apply (simp) + apply (frule sorted_spvec_cons1, simp) + apply (simp add: sorted_spvec.simps) + apply (case_tac list) + apply (simp_all) + done + +lemma sorted_spvec_abs_spvec: + "sorted_spvec v \ sorted_spvec (abs_spvec v)" + apply (induct v) + apply (simp) + apply (frule sorted_spvec_cons1, simp) + apply (simp add: sorted_spvec.simps) + apply (case_tac list) + apply (simp_all) + done + defs smult_spvec_def: "smult_spvec y arr == map (% a. (fst a, y * snd a)) arr" @@ -542,9 +595,6 @@ else (le_spvec(snd a, snd b) & le_spmat (as, bs))))" -lemma spec2: "! j i. P j i \ P j i" by blast -lemma neg_imp: "(\ Q \ \ P) \ P \ Q" by blast - constdefs disj_matrices :: "('a::zero) matrix \ 'a matrix \ bool" "disj_matrices A B == (! j i. (Rep_matrix A j i \ 0) \ (Rep_matrix B j i = 0)) & (! j i. (Rep_matrix B j i \ 0) \ (Rep_matrix A j i = 0))" @@ -597,22 +647,6 @@ (B + A <= C) = (A <= C & (B::('a::lordered_ab_group) matrix) <= 0)" by (auto simp add: disj_matrices_add[of B A 0 C,simplified] disj_matrices_commute) -lemma singleton_le_zero[simp]: "(singleton_matrix j i x <= 0) = (x <= (0::'a::{order,zero}))" - apply (auto) - apply (simp add: le_matrix_def) - apply (drule_tac j=j and i=i in spec2) - apply (simp) - apply (simp add: le_matrix_def) - done - -lemma singleton_ge_zero[simp]: "(0 <= singleton_matrix j i x) = ((0::'a::{order,zero}) <= x)" - apply (auto) - apply (simp add: le_matrix_def) - apply (drule_tac j=j and i=i in spec2) - apply (simp) - apply (simp add: le_matrix_def) - done - lemma disj_sparse_row_singleton: "i <= j \ sorted_spvec((j,y)#v) \ disj_matrices (sparse_row_vector v) (singleton_matrix 0 i x)" apply (simp add: disj_matrices_def) apply (rule conjI) @@ -641,27 +675,6 @@ lemma disj_singleton_matrices[simp]: "disj_matrices (singleton_matrix j i x) (singleton_matrix u v y) = (j \ u | i \ v | x = 0 | y = 0)" by (auto simp add: disj_matrices_def) -lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \ (sorted_spvec b) \ (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)" - apply (rule le_spvec.induct) - 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) - apply (rule conjI, intro strip) - apply (simp add: sorted_spvec_cons1) - apply (subst disj_matrices_add_x_le) - apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute) - apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute) - apply (simp, blast) - apply (intro strip, rule conjI, intro strip) - apply (simp add: sorted_spvec_cons1) - apply (subst disj_matrices_add_le_x) - 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) - apply (blast) - apply (intro strip) - apply (simp add: sorted_spvec_cons1) - apply (case_tac "a=aa", simp_all) - apply (subst disj_matrices_add) - apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute) - done - lemma disj_move_sparse_vec_mat[simplified disj_matrices_commute]: "j <= a \ sorted_spvec((a,c)#as) \ disj_matrices (move_matrix (sparse_row_vector b) (int j) i) (sparse_row_matrix as)" apply (auto simp add: neg_def disj_matrices_def) @@ -685,24 +698,28 @@ apply (rule nrows, rule order_trans[of _ 1], simp, drule nrows_notzero, drule less_le_trans[OF _ nrows_spvec], arith)+ done -lemma move_matrix_le_zero[simp]: "0 <= j \ 0 <= i \ (move_matrix A j i <= 0) = (A <= (0::('a::{order,zero}) matrix))" - apply (auto simp add: le_matrix_def neg_def) - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) - apply (auto) +lemma le_spvec_iff_sparse_row_le[rule_format]: "(sorted_spvec a) \ (sorted_spvec b) \ (le_spvec (a,b)) = (sparse_row_vector a <= sparse_row_vector b)" + apply (rule le_spvec.induct) + 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) + apply (rule conjI, intro strip) + apply (simp add: sorted_spvec_cons1) + apply (subst disj_matrices_add_x_le) + apply (simp add: disj_sparse_row_singleton[OF less_imp_le] disj_matrices_x_add disj_matrices_commute) + apply (simp add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute) + apply (simp, blast) + apply (intro strip, rule conjI, intro strip) + apply (simp add: sorted_spvec_cons1) + apply (subst disj_matrices_add_le_x) + 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) + apply (blast) + apply (intro strip) + apply (simp add: sorted_spvec_cons1) + apply (case_tac "a=aa", simp_all) + apply (subst disj_matrices_add) + apply (simp_all add: disj_sparse_row_singleton[OF order_refl] disj_matrices_commute) done -lemma move_matrix_zero_le[simp]: "0 <= j \ 0 <= i \ (0 <= move_matrix A j i) = ((0::('a::{order,zero}) matrix) <= A)" - apply (auto simp add: le_matrix_def neg_def) - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) - apply (auto) - done - -lemma move_matrix_le_move_matrix_iff[simp]: "0 <= j \ 0 <= i \ (move_matrix A j i <= move_matrix B j i) = (A <= (B::('a::{order,zero}) matrix))" - apply (auto simp add: le_matrix_def neg_def) - apply (drule_tac j="ja+(nat j)" and i="ia+(nat i)" in spec2) - apply (auto) - done - lemma le_spvec_empty2_sparse_row[rule_format]: "(sorted_spvec b) \ (le_spvec (b,[]) = (sparse_row_vector b <= 0))" apply (induct b) apply (simp_all add: sorted_spvec_cons1) @@ -752,39 +769,169 @@ apply (simp add: sorted_spvec_cons1 le_spvec_iff_sparse_row_le) done -term smult_spvec -term addmult_spvec -term add_spvec -term mult_spvec_spmat -term mult_spmat -term add_spmat -term le_spvec -term le_spmat -term sorted_spvec -term sorted_spmat +ML {* simp_depth_limit := 999 *} + +consts + abs_spmat :: "('a::lordered_ring) spmat \ 'a spmat" + minus_spmat :: "('a::lordered_ring) spmat \ 'a spmat" + +primrec + "abs_spmat [] = []" + "abs_spmat (a#as) = (fst a, abs_spvec (snd a))#(abs_spmat as)" + +primrec + "minus_spmat [] = []" + "minus_spmat (a#as) = (fst a, minus_spvec (snd a))#(minus_spmat as)" + +lemma sparse_row_matrix_minus: + "sparse_row_matrix (minus_spmat A) = - (sparse_row_matrix A)" + apply (induct A) + apply (simp_all add: sparse_row_vector_minus sparse_row_matrix_cons) + apply (subst Rep_matrix_inject[symmetric]) + apply (rule ext)+ + apply simp + done -thm sparse_row_mult_spmat -thm sparse_row_add_spmat -thm le_spmat_iff_sparse_row_le +lemma Rep_sparse_row_vector_zero: "x \ 0 \ Rep_matrix (sparse_row_vector v) x y = 0" +proof - + assume x:"x \ 0" + have r:"nrows (sparse_row_vector v) <= Suc 0" by (rule nrows_spvec) + show ?thesis + apply (rule nrows) + apply (subgoal_tac "Suc 0 <= x") + apply (insert r) + apply (simp only:) + apply (insert x) + apply arith + done +qed + +lemma sparse_row_matrix_abs: + "sorted_spvec A \ sorted_spmat A \ sparse_row_matrix (abs_spmat A) = abs (sparse_row_matrix A)" + apply (induct A) + apply (simp_all add: sparse_row_vector_abs sparse_row_matrix_cons) + apply (frule_tac sorted_spvec_cons1, simp) + apply (subst Rep_matrix_inject[symmetric]) + apply (rule ext)+ + apply auto + apply (case_tac "x=a") + apply (simp) + apply (subst sorted_sparse_row_matrix_zero) + apply auto + apply (subst Rep_sparse_row_vector_zero) + apply (simp_all add: neg_def) + done + +lemma sorted_spvec_minus_spmat: "sorted_spvec A \ sorted_spvec (minus_spmat A)" + apply (induct A) + apply (simp) + apply (frule sorted_spvec_cons1, simp) + apply (simp add: sorted_spvec.simps) + apply (case_tac list) + apply (simp_all) + done + +lemma sorted_spvec_abs_spmat: "sorted_spvec A \ sorted_spvec (abs_spmat A)" + apply (induct A) + apply (simp) + apply (frule sorted_spvec_cons1, simp) + apply (simp add: sorted_spvec.simps) + apply (case_tac list) + apply (simp_all) + done + +lemma sorted_spmat_minus_spmat: "sorted_spmat A \ sorted_spmat (minus_spmat A)" + apply (induct A) + apply (simp_all add: sorted_spvec_minus_spvec) + done + +lemma sorted_spmat_abs_spmat: "sorted_spmat A \ sorted_spmat (abs_spmat A)" + apply (induct A) + apply (simp_all add: sorted_spvec_abs_spvec) + done -thm sorted_spvec_mult_spmat -thm sorted_spmat_mult_spmat -thm sorted_spvec_add_spmat -thm sorted_spmat_add_spmat +constdefs + diff_spmat :: "('a::lordered_ring) spmat \ 'a spmat \ 'a spmat" + "diff_spmat A B == add_spmat (A, minus_spmat B)" + +lemma sorted_spmat_diff_spmat: "sorted_spmat A \ sorted_spmat B \ sorted_spmat (diff_spmat A B)" + by (simp add: diff_spmat_def sorted_spmat_minus_spmat sorted_spmat_add_spmat) + +lemma sorted_spvec_diff_spmat: "sorted_spvec A \ sorted_spvec B \ sorted_spvec (diff_spmat A B)" + by (simp add: diff_spmat_def sorted_spvec_minus_spmat sorted_spvec_add_spmat) + +lemma sparse_row_diff_spmat: "sparse_row_matrix (diff_spmat A B ) = (sparse_row_matrix A) - (sparse_row_matrix B)" + by (simp add: diff_spmat_def sparse_row_add_spmat sparse_row_matrix_minus) + +constdefs + sorted_sparse_matrix :: "'a spmat \ bool" + "sorted_sparse_matrix A == (sorted_spvec A) & (sorted_spmat A)" + +lemma sorted_sparse_matrix_imp_spvec: "sorted_sparse_matrix A \ sorted_spvec A" + by (simp add: sorted_sparse_matrix_def) + +lemma sorted_sparse_matrix_imp_spmat: "sorted_sparse_matrix A \ sorted_spmat A" + by (simp add: sorted_sparse_matrix_def) + +lemmas sparse_row_matrix_op_simps = + sorted_sparse_matrix_imp_spmat sorted_sparse_matrix_imp_spvec + sparse_row_add_spmat sorted_spvec_add_spmat sorted_spmat_add_spmat + sparse_row_diff_spmat sorted_spvec_diff_spmat sorted_spmat_diff_spmat + sparse_row_matrix_minus sorted_spvec_minus_spmat sorted_spmat_minus_spmat + sparse_row_mult_spmat sorted_spvec_mult_spmat sorted_spmat_mult_spmat + sparse_row_matrix_abs sorted_spvec_abs_spmat sorted_spmat_abs_spmat + le_spmat_iff_sparse_row_le + +lemma zero_eq_Numeral0: "(0::_::number_ring) = Numeral0" by simp -thm smult_spvec_empty -thm smult_spvec_cons -thm addmult_spvec.simps -thm add_spvec.simps -thm add_spmat.simps -thm mult_spvec_spmat.simps -thm mult_spmat.simps -thm le_spvec.simps -thm le_spmat.simps -thm sorted_spvec.simps -thm sorted_spmat.simps +lemmas sparse_row_matrix_arith_simps[simplified zero_eq_Numeral0] = + mult_spmat.simps mult_spvec_spmat.simps + addmult_spvec.simps + smult_spvec_empty smult_spvec_cons + add_spmat.simps add_spvec.simps + minus_spmat.simps minus_spvec.simps + abs_spmat.simps abs_spvec.simps + diff_spmat_def + le_spmat.simps le_spvec.simps + +lemmas sorted_sp_simps = + sorted_spvec.simps + sorted_spmat.simps + sorted_sparse_matrix_def + +lemma bool1: "(\ True) = False" by blast +lemma bool2: "(\ False) = True" by blast +lemma bool3: "((P\bool) \ True) = P" by blast +lemma bool4: "(True \ (P\bool)) = P" by blast +lemma bool5: "((P\bool) \ False) = False" by blast +lemma bool6: "(False \ (P\bool)) = False" by blast +lemma bool7: "((P\bool) \ True) = True" by blast +lemma bool8: "(True \ (P\bool)) = True" by blast +lemma bool9: "((P\bool) \ False) = P" by blast +lemma bool10: "(False \ (P\bool)) = P" by blast +lemmas boolarith = bool1 bool2 bool3 bool4 bool5 bool6 bool7 bool8 bool9 bool10 + +lemma if_case_eq: "(if b then x else y) = (case b of True => x | False => y)" by simp + +lemma spm_linprog_dual_estimate_1: + assumes + "sorted_sparse_matrix A1" + "sorted_sparse_matrix A2" + "sorted_sparse_matrix c1" + "sorted_sparse_matrix c2" + "sorted_sparse_matrix y" + "sorted_spvec b" + "sorted_spvec r" + "le_spmat ([], y)" + "A * x \ sparse_row_matrix (b::('a::lordered_ring) spmat)" + "sparse_row_matrix A1 <= A" + "A <= sparse_row_matrix A2" + "sparse_row_matrix c1 <= c" + "c <= sparse_row_matrix c2" + "abs x \ sparse_row_matrix r" + shows + "c * x \ sparse_row_matrix (add_spmat (mult_spmat y b, mult_spmat (add_spmat (add_spmat (mult_spmat y (diff_spmat A2 A1), + abs_spmat (diff_spmat (mult_spmat y A1) c1)), diff_spmat c2 c1)) r))" + by (insert prems, simp add: sparse_row_matrix_op_simps linprog_dual_estimate_1[where A=A]) end - - - diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/codegen_prep.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/codegen_prep.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,110 @@ +structure codegen_prep = +struct + +exception Prepare_thms of string; + +fun is_meta_eq th = + let + fun check ((Const ("==", _)) $ _ $ _) = true + | check _ = false + in + check (concl_of th) + end + +fun printty ty = + let + fun immerse s [] = [] + | immerse s [x] = [x] + | immerse s (x::xs) = x::s::(immerse s xs) + + fun py t = + let + val (name, args) = dest_Type t + val args' = map printty args + in + concat (immerse "_" (name::args')) + end + + val (args, res) = strip_type ty + val tystrs = map py (args@[res]) + + val s = "'"^(concat (immerse "_" tystrs))^"'" + in + s + end + +fun head_name th = + let val s = + let + val h = fst (strip_comb (hd (snd (strip_comb (concl_of th))))) + val n = fst (dest_Const h) + val ty = snd (dest_Const h) + val hn = fst (dest_Const h) + in + hn^"_"^(printty ty) handle _ => (writeln ("warning: polymorphic "^hn); hn) + end + in + s + end + +val mangle_id = + let + fun tr #"=" = "_eq_" + | tr #"." = "_" + | tr #" " = "_" + | tr #"<" = "_le_" + | tr #">" = "_ge_" + | tr #"(" = "_bro_" + | tr #")" = "_brc_" + | tr #"+" = "_plus_" + | tr #"*" = "_mult_" + | tr #"-" = "_minus_" + | tr #"|" = "_or_" + | tr #"&" = "_and_" + | tr x = Char.toString x + fun m s = "simp_"^(String.translate tr s) + in + m + end + +fun group f [] = [] + | group f (x::xs) = + let + val g = group f xs + val key = f x + in + case assoc (g, key) of + None => (key, [x])::g + | Some l => overwrite (g, (key, x::l)) + end + +fun prepare_thms ths = + let + val ths = (filter is_meta_eq ths) @ (map (standard o mk_meta_eq) (filter (not o is_meta_eq) ths)) + val _ = if forall Thm.no_prems ths then () else raise (Prepare_thms "premisse found") + val thmgroups = group head_name ths + val test_clashgroups = group (fn (a,b) => mangle_id a) thmgroups + val _ = if (length thmgroups <> length test_clashgroups) then raise (Prepare_thms "clash after name mangling") else () + + fun prep (name, ths) = + let + val m = mangle_id name + + in + (m, true, ths) + end + + val thmgroups = map prep thmgroups + in + (thmgroups) + end + +fun writestr filename s = + let + val f = TextIO.openOut filename + val _ = TextIO.output(f, s) + val _ = TextIO.closeOut f + in + () + end +end diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/conv.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/conv.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,26 @@ +exception Fail_conv; + +fun orelsec conv1 conv2 ct = conv1 ct handle Fail_conv => conv2 ct + +val allc = Thm.reflexive + +fun thenc conv1 conv2 ct = + let + fun rhs_of t = snd (Thm.dest_comb (strip_imp_concl (cprop_of t))) + + val eq = conv1 ct + in + Thm.transitive eq (conv2 (rhs_of eq)) + end + +fun subc conv ct = + case term_of ct of + _ $ _ => + let + val (ct1, ct2) = Thm.dest_comb ct + in + Thm.combination (conv ct1) (conv ct2) + end + | _ => allc ct + +fun botc conv ct = thenc (subc (botc conv)) (orelsec conv allc) ct \ No newline at end of file diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/eq_codegen.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/eq_codegen.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,493 @@ +fun inst_cterm inst ct = fst (Drule.dest_equals + (Thm.cprop_of (Thm.instantiate inst (reflexive ct)))); +fun tyinst_cterm tyinst = inst_cterm (tyinst, []); + +val bla = ref ([] : term list); + +(******************************************************) +(* Code generator for equational proofs *) +(******************************************************) +fun my_mk_meta_eq thm = + let + val (_, eq) = Thm.dest_comb (cprop_of thm); + val (ct, rhs) = Thm.dest_comb eq; + val (_, lhs) = Thm.dest_comb ct + in Thm.implies_elim (Drule.instantiate' [Some (ctyp_of_term lhs)] + [Some lhs, Some rhs] eq_reflection) thm + end; + +structure SimprocsCodegen = +struct + +val simp_thms = ref ([] : thm list); + +fun parens b = if b then Pretty.enclose "(" ")" else Pretty.block; + +fun gen_mk_val f xs ps = Pretty.block ([Pretty.str "val ", + f (length xs > 1) (flat + (separate [Pretty.str ",", Pretty.brk 1] (map (single o Pretty.str) xs))), + Pretty.str " =", Pretty.brk 1] @ ps @ [Pretty.str ";"]); + +val mk_val = gen_mk_val parens; +val mk_vall = gen_mk_val (K (Pretty.enclose "[" "]")); + +fun rename s = if s mem ThmDatabase.ml_reserved then s ^ "'" else s; + +fun mk_decomp_name (Var ((s, i), _)) = rename (if i=0 then s else s ^ string_of_int i) + | mk_decomp_name (Const (s, _)) = rename (Codegen.mk_id (Sign.base_name s)) + | mk_decomp_name _ = "ct"; + +fun decomp_term_code cn ((vs, bs, ps), (v, t)) = + if exists (equal t o fst) bs then (vs, bs, ps) + else (case t of + Var _ => (vs, bs @ [(t, v)], ps) + | Const _ => (vs, if cn then bs @ [(t, v)] else bs, ps) + | Bound _ => (vs, bs, ps) + | Abs (s, T, t) => + let + val v1 = variant vs s; + val v2 = variant (v1 :: vs) (mk_decomp_name t) + in + decomp_term_code cn ((v1 :: v2 :: vs, + bs @ [(Free (s, T), v1)], + ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_abs", Pretty.brk 1, + Pretty.str "None", Pretty.brk 1, Pretty.str v]]), (v2, t)) + end + | t $ u => + let + val v1 = variant vs (mk_decomp_name t); + val v2 = variant (v1 :: vs) (mk_decomp_name u); + val (vs', bs', ps') = decomp_term_code cn ((v1 :: v2 :: vs, bs, + ps @ [mk_val [v1, v2] [Pretty.str "Thm.dest_comb", Pretty.brk 1, + Pretty.str v]]), (v1, t)); + val (vs'', bs'', ps'') = decomp_term_code cn ((vs', bs', ps'), (v2, u)) + in + if bs'' = bs then (vs, bs, ps) else (vs'', bs'', ps'') + end); + +val strip_tv = implode o tl o explode; + +fun mk_decomp_tname (TVar ((s, i), _)) = + strip_tv ((if i=0 then s else s ^ string_of_int i) ^ "T") + | mk_decomp_tname (Type (s, _)) = Codegen.mk_id (Sign.base_name s) ^ "T" + | mk_decomp_tname _ = "cT"; + +fun decomp_type_code ((vs, bs, ps), (v, TVar (ixn, _))) = + if exists (equal ixn o fst) bs then (vs, bs, ps) + else (vs, bs @ [(ixn, v)], ps) + | decomp_type_code ((vs, bs, ps), (v, Type (_, Ts))) = + let + val vs' = variantlist (map mk_decomp_tname Ts, vs); + val (vs'', bs', ps') = + foldl decomp_type_code ((vs @ vs', bs, ps @ + [mk_vall vs' [Pretty.str "Thm.dest_ctyp", Pretty.brk 1, + Pretty.str v]]), vs' ~~ Ts) + in + if bs' = bs then (vs, bs, ps) else (vs'', bs', ps') + end; + +fun gen_mk_bindings s dest decomp ((vs, bs, ps), (v, x)) = + let + val s' = variant vs s; + val (vs', bs', ps') = decomp ((s' :: vs, bs, ps @ + [mk_val [s'] (dest v)]), (s', x)) + in + if bs' = bs then (vs, bs, ps) else (vs', bs', ps') + end; + +val mk_term_bindings = gen_mk_bindings "ct" + (fn s => [Pretty.str "cprop_of", Pretty.brk 1, Pretty.str s]) + (decomp_term_code true); + +val mk_type_bindings = gen_mk_bindings "cT" + (fn s => [Pretty.str "Thm.ctyp_of_term", Pretty.brk 1, Pretty.str s]) + decomp_type_code; + +fun pretty_pattern b (Const (s, _)) = Pretty.block [Pretty.str "Const", + Pretty.brk 1, Pretty.str ("(\"" ^ s ^ "\", _)")] + | pretty_pattern b (t as _ $ _) = parens b + (flat (separate [Pretty.str " $", Pretty.brk 1] + (map (single o pretty_pattern true) (op :: (strip_comb t))))) + | pretty_pattern b _ = Pretty.str "_"; + +fun term_consts' t = foldl_aterms + (fn (cs, c as Const _) => c ins cs | (cs, _) => cs) ([], t); + +fun mk_apps s b p [] = p + | mk_apps s b p (q :: qs) = + mk_apps s b (parens (b orelse not (null qs)) + [Pretty.str s, Pretty.brk 1, p, Pretty.brk 1, q]) qs; + +fun mk_refleq eq ct = mk_val [eq] [Pretty.str ("Thm.reflexive " ^ ct)]; + +fun mk_tyinst ((s, i), s') = + Pretty.block [Pretty.str ("((" ^ quote s ^ ","), Pretty.brk 1, + Pretty.str (string_of_int i ^ "),"), Pretty.brk 1, + Pretty.str (s' ^ ")")]; + +fun inst_ty b ty_bs t s = (case term_tvars t of + [] => Pretty.str s + | Ts => parens b [Pretty.str "tyinst_cterm", Pretty.brk 1, + Pretty.list "[" "]" (map (fn (ixn, _) => mk_tyinst + (ixn, the (assoc (ty_bs, ixn)))) Ts), + Pretty.brk 1, Pretty.str s]); + +fun mk_cterm_code b ty_bs ts xs (vals, t $ u) = + let + val (vals', p1) = mk_cterm_code true ty_bs ts xs (vals, t); + val (vals'', p2) = mk_cterm_code true ty_bs ts xs (vals', u) + in + (vals'', parens b [Pretty.str "Thm.capply", Pretty.brk 1, + p1, Pretty.brk 1, p2]) + end + | mk_cterm_code b ty_bs ts xs (vals, Abs (s, T, t)) = + let + val u = Free (s, T); + val Some s' = assoc (ts, u); + val p = Pretty.str s'; + val (vals', p') = mk_cterm_code true ty_bs ts (p :: xs) + (if null (typ_tvars T) then vals + else vals @ [(u, (("", s'), [mk_val [s'] [inst_ty true ty_bs u s']]))], t) + in (vals', + parens b [Pretty.str "Thm.cabs", Pretty.brk 1, p, Pretty.brk 1, p']) + end + | mk_cterm_code b ty_bs ts xs (vals, Bound i) = (vals, nth_elem (i, xs)) + | mk_cterm_code b ty_bs ts xs (vals, t) = (case assoc (vals, t) of + None => + let val Some s = assoc (ts, t) + in (if is_Const t andalso not (null (term_tvars t)) then + vals @ [(t, (("", s), [mk_val [s] [inst_ty true ty_bs t s]]))] + else vals, Pretty.str s) + end + | Some ((_, s), _) => (vals, Pretty.str s)); + +fun get_cases sg = + Symtab.foldl (fn (tab, (k, {case_rewrites, ...})) => Symtab.update_new + ((fst (dest_Const (head_of (fst (HOLogic.dest_eq (HOLogic.dest_Trueprop + (prop_of (hd case_rewrites))))))), map my_mk_meta_eq case_rewrites), tab)) + (Symtab.empty, DatatypePackage.get_datatypes_sg sg); + +fun decomp_case th = + let + val (lhs, _) = Logic.dest_equals (prop_of th); + val (f, ts) = strip_comb lhs; + val (us, u) = split_last ts; + val (Const (s, _), vs) = strip_comb u + in (us, s, vs, u) end; + +fun rename vs t = + let + fun mk_subst ((vs, subs), Var ((s, i), T)) = + let val s' = variant vs s + in if s = s' then (vs, subs) + else (s' :: vs, ((s, i), Var ((s', i), T)) :: subs) + end; + val (vs', subs) = foldl mk_subst ((vs, []), term_vars t) + in (vs', subst_Vars subs t) end; + +fun is_instance sg t u = t = subst_TVars_Vartab + (Type.typ_match (Sign.tsig_of sg) (Vartab.empty, + (fastype_of u, fastype_of t))) u handle Type.TYPE_MATCH => false; + +(* +fun lookup sg fs t = apsome snd (Library.find_first + (is_instance sg t o fst) fs); +*) + +fun lookup sg fs t = (case Library.find_first (is_instance sg t o fst) fs of + None => (bla := (t ins !bla); None) + | Some (_, x) => Some x); + +fun unint sg fs t = forall (is_none o lookup sg fs) (term_consts' t); + +fun mk_let s i xs ys = + Pretty.blk (0, [Pretty.blk (i, separate Pretty.fbrk (Pretty.str s :: xs)), + Pretty.fbrk, + Pretty.blk (i, ([Pretty.str "in", Pretty.fbrk] @ ys)), + Pretty.fbrk, Pretty.str "end"]); + +(*****************************************************************************) +(* Generate bindings for simplifying term t *) +(* mkeq: whether to generate reflexivity theorem for uninterpreted terms *) +(* fs: interpreted functions *) +(* ts: atomic terms *) +(* vs: used identifiers *) +(* vals: list of bindings of the form ((eq, ct), ps) where *) +(* eq: name of equational theorem *) +(* ct: name of simplified cterm *) +(* ps: ML code for creating the above two items *) +(*****************************************************************************) + +fun mk_simpl_code sg case_tab mkeq fs ts ty_bs thm_bs ((vs, vals), t) = + (case assoc (vals, t) of + Some ((eq, ct), ps) => (* binding already generated *) + if mkeq andalso eq="" then + let val eq' = variant vs "eq" + in ((eq' :: vs, overwrite (vals, + (t, ((eq', ct), ps @ [mk_refleq eq' ct])))), (eq', ct)) + end + else ((vs, vals), (eq, ct)) + | None => (case assoc (ts, t) of + Some v => (* atomic term *) + let val xs = if not (null (term_tvars t)) andalso is_Const t then + [mk_val [v] [inst_ty false ty_bs t v]] else [] + in + if mkeq then + let val eq = variant vs "eq" + in ((eq :: vs, vals @ + [(t, ((eq, v), xs @ [mk_refleq eq v]))]), (eq, v)) + end + else ((vs, if null xs then vals else vals @ + [(t, (("", v), xs))]), ("", v)) + end + | None => (* complex term *) + let val (f as Const (cname, _), us) = strip_comb t + in case Symtab.lookup (case_tab, cname) of + Some cases => (* case expression *) + let + val (us', u) = split_last us; + val b = unint sg fs u; + val ((vs1, vals1), (eq, ct)) = + mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs ((vs, vals), u); + val xs = variantlist (replicate (length us') "f", vs1); + val (vals2, ps) = foldl_map + (mk_cterm_code false ty_bs ts []) (vals1, us'); + val fvals = map (fn (x, p) => mk_val [x] [p]) (xs ~~ ps); + val uT = fastype_of u; + val (us'', _, _, u') = decomp_case (hd cases); + val (vs2, ty_bs', ty_vals) = mk_type_bindings + (mk_type_bindings ((vs1 @ xs, [], []), + (hd xs, fastype_of (hd us''))), (ct, fastype_of u')); + val insts1 = map mk_tyinst ty_bs'; + val i = length vals2; + + fun mk_case_code ((vs, vals), (f, (name, eqn))) = + let + val (fvs, cname, cvs, _) = decomp_case eqn; + val Ts = binder_types (fastype_of f); + val ys = variantlist (map (fst o fst o dest_Var) cvs, vs); + val cvs' = map Var (map (rpair 0) ys ~~ Ts); + val rs = cvs' ~~ cvs; + val lhs = list_comb (Const (cname, Ts ---> uT), cvs'); + val rhs = foldl betapply (f, cvs'); + val (vs', tm_bs, tm_vals) = decomp_term_code false + ((vs @ ys, [], []), (ct, lhs)); + val ((vs'', all_vals), (eq', ct')) = mk_simpl_code sg case_tab + false fs (tm_bs @ ts) ty_bs thm_bs ((vs', vals), rhs); + val (old_vals, eq_vals) = splitAt (i, all_vals); + val vs''' = vs @ filter (fn v => exists + (fn (_, ((v', _), _)) => v = v') old_vals) (vs'' \\ vs'); + val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(", + inst_ty false ty_bs' t (the (assoc (thm_bs, t))), Pretty.str ",", + Pretty.brk 1, Pretty.str (s ^ ")")]) ((fvs ~~ xs) @ + (map (fn (v, s) => (the (assoc (rs, v)), s)) tm_bs)); + val eq'' = if null insts1 andalso null insts2 then Pretty.str name + else parens (eq' <> "") [Pretty.str + (if null cvs then "Thm.instantiate" else "Drule.instantiate"), + Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1, + Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2, + Pretty.str ")", Pretty.brk 1, Pretty.str name]; + val eq''' = if eq' = "" then eq'' else + Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1, + eq'', Pretty.brk 1, Pretty.str eq'] + in + ((vs''', old_vals), Pretty.block [pretty_pattern false lhs, + Pretty.str " =>", + Pretty.brk 1, mk_let "let" 2 (tm_vals @ flat (map (snd o snd) eq_vals)) + [Pretty.str ("(" ^ ct' ^ ","), Pretty.brk 1, eq''', Pretty.str ")"]]) + end; + + val case_names = map (fn i => Sign.base_name cname ^ "_" ^ + string_of_int i) (1 upto length cases); + val ((vs3, vals3), case_ps) = foldl_map mk_case_code + ((vs2, vals2), us' ~~ (case_names ~~ cases)); + val eq' = variant vs3 "eq"; + val ct' = variant (eq' :: vs3) "ct"; + val eq'' = variant (eq' :: ct' :: vs3) "eq"; + val case_vals = + fvals @ ty_vals @ + [mk_val [ct', eq'] ([Pretty.str "(case", Pretty.brk 1, + Pretty.str ("term_of " ^ ct ^ " of"), Pretty.brk 1] @ + flat (separate [Pretty.brk 1, Pretty.str "| "] + (map single case_ps)) @ [Pretty.str ")"])] + in + if b then + ((eq' :: ct' :: vs3, vals3 @ + [(t, ((eq', ct'), case_vals))]), (eq', ct')) + else + let val ((vs4, vals4), (_, ctcase)) = mk_simpl_code sg case_tab false + fs ts ty_bs thm_bs ((eq' :: eq'' :: ct' :: vs3, vals3), f) + in + ((vs4, vals4 @ [(t, ((eq'', ct'), case_vals @ + [mk_val [eq''] [Pretty.str "Thm.transitive", Pretty.brk 1, + Pretty.str "(Thm.combination", Pretty.brk 1, + Pretty.str "(Thm.reflexive", Pretty.brk 1, + mk_apps "Thm.capply" true (Pretty.str ctcase) + (map Pretty.str xs), + Pretty.str ")", Pretty.brk 1, Pretty.str (eq ^ ")"), + Pretty.brk 1, Pretty.str eq']]))]), (eq'', ct')) + end + end + + | None => + let + val b = forall (unint sg fs) us; + val (q, eqs) = foldl_map + (mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs) ((vs, vals), us); + val ((vs', vals'), (eqf, ctf)) = if is_some (lookup sg fs f) andalso b + then (q, ("", "")) + else mk_simpl_code sg case_tab (not b) fs ts ty_bs thm_bs (q, f); + val ct = variant vs' "ct"; + val eq = variant (ct :: vs') "eq"; + val ctv = mk_val [ct] [mk_apps "Thm.capply" false + (Pretty.str ctf) (map (Pretty.str o snd) eqs)]; + fun combp b = mk_apps "Thm.combination" b + (Pretty.str eqf) (map (Pretty.str o fst) eqs) + in + case (lookup sg fs f, b) of + (None, true) => (* completely uninterpreted *) + if mkeq then ((ct :: eq :: vs', vals' @ + [(t, ((eq, ct), [ctv, mk_refleq eq ct]))]), (eq, ct)) + else ((ct :: vs', vals' @ [(t, (("", ct), [ctv]))]), ("", ct)) + | (None, false) => (* function uninterpreted *) + ((eq :: ct :: vs', vals' @ + [(t, ((eq, ct), [ctv, mk_val [eq] [combp false]]))]), (eq, ct)) + | (Some (s, _, _), true) => (* arguments uninterpreted *) + ((eq :: ct :: vs', vals' @ + [(t, ((eq, ct), [mk_val [ct, eq] (separate (Pretty.brk 1) + (Pretty.str s :: map (Pretty.str o snd) eqs))]))]), (eq, ct)) + | (Some (s, _, _), false) => (* function and arguments interpreted *) + let val eq' = variant (eq :: ct :: vs') "eq" + in ((eq' :: eq :: ct :: vs', vals' @ [(t, ((eq', ct), + [mk_val [ct, eq] (separate (Pretty.brk 1) + (Pretty.str s :: map (Pretty.str o snd) eqs)), + mk_val [eq'] [Pretty.str "Thm.transitive", Pretty.brk 1, + combp true, Pretty.brk 1, Pretty.str eq]]))]), (eq', ct)) + end + end + end)); + +fun lhs_of thm = fst (Logic.dest_equals (prop_of thm)); +fun rhs_of thm = snd (Logic.dest_equals (prop_of thm)); + +fun mk_funs_code sg case_tab fs fs' = + let + val case_thms = mapfilter (fn s => (case Symtab.lookup (case_tab, s) of + None => None + | Some thms => Some (unsuffix "_case" (Sign.base_name s) ^ ".cases", + map (fn i => Sign.base_name s ^ "_" ^ string_of_int i) + (1 upto length thms) ~~ thms))) + (foldr add_term_consts (map (prop_of o snd) + (flat (map (#3 o snd) fs')), [])); + val case_vals = map (fn (s, cs) => mk_vall (map fst cs) + [Pretty.str "map my_mk_meta_eq", Pretty.brk 1, + Pretty.str ("(thms \"" ^ s ^ "\")")]) case_thms; + val (vs, thm_bs, thm_vals) = foldl mk_term_bindings (([], [], []), + flat (map (map (apsnd prop_of) o #3 o snd) fs') @ + map (apsnd prop_of) (flat (map snd case_thms))); + + fun mk_fun_code (prfx, (fname, d, eqns)) = + let + val (f, ts) = strip_comb (lhs_of (snd (hd eqns))); + val args = variantlist (replicate (length ts) "ct", vs); + val (vs', ty_bs, ty_vals) = foldl mk_type_bindings + ((vs @ args, [], []), args ~~ map fastype_of ts); + val insts1 = map mk_tyinst ty_bs; + + fun mk_eqn_code (name, eqn) = + let + val (_, argts) = strip_comb (lhs_of eqn); + val (vs'', tm_bs, tm_vals) = foldl (decomp_term_code false) + ((vs', [], []), args ~~ argts); + val ((vs''', eq_vals), (eq, ct)) = mk_simpl_code sg case_tab false fs + (tm_bs @ filter_out (is_Var o fst) thm_bs) ty_bs thm_bs + ((vs'', []), rhs_of eqn); + val insts2 = map (fn (t, s) => Pretty.block [Pretty.str "(", + inst_ty false ty_bs t (the (assoc (thm_bs, t))), Pretty.str ",", Pretty.brk 1, + Pretty.str (s ^ ")")]) tm_bs + val eq' = if null insts1 andalso null insts2 then Pretty.str name + else parens (eq <> "") [Pretty.str "Thm.instantiate", + Pretty.brk 1, Pretty.str "(", Pretty.list "[" "]" insts1, + Pretty.str ",", Pretty.brk 1, Pretty.list "[" "]" insts2, + Pretty.str ")", Pretty.brk 1, Pretty.str name]; + val eq'' = if eq = "" then eq' else + Pretty.block [Pretty.str "Thm.transitive", Pretty.brk 1, + eq', Pretty.brk 1, Pretty.str eq] + in + Pretty.block [parens (length argts > 1) + (Pretty.commas (map (pretty_pattern false) argts)), + Pretty.str " =>", + Pretty.brk 1, mk_let "let" 2 (ty_vals @ tm_vals @ flat (map (snd o snd) eq_vals)) + [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1, eq'', Pretty.str ")"]] + end; + + val default = if d then + let + val Some s = assoc (thm_bs, f); + val ct = variant vs' "ct" + in [Pretty.brk 1, Pretty.str "handle", Pretty.brk 1, + Pretty.str "Match =>", Pretty.brk 1, mk_let "let" 2 + (ty_vals @ (if null (term_tvars f) then [] else + [mk_val [s] [inst_ty false ty_bs f s]]) @ + [mk_val [ct] [mk_apps "Thm.capply" false (Pretty.str s) + (map Pretty.str args)]]) + [Pretty.str ("(" ^ ct ^ ","), Pretty.brk 1, + Pretty.str "Thm.reflexive", Pretty.brk 1, Pretty.str (ct ^ ")")]] + end + else [] + in + ("and ", Pretty.block (separate (Pretty.brk 1) + (Pretty.str (prfx ^ fname) :: map Pretty.str args) @ + [Pretty.str " =", Pretty.brk 1, Pretty.str "(case", Pretty.brk 1, + Pretty.list "(" ")" (map (fn s => Pretty.str ("term_of " ^ s)) args), + Pretty.str " of", Pretty.brk 1] @ + flat (separate [Pretty.brk 1, Pretty.str "| "] + (map (single o mk_eqn_code) eqns)) @ [Pretty.str ")"] @ default)) + end; + + val (_, decls) = foldl_map mk_fun_code ("fun ", map snd fs') + in + mk_let "local" 2 (case_vals @ thm_vals) (separate Pretty.fbrk decls) + end; + +fun mk_simprocs_code sg eqns = + let + val case_tab = get_cases sg; + fun get_head th = head_of (fst (Logic.dest_equals (prop_of th))); + fun attach_term (x as (_, _, (_, th) :: _)) = (get_head th, x); + val eqns' = map attach_term eqns; + fun mk_node (s, _, (_, th) :: _) = (s, get_head th); + fun mk_edges (s, _, ths) = map (pair s) (distinct + (mapfilter (fn t => apsome #1 (lookup sg eqns' t)) + (flat (map (term_consts' o prop_of o snd) ths)))); + val gr = foldr (uncurry Graph.add_edge) + (map (pair "" o #1) eqns @ flat (map mk_edges eqns), + foldr (uncurry Graph.new_node) + (("", Bound 0) :: map mk_node eqns, Graph.empty)); + val keys = rev (Graph.all_succs gr [""] \ ""); + fun gr_ord (x :: _, y :: _) = + int_ord (find_index (equal x) keys, find_index (equal y) keys); + val scc = map (fn xs => filter (fn (_, (s, _, _)) => s mem xs) eqns') + (sort gr_ord (Graph.strong_conn gr \ [""])); + in + flat (separate [Pretty.str ";", Pretty.fbrk, Pretty.str " ", Pretty.fbrk] + (map (fn eqns'' => [mk_funs_code sg case_tab eqns' eqns'']) scc)) @ + [Pretty.str ";", Pretty.fbrk] + end; + +fun use_simprocs_code sg eqns = + let + fun attach_name (i, x) = (i+1, ("simp_thm_" ^ string_of_int i, x)); + fun attach_names (i, (s, b, eqs)) = + let val (i', eqs') = foldl_map attach_name (i, eqs) + in (i', (s, b, eqs')) end; + val (_, eqns') = foldl_map attach_names (1, eqns); + val (names, thms) = split_list (flat (map #3 eqns')); + val s = setmp print_mode [] Pretty.string_of + (mk_let "local" 2 [mk_vall names [Pretty.str "!SimprocsCodegen.simp_thms"]] + (mk_simprocs_code sg eqns')) + in + (simp_thms := thms; use_text Context.ml_output false s) + end; + +end; diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Matrix/fspmlp.ML --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/HOL/Matrix/fspmlp.ML Fri Sep 03 17:10:36 2004 +0200 @@ -0,0 +1,303 @@ +signature FSPMLP = +sig + type linprog + + val y : linprog -> cterm + val A : linprog -> cterm * cterm + val b : linprog -> cterm + val c : linprog -> cterm * cterm + val r : linprog -> cterm + + exception Load of string + + val load : string -> int -> bool -> linprog +end + +structure fspmlp : FSPMLP = +struct + +type linprog = cterm * (cterm * cterm) * cterm * (cterm * cterm) * cterm + +fun y (c1, c2, c3, c4, c5) = c1 +fun A (c1, c2, c3, c4, c5) = c2 +fun b (c1, c2, c3, c4, c5) = c3 +fun c (c1, c2, c3, c4, c5) = c4 +fun r (c1, c2, c3, c4, c5) = c5 + +structure CplexFloatSparseMatrixConverter = +MAKE_CPLEX_MATRIX_CONVERTER(structure cplex = Cplex and matrix_builder = FloatSparseMatrixBuilder); + +datatype bound_type = LOWER | UPPER + +fun intbound_ord ((i1, b1),(i2,b2)) = + if i1 < i2 then LESS + else if i1 = i2 then + (if b1 = b2 then EQUAL else if b1=LOWER then LESS else GREATER) + else GREATER + +structure Inttab = TableFun(type key = int val ord = (rev_order o int_ord)); + +structure VarGraph = TableFun(type key = int*bound_type val ord = intbound_ord); +(* key -> (float option) * (int -> (float * (((float * float) * key) list)))) *) +(* dest_key -> (sure_bound * (row_index -> (row_bound * (((coeff_lower * coeff_upper) * src_key) list)))) *) + +exception Internal of string; + +fun add_row_bound g dest_key row_index row_bound = + let + val x = + case VarGraph.lookup (g, dest_key) of + None => (None, Inttab.update ((row_index, (row_bound, [])), Inttab.empty)) + | Some (sure_bound, f) => + (sure_bound, + case Inttab.lookup (f, row_index) of + None => Inttab.update ((row_index, (row_bound, [])), f) + | Some _ => raise (Internal "add_row_bound")) + in + VarGraph.update ((dest_key, x), g) + end + +fun update_sure_bound g (key as (_, btype)) bound = + let + val x = + case VarGraph.lookup (g, key) of + None => (Some bound, Inttab.empty) + | Some (None, f) => (Some bound, f) + | Some (Some old_bound, f) => + (Some ((case btype of + UPPER => FloatArith.min + | LOWER => FloatArith.max) + old_bound bound), f) + in + VarGraph.update ((key, x), g) + end + +fun get_sure_bound g key = + case VarGraph.lookup (g, key) of + None => None + | Some (sure_bound, _) => sure_bound + +(*fun get_row_bound g key row_index = + case VarGraph.lookup (g, key) of + None => None + | Some (sure_bound, f) => + (case Inttab.lookup (f, row_index) of + None => None + | Some (row_bound, _) => (sure_bound, row_bound))*) + +fun add_edge g src_key dest_key row_index coeff = + case VarGraph.lookup (g, dest_key) of + None => raise (Internal "add_edge: dest_key not found") + | Some (sure_bound, f) => + (case Inttab.lookup (f, row_index) of + None => raise (Internal "add_edge: row_index not found") + | Some (row_bound, sources) => + VarGraph.update ((dest_key, (sure_bound, Inttab.update ((row_index, (row_bound, (coeff, src_key) :: sources)), f))), g)) + +fun split_graph g = + let + fun split ((r1, r2), (key, (sure_bound, _))) = + case sure_bound of + None => (r1, r2) + | Some bound => + (case key of + (u, UPPER) => (r1, Inttab.update ((u, bound), r2)) + | (u, LOWER) => (Inttab.update ((u, bound), r1), r2)) + in + VarGraph.foldl split ((Inttab.empty, Inttab.empty), g) + end + +fun it2list t = + let + fun tolist (l, a) = a::l + in + Inttab.foldl tolist ([], t) + end + +(* If safe is true, termination is guaranteed, but the sure bounds may be not optimal (relative to the algorithm). + If safe is false, termination is not guaranteed, but on termination the sure bounds are optimal (relative to the algorithm) *) +fun propagate_sure_bounds safe names g = + let + (* returns None if no new sure bound could be calculated, otherwise the new sure bound is returned *) + fun calc_sure_bound_from_sources g (key as (_, btype)) = + let + fun mult_upper x (lower, upper) = + if FloatArith.is_negative x then + FloatArith.mul x lower + else + FloatArith.mul x upper + + fun mult_lower x (lower, upper) = + if FloatArith.is_negative x then + FloatArith.mul x upper + else + FloatArith.mul x lower + + val mult_btype = case btype of UPPER => mult_upper | LOWER => mult_lower + + fun calc_sure_bound (sure_bound, (row_index, (row_bound, sources))) = + let + fun add_src_bound (sum, (coeff, src_key)) = + case sum of + None => None + | Some x => + (case get_sure_bound g src_key of + None => None + | Some src_sure_bound => Some (FloatArith.add x (mult_btype src_sure_bound coeff))) + in + case foldl add_src_bound (Some row_bound, sources) of + None => sure_bound + | new_sure_bound as (Some new_bound) => + (case sure_bound of + None => new_sure_bound + | Some old_bound => + Some (case btype of + UPPER => FloatArith.min old_bound new_bound + | LOWER => FloatArith.max old_bound new_bound)) + end + in + case VarGraph.lookup (g, key) of + None => None + | Some (sure_bound, f) => + let + val x = Inttab.foldl calc_sure_bound (sure_bound, f) + in + if x = sure_bound then None else x + end + end + + fun propagate ((g, b), (key, _)) = + case calc_sure_bound_from_sources g key of + None => (g,b) + | Some bound => (update_sure_bound g key bound, + if safe then + case get_sure_bound g key of + None => true + | _ => b + else + true) + + val (g, b) = VarGraph.foldl propagate ((g, false), g) + in + if b then propagate_sure_bounds safe names g else g + end + +exception Load of string; + +fun calcr safe_propagation xlen names prec A b = + let + val empty = Inttab.empty + + fun instab t i x = Inttab.update ((i,x), t) + + fun isnegstr x = String.isPrefix "-" x + fun negstr x = if isnegstr x then String.extract (x, 1, NONE) else "-"^x + + fun test_1 (lower, upper) = + if lower = upper then + (if FloatArith.is_equal lower (IntInf.fromInt ~1, FloatArith.izero) then ~1 + else if FloatArith.is_equal lower (IntInf.fromInt 1, FloatArith.izero) then 1 + else 0) + else 0 + + fun calcr (g, (row_index, a)) = + let + val b = FloatSparseMatrixBuilder.v_elem_at b row_index + val (_, b2) = ExactFloatingPoint.approx_decstr_by_bin prec (case b of None => "0" | Some b => b) + val approx_a = FloatSparseMatrixBuilder.v_fold (fn (l, (i,s)) => + (i, ExactFloatingPoint.approx_decstr_by_bin prec s)::l) [] a + + fun fold_dest_nodes (g, (dest_index, dest_value)) = + let + val dest_test = test_1 dest_value + in + if dest_test = 0 then + g + else let + val (dest_key as (_, dest_btype), row_bound) = + if dest_test = ~1 then + ((dest_index, LOWER), FloatArith.neg b2) + else + ((dest_index, UPPER), b2) + + fun fold_src_nodes (g, (src_index, src_value as (src_lower, src_upper))) = + if src_index = dest_index then g + else + let + val coeff = case dest_btype of + UPPER => (FloatArith.neg src_upper, FloatArith.neg src_lower) + | LOWER => src_value + in + if FloatArith.is_negative src_lower then + add_edge g (src_index, UPPER) dest_key row_index coeff + else + add_edge g (src_index, LOWER) dest_key row_index coeff + end + in + foldl fold_src_nodes ((add_row_bound g dest_key row_index row_bound), approx_a) + end + end + in + case approx_a of + [] => g + | [(u, a)] => + let + val atest = test_1 a + in + if atest = ~1 then + update_sure_bound g (u, LOWER) (FloatArith.neg b2) + else if atest = 1 then + update_sure_bound g (u, UPPER) b2 + else + g + end + | _ => foldl fold_dest_nodes (g, approx_a) + end + + val g = FloatSparseMatrixBuilder.m_fold calcr VarGraph.empty A + val g = propagate_sure_bounds safe_propagation names g + + val (r1, r2) = split_graph g + + fun abs_estimate i r1 r2 = + if i = 0 then FloatSparseMatrixBuilder.empty_spmat + else + let + val index = xlen-i + val r = abs_estimate (i-1) r1 r2 + val b1 = case Inttab.lookup (r1, index) of None => raise (Load ("x-value not bounded from below: "^(names index))) | Some x => x + val b2 = case Inttab.lookup (r2, index) of None => raise (Load ("x-value not bounded from above: "^(names index))) | Some x => x + val abs_max = FloatArith.max (FloatArith.neg (FloatArith.negative_part b1)) (FloatArith.positive_part b2) + val vec = FloatSparseMatrixBuilder.cons_spvec (FloatSparseMatrixBuilder.mk_spvec_entry 0 abs_max) FloatSparseMatrixBuilder.empty_spvec + in + FloatSparseMatrixBuilder.cons_spmat (FloatSparseMatrixBuilder.mk_spmat_entry index vec) r + end + in + FloatSparseMatrixBuilder.sign_term (abs_estimate xlen r1 r2) + end + +fun load filename prec safe_propagation = + let + val prog = Cplex.load_cplexFile filename + val prog = Cplex.elim_nonfree_bounds prog + val prog = Cplex.relax_strict_ineqs prog + val (maximize, c, A, b, (xlen, names, _)) = CplexFloatSparseMatrixConverter.convert_prog prog + val r = calcr safe_propagation xlen names prec A b + val _ = if maximize then () else raise Load "sorry, cannot handle minimization problems" + val (dualprog, indexof) = FloatSparseMatrixBuilder.dual_cplexProg c A b + val results = Cplex.solve dualprog + val (optimal,v) = CplexFloatSparseMatrixConverter.convert_results results indexof + val A = FloatSparseMatrixBuilder.cut_matrix v None A + fun id x = x + val v = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 v + val b = FloatSparseMatrixBuilder.transpose_matrix (FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 b) + val c = FloatSparseMatrixBuilder.set_vector FloatSparseMatrixBuilder.empty_matrix 0 c + val (y1, _) = FloatSparseMatrixBuilder.approx_matrix prec FloatArith.positive_part v + val A = FloatSparseMatrixBuilder.approx_matrix prec id A + val (_,b2) = FloatSparseMatrixBuilder.approx_matrix prec id b + val c = FloatSparseMatrixBuilder.approx_matrix prec id c + in + (y1, A, b2, c, r) + end handle CplexFloatSparseMatrixConverter.Converter s => (raise (Load ("Converter: "^s))) + +end \ No newline at end of file diff -r e7616269fdca -r 5f621aa35c25 src/HOL/OrderedGroup.thy --- a/src/HOL/OrderedGroup.thy Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/OrderedGroup.thy Fri Sep 03 17:10:36 2004 +0200 @@ -842,6 +842,42 @@ lemma minus_add_cancel: "-(a::'a::ab_group_add) + (a + b) = b" by (simp add: add_assoc[symmetric]) +lemma le_add_right_mono: + assumes + "a <= b + (c::'a::pordered_ab_group_add)" + "c <= d" + shows "a <= b + d" + apply (rule_tac order_trans[where y = "b+c"]) + apply (simp_all add: prems) + done + +lemmas group_eq_simps = + mult_ac + add_ac + add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2 + diff_eq_eq eq_diff_eq + +lemma estimate_by_abs: +"a + b <= (c::'a::lordered_ab_group_abs) \ a <= c + abs b" +proof - + assume 1: "a+b <= c" + have 2: "a <= c+(-b)" + apply (insert 1) + apply (drule_tac add_right_mono[where c="-b"]) + apply (simp add: group_eq_simps) + done + have 3: "(-b) <= abs b" by (rule abs_ge_minus_self) + show ?thesis by (rule le_add_right_mono[OF 2 3]) +qed + +lemma abs_of_ge_0: "0 <= (y::'a::lordered_ab_group_abs) \ abs y = y" +proof - + assume 1:"0 <= y" + have 2:"-y <= 0" by (simp add: 1) + from 1 2 have 3:"-y <= y" by (simp only:) + show ?thesis by (simp add: abs_lattice ge_imp_join_eq[OF 3]) +qed + ML {* val add_zero_left = thm"add_0"; val add_zero_right = thm"add_0_right"; diff -r e7616269fdca -r 5f621aa35c25 src/HOL/Ring_and_Field.thy --- a/src/HOL/Ring_and_Field.thy Fri Sep 03 10:27:05 2004 +0200 +++ b/src/HOL/Ring_and_Field.thy Fri Sep 03 17:10:36 2004 +0200 @@ -546,14 +546,14 @@ done text{*This list of rewrites decides ring equalities by ordered rewriting.*} -lemmas ring_eq_simps = - mult_ac +lemmas ring_eq_simps = +(* mult_ac*) left_distrib right_distrib left_diff_distrib right_diff_distrib - add_ac + group_eq_simps +(* add_ac add_diff_eq diff_add_eq diff_diff_eq diff_diff_eq2 - diff_eq_eq eq_diff_eq + diff_eq_eq eq_diff_eq *) - subsection {* Fields *} lemma right_inverse [simp]: @@ -1503,6 +1503,90 @@ text{*Moving this up spoils many proofs using @{text mult_le_cancel_right}*} declare times_divide_eq_left [simp] +lemma linprog_dual_estimate: + assumes + "A * x \ (b::'a::lordered_ring)" + "0 \ y" + "abs (A - A') \ \A" + "b \ b'" + "abs (c - c') \ \c" + "abs x \ r" + shows + "c * x \ y * b' + (y * \A + abs (y * A' - c') + \c) * r" +proof - + from prems have 1: "y * b <= y * b'" by (simp add: mult_left_mono) + from prems have 2: "y * (A * x) <= y * b" by (simp add: mult_left_mono) + have 3: "y * (A * x) = c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x" by (simp add: ring_eq_simps) + from 1 2 3 have 4: "c * x + (y * (A - A') + (y * A' - c') + (c'-c)) * x <= y * b'" by simp + have 5: "c * x <= y * b' + abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)" + by (simp only: 4 estimate_by_abs) + have 6: "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <= abs (y * (A - A') + (y * A' - c') + (c'-c)) * abs x" + by (simp add: abs_le_mult) + 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" + by (simp add: abs_triangle_ineq mult_right_mono) + 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" + by (simp add: abs_triangle_ineq mult_right_mono) + 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" + by (simp add: abs_le_mult mult_right_mono) + have 10: "c'-c = -(c-c')" by (simp add: ring_eq_simps) + have 11: "abs (c'-c) = abs (c-c')" + by (subst 10, subst abs_minus_cancel, simp) + 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') + \c) * abs x" + by (simp add: 11 prems mult_right_mono) + have 13: "(abs y * abs (A-A') + abs (y*A'-c') + \c) * abs x <= (abs y * \A + abs (y*A'-c') + \c) * abs x" + by (simp add: prems mult_right_mono mult_left_mono) + have r: "(abs y * \A + abs (y*A'-c') + \c) * abs x <= (abs y * \A + abs (y*A'-c') + \c) * r" + apply (rule mult_left_mono) + apply (simp add: prems) + apply (rule_tac add_mono[of "0::'a" _ "0", simplified])+ + apply (rule mult_left_mono[of "0" "\A", simplified]) + apply (simp_all) + apply (rule order_trans[where y="abs (A-A')"], simp_all add: prems) + apply (rule order_trans[where y="abs (c-c')"], simp_all add: prems) + done + from 6 7 8 9 12 13 r have 14:" abs((y * (A - A') + (y * A' - c') + (c'-c)) * x) <=(abs y * \A + abs (y*A'-c') + \c) * r" + by (simp) + show ?thesis + apply (rule_tac le_add_right_mono[of _ _ "abs((y * (A - A') + (y * A' - c') + (c'-c)) * x)"]) + apply (simp_all add: 5 14[simplified abs_of_ge_0[of y, simplified prems]]) + done +qed + +lemma le_ge_imp_abs_diff_1: + assumes + "A1 <= (A::'a::lordered_ring)" + "A <= A2" + shows "abs (A-A1) <= A2-A1" +proof - + have "0 <= A - A1" + proof - + have 1: "A - A1 = A + (- A1)" by simp + show ?thesis by (simp only: 1 add_right_mono[of A1 A "-A1", simplified, simplified prems]) + qed + then have "abs (A-A1) = A-A1" by (rule abs_of_ge_0) + with prems show "abs (A-A1) <= (A2-A1)" by simp +qed + +lemma linprog_dual_estimate_1: + assumes + "A * x \ (b::'a::lordered_ring)" + "0 \ y" + "A1 <= A" + "A <= A2" + "c1 <= c" + "c <= c2" + "abs x \ r" + shows + "c * x \ y * b + (y * (A2 - A1) + abs (y * A1 - c1) + (c2 - c1)) * r" +proof - + from prems have delta_A: "abs (A-A1) <= (A2-A1)" by (simp add: le_ge_imp_abs_diff_1) + from prems have delta_c: "abs (c-c1) <= (c2-c1)" by (simp add: le_ge_imp_abs_diff_1) + show ?thesis + apply (rule_tac linprog_dual_estimate) + apply (auto intro: delta_A delta_c simp add: prems) + done +qed + ML {* val left_distrib = thm "left_distrib"; val right_distrib = thm "right_distrib";