src/HOL/Library/normarith.ML
author boehmes
Wed Aug 26 11:40:28 2009 +0200 (2009-08-26)
changeset 32402 5731300da417
parent 31446 2d91b2416de8
child 32645 1cc5b24f5a01
permissions -rw-r--r--
added further conversions and conversionals
     1 (* Title:      Library/normarith.ML
     2    Author:     Amine Chaieb, University of Cambridge
     3    Description: A simple decision procedure for linear problems in euclidean space
     4 *)
     5 
     6   (* Now the norm procedure for euclidean spaces *)
     7 
     8 
     9 signature NORM_ARITH = 
    10 sig
    11  val norm_arith : Proof.context -> conv
    12  val norm_arith_tac : Proof.context -> int -> tactic
    13 end
    14 
    15 structure NormArith : NORM_ARITH = 
    16 struct
    17 
    18  open Conv Thm;
    19  val bool_eq = op = : bool *bool -> bool
    20   fun dest_ratconst t = case term_of t of
    21    Const(@{const_name divide}, _)$a$b => Rat.rat_of_quotient(HOLogic.dest_number a |> snd, HOLogic.dest_number b |> snd)
    22  | Const(@{const_name inverse}, _)$a => Rat.rat_of_quotient(1, HOLogic.dest_number a |> snd)
    23  | _ => Rat.rat_of_int (HOLogic.dest_number (term_of t) |> snd)
    24  fun is_ratconst t = can dest_ratconst t
    25  fun augment_norm b t acc = case term_of t of 
    26      Const(@{const_name norm}, _) $ _ => insert (eq_pair bool_eq (op aconvc)) (b,dest_arg t) acc
    27    | _ => acc
    28  fun find_normedterms t acc = case term_of t of
    29     @{term "op + :: real => _"}$_$_ =>
    30             find_normedterms (dest_arg1 t) (find_normedterms (dest_arg t) acc)
    31       | @{term "op * :: real => _"}$_$n =>
    32             if not (is_ratconst (dest_arg1 t)) then acc else
    33             augment_norm (dest_ratconst (dest_arg1 t) >=/ Rat.zero) 
    34                       (dest_arg t) acc
    35       | _ => augment_norm true t acc 
    36 
    37  val cterm_lincomb_neg = Ctermfunc.mapf Rat.neg
    38  fun cterm_lincomb_cmul c t = 
    39     if c =/ Rat.zero then Ctermfunc.undefined else Ctermfunc.mapf (fn x => x */ c) t
    40  fun cterm_lincomb_add l r = Ctermfunc.combine (curry op +/) (fn x => x =/ Rat.zero) l r
    41  fun cterm_lincomb_sub l r = cterm_lincomb_add l (cterm_lincomb_neg r)
    42  fun cterm_lincomb_eq l r = Ctermfunc.is_undefined (cterm_lincomb_sub l r)
    43 
    44  val int_lincomb_neg = Intfunc.mapf Rat.neg
    45  fun int_lincomb_cmul c t = 
    46     if c =/ Rat.zero then Intfunc.undefined else Intfunc.mapf (fn x => x */ c) t
    47  fun int_lincomb_add l r = Intfunc.combine (curry op +/) (fn x => x =/ Rat.zero) l r
    48  fun int_lincomb_sub l r = int_lincomb_add l (int_lincomb_neg r)
    49  fun int_lincomb_eq l r = Intfunc.is_undefined (int_lincomb_sub l r)
    50 
    51 fun vector_lincomb t = case term_of t of 
    52    Const(@{const_name plus}, _) $ _ $ _ =>
    53     cterm_lincomb_add (vector_lincomb (dest_arg1 t)) (vector_lincomb (dest_arg t))
    54  | Const(@{const_name minus}, _) $ _ $ _ =>
    55     cterm_lincomb_sub (vector_lincomb (dest_arg1 t)) (vector_lincomb (dest_arg t))
    56  | Const(@{const_name scaleR}, _)$_$_ =>
    57     cterm_lincomb_cmul (dest_ratconst (dest_arg1 t)) (vector_lincomb (dest_arg t))
    58  | Const(@{const_name uminus}, _)$_ =>
    59      cterm_lincomb_neg (vector_lincomb (dest_arg t))
    60 (* FIXME: how should we handle numerals?
    61  | Const(@ {const_name vec},_)$_ => 
    62    let 
    63      val b = ((snd o HOLogic.dest_number o term_of o dest_arg) t = 0 
    64                handle TERM _=> false)
    65    in if b then Ctermfunc.onefunc (t,Rat.one)
    66       else Ctermfunc.undefined
    67    end
    68 *)
    69  | _ => Ctermfunc.onefunc (t,Rat.one)
    70 
    71  fun vector_lincombs ts =
    72   fold_rev 
    73    (fn t => fn fns => case AList.lookup (op aconvc) fns t of
    74      NONE => 
    75        let val f = vector_lincomb t 
    76        in case find_first (fn (_,f') => cterm_lincomb_eq f f') fns of
    77            SOME (_,f') => (t,f') :: fns
    78          | NONE => (t,f) :: fns 
    79        end
    80    | SOME _ => fns) ts []
    81 
    82 fun replacenegnorms cv t = case term_of t of 
    83   @{term "op + :: real => _"}$_$_ => binop_conv (replacenegnorms cv) t
    84 | @{term "op * :: real => _"}$_$_ => 
    85     if dest_ratconst (dest_arg1 t) </ Rat.zero then arg_conv cv t else reflexive t
    86 | _ => reflexive t
    87 fun flip v eq = 
    88   if Ctermfunc.defined eq v 
    89   then Ctermfunc.update (v, Rat.neg (Ctermfunc.apply eq v)) eq else eq
    90 fun allsubsets s = case s of 
    91   [] => [[]]
    92 |(a::t) => let val res = allsubsets t in
    93                map (cons a) res @ res end
    94 fun evaluate env lin =
    95  Intfunc.fold (fn (x,c) => fn s => s +/ c */ (Intfunc.apply env x)) 
    96    lin Rat.zero
    97 
    98 fun solve (vs,eqs) = case (vs,eqs) of
    99   ([],[]) => SOME (Intfunc.onefunc (0,Rat.one))
   100  |(_,eq::oeqs) => 
   101    (case filter (member (op =) vs) (Intfunc.dom eq) of (*FIXME use find_first here*)
   102      [] => NONE
   103     | v::_ => 
   104        if Intfunc.defined eq v 
   105        then 
   106         let 
   107          val c = Intfunc.apply eq v
   108          val vdef = int_lincomb_cmul (Rat.neg (Rat.inv c)) eq
   109          fun eliminate eqn = if not (Intfunc.defined eqn v) then eqn 
   110                              else int_lincomb_add (int_lincomb_cmul (Intfunc.apply eqn v) vdef) eqn
   111         in (case solve (vs \ v,map eliminate oeqs) of
   112             NONE => NONE
   113           | SOME soln => SOME (Intfunc.update (v, evaluate soln (Intfunc.undefine v vdef)) soln))
   114         end
   115        else NONE)
   116 
   117 fun combinations k l = if k = 0 then [[]] else
   118  case l of 
   119   [] => []
   120 | h::t => map (cons h) (combinations (k - 1) t) @ combinations k t
   121 
   122 
   123 fun forall2 p l1 l2 = case (l1,l2) of 
   124    ([],[]) => true
   125  | (h1::t1,h2::t2) => p h1 h2 andalso forall2 p t1 t2
   126  | _ => false;
   127 
   128 
   129 fun vertices vs eqs =
   130  let 
   131   fun vertex cmb = case solve(vs,cmb) of
   132     NONE => NONE
   133    | SOME soln => SOME (map (fn v => Intfunc.tryapplyd soln v Rat.zero) vs)
   134   val rawvs = map_filter vertex (combinations (length vs) eqs)
   135   val unset = filter (forall (fn c => c >=/ Rat.zero)) rawvs 
   136  in fold_rev (insert (uncurry (forall2 (curry op =/)))) unset [] 
   137  end 
   138 
   139 fun subsumes l m = forall2 (fn x => fn y => Rat.abs x <=/ Rat.abs y) l m 
   140 
   141 fun subsume todo dun = case todo of
   142  [] => dun
   143 |v::ovs => 
   144    let val dun' = if exists (fn w => subsumes w v) dun then dun
   145                   else v::(filter (fn w => not(subsumes v w)) dun) 
   146    in subsume ovs dun' 
   147    end;
   148 
   149 fun match_mp PQ P = P RS PQ;
   150 
   151 fun cterm_of_rat x = 
   152 let val (a, b) = Rat.quotient_of_rat x
   153 in 
   154  if b = 1 then Numeral.mk_cnumber @{ctyp "real"} a
   155   else Thm.capply (Thm.capply @{cterm "op / :: real => _"} 
   156                    (Numeral.mk_cnumber @{ctyp "real"} a))
   157         (Numeral.mk_cnumber @{ctyp "real"} b)
   158 end;
   159 
   160 fun norm_cmul_rule c th = instantiate' [] [SOME (cterm_of_rat c)] (th RS @{thm norm_cmul_rule_thm});
   161 
   162 fun norm_add_rule th1 th2 = [th1, th2] MRS @{thm norm_add_rule_thm};
   163 
   164   (* I think here the static context should be sufficient!! *)
   165 fun inequality_canon_rule ctxt = 
   166  let 
   167   (* FIXME : Should be computed statically!! *)
   168   val real_poly_conv = 
   169     Normalizer.semiring_normalize_wrapper ctxt
   170      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
   171  in fconv_rule (arg_conv ((rewr_conv @{thm ge_iff_diff_ge_0}) then_conv arg_conv (field_comp_conv then_conv real_poly_conv)))
   172 end;
   173 
   174  fun absc cv ct = case term_of ct of 
   175  Abs (v,_, _) => 
   176   let val (x,t) = Thm.dest_abs (SOME v) ct
   177   in Thm.abstract_rule ((fst o dest_Free o term_of) x) x (cv t)
   178   end
   179  | _ => all_conv ct;
   180 
   181 fun sub_conv cv ct = (comb_conv cv else_conv absc cv) ct;
   182 fun botc1 conv ct = 
   183   ((sub_conv (botc1 conv)) then_conv (conv else_conv all_conv)) ct;
   184 
   185  fun rewrs_conv eqs ct = first_conv (map rewr_conv eqs) ct;
   186  val apply_pth1 = rewr_conv @{thm pth_1};
   187  val apply_pth2 = rewr_conv @{thm pth_2};
   188  val apply_pth3 = rewr_conv @{thm pth_3};
   189  val apply_pth4 = rewrs_conv @{thms pth_4};
   190  val apply_pth5 = rewr_conv @{thm pth_5};
   191  val apply_pth6 = rewr_conv @{thm pth_6};
   192  val apply_pth7 = rewrs_conv @{thms pth_7};
   193  val apply_pth8 = rewr_conv @{thm pth_8} then_conv arg1_conv field_comp_conv then_conv (try_conv (rewr_conv (mk_meta_eq @{thm scaleR_zero_left})));
   194  val apply_pth9 = rewrs_conv @{thms pth_9} then_conv arg1_conv (arg1_conv field_comp_conv);
   195  val apply_ptha = rewr_conv @{thm pth_a};
   196  val apply_pthb = rewrs_conv @{thms pth_b};
   197  val apply_pthc = rewrs_conv @{thms pth_c};
   198  val apply_pthd = try_conv (rewr_conv @{thm pth_d});
   199 
   200 fun headvector t = case t of 
   201   Const(@{const_name plus}, _)$
   202    (Const(@{const_name scaleR}, _)$l$v)$r => v
   203  | Const(@{const_name scaleR}, _)$l$v => v
   204  | _ => error "headvector: non-canonical term"
   205 
   206 fun vector_cmul_conv ct =
   207    ((apply_pth5 then_conv arg1_conv field_comp_conv) else_conv
   208     (apply_pth6 then_conv binop_conv vector_cmul_conv)) ct
   209 
   210 fun vector_add_conv ct = apply_pth7 ct 
   211  handle CTERM _ => 
   212   (apply_pth8 ct 
   213    handle CTERM _ => 
   214     (case term_of ct of 
   215      Const(@{const_name plus},_)$lt$rt =>
   216       let 
   217        val l = headvector lt 
   218        val r = headvector rt
   219       in (case TermOrd.fast_term_ord (l,r) of
   220          LESS => (apply_pthb then_conv arg_conv vector_add_conv 
   221                   then_conv apply_pthd) ct
   222         | GREATER => (apply_pthc then_conv arg_conv vector_add_conv 
   223                      then_conv apply_pthd) ct 
   224         | EQUAL => (apply_pth9 then_conv 
   225                 ((apply_ptha then_conv vector_add_conv) else_conv 
   226               arg_conv vector_add_conv then_conv apply_pthd)) ct)
   227       end
   228      | _ => reflexive ct))
   229 
   230 fun vector_canon_conv ct = case term_of ct of
   231  Const(@{const_name plus},_)$_$_ =>
   232   let 
   233    val ((p,l),r) = Thm.dest_comb ct |>> Thm.dest_comb
   234    val lth = vector_canon_conv l 
   235    val rth = vector_canon_conv r
   236    val th = Drule.binop_cong_rule p lth rth 
   237   in fconv_rule (arg_conv vector_add_conv) th end
   238 
   239 | Const(@{const_name scaleR}, _)$_$_ =>
   240   let 
   241    val (p,r) = Thm.dest_comb ct
   242    val rth = Drule.arg_cong_rule p (vector_canon_conv r) 
   243   in fconv_rule (arg_conv (apply_pth4 else_conv vector_cmul_conv)) rth
   244   end
   245 
   246 | Const(@{const_name minus},_)$_$_ => (apply_pth2 then_conv vector_canon_conv) ct
   247 
   248 | Const(@{const_name uminus},_)$_ => (apply_pth3 then_conv vector_canon_conv) ct
   249 
   250 (* FIXME
   251 | Const(@{const_name vec},_)$n => 
   252   let val n = Thm.dest_arg ct
   253   in if is_ratconst n andalso not (dest_ratconst n =/ Rat.zero) 
   254      then reflexive ct else apply_pth1 ct
   255   end
   256 *)
   257 | _ => apply_pth1 ct
   258 
   259 fun norm_canon_conv ct = case term_of ct of
   260   Const(@{const_name norm},_)$_ => arg_conv vector_canon_conv ct
   261  | _ => raise CTERM ("norm_canon_conv", [ct])
   262 
   263 fun fold_rev2 f [] [] z = z
   264  | fold_rev2 f (x::xs) (y::ys) z = f x y (fold_rev2 f xs ys z)
   265  | fold_rev2 f _ _ _ = raise UnequalLengths;
   266 
   267 fun int_flip v eq = 
   268   if Intfunc.defined eq v 
   269   then Intfunc.update (v, Rat.neg (Intfunc.apply eq v)) eq else eq;
   270 
   271 local
   272  val pth_zero = @{thm norm_zero}
   273  val tv_n = (ctyp_of_term o dest_arg o dest_arg1 o dest_arg o cprop_of)
   274              pth_zero
   275  val concl = dest_arg o cprop_of 
   276  fun real_vector_combo_prover ctxt translator (nubs,ges,gts) = 
   277   let 
   278    (* FIXME: Should be computed statically!!*)
   279    val real_poly_conv = 
   280       Normalizer.semiring_normalize_wrapper ctxt
   281        (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
   282    val sources = map (dest_arg o dest_arg1 o concl) nubs
   283    val rawdests = fold_rev (find_normedterms o dest_arg o concl) (ges @ gts) [] 
   284    val _ = if not (forall fst rawdests) then error "real_vector_combo_prover: Sanity check" 
   285            else ()
   286    val dests = distinct (op aconvc) (map snd rawdests)
   287    val srcfuns = map vector_lincomb sources
   288    val destfuns = map vector_lincomb dests 
   289    val vvs = fold_rev (curry (gen_union op aconvc) o Ctermfunc.dom) (srcfuns @ destfuns) []
   290    val n = length srcfuns
   291    val nvs = 1 upto n
   292    val srccombs = srcfuns ~~ nvs
   293    fun consider d =
   294     let 
   295      fun coefficients x =
   296       let 
   297        val inp = if Ctermfunc.defined d x then Intfunc.onefunc (0, Rat.neg(Ctermfunc.apply d x))
   298                       else Intfunc.undefined 
   299       in fold_rev (fn (f,v) => fn g => if Ctermfunc.defined f x then Intfunc.update (v, Ctermfunc.apply f x) g else g) srccombs inp 
   300       end
   301      val equations = map coefficients vvs
   302      val inequalities = map (fn n => Intfunc.onefunc (n,Rat.one)) nvs
   303      fun plausiblevertices f =
   304       let 
   305        val flippedequations = map (fold_rev int_flip f) equations
   306        val constraints = flippedequations @ inequalities
   307        val rawverts = vertices nvs constraints
   308        fun check_solution v =
   309         let 
   310           val f = fold_rev2 (curry Intfunc.update) nvs v (Intfunc.onefunc (0, Rat.one))
   311         in forall (fn e => evaluate f e =/ Rat.zero) flippedequations
   312         end
   313        val goodverts = filter check_solution rawverts
   314        val signfixups = map (fn n => if n mem_int  f then ~1 else 1) nvs 
   315       in map (map2 (fn s => fn c => Rat.rat_of_int s */ c) signfixups) goodverts
   316       end
   317      val allverts = fold_rev append (map plausiblevertices (allsubsets nvs)) [] 
   318     in subsume allverts []
   319     end
   320    fun compute_ineq v =
   321     let 
   322      val ths = map_filter (fn (v,t) => if v =/ Rat.zero then NONE 
   323                                      else SOME(norm_cmul_rule v t))
   324                             (v ~~ nubs) 
   325      fun end_itlist f xs = split_last xs |> uncurry (fold_rev f)
   326     in inequality_canon_rule ctxt (end_itlist norm_add_rule ths)
   327     end
   328    val ges' = map_filter (try compute_ineq) (fold_rev (append o consider) destfuns []) @
   329                  map (inequality_canon_rule ctxt) nubs @ ges
   330    val zerodests = filter
   331         (fn t => null (Ctermfunc.dom (vector_lincomb t))) (map snd rawdests)
   332 
   333   in RealArith.real_linear_prover translator
   334         (map (fn t => instantiate ([(tv_n, ctyp_of_term t)],[]) pth_zero)
   335             zerodests,
   336         map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv
   337                        arg_conv (arg_conv real_poly_conv))) ges',
   338         map (fconv_rule (try_conv (More_Conv.top_sweep_conv (K norm_canon_conv) ctxt) then_conv 
   339                        arg_conv (arg_conv real_poly_conv))) gts)
   340   end
   341 in val real_vector_combo_prover = real_vector_combo_prover
   342 end;
   343 
   344 local
   345  val pth = @{thm norm_imp_pos_and_ge}
   346  val norm_mp = match_mp pth
   347  val concl = dest_arg o cprop_of
   348  fun conjunct1 th = th RS @{thm conjunct1}
   349  fun conjunct2 th = th RS @{thm conjunct2}
   350  fun C f x y = f y x
   351 fun real_vector_ineq_prover ctxt translator (ges,gts) = 
   352  let 
   353 (*   val _ = error "real_vector_ineq_prover: pause" *)
   354   val ntms = fold_rev find_normedterms (map (dest_arg o concl) (ges @ gts)) []
   355   val lctab = vector_lincombs (map snd (filter (not o fst) ntms))
   356   val (fxns, ctxt') = Variable.variant_fixes (replicate (length lctab) "x") ctxt
   357   fun instantiate_cterm' ty tms = Drule.cterm_rule (Drule.instantiate' ty tms)
   358   fun mk_norm t = capply (instantiate_cterm' [SOME (ctyp_of_term t)] [] @{cpat "norm :: (?'a :: real_normed_vector) => real"}) t
   359   fun mk_equals l r = capply (capply (instantiate_cterm' [SOME (ctyp_of_term l)] [] @{cpat "op == :: ?'a =>_"}) l) r
   360   val asl = map2 (fn (t,_) => fn n => assume (mk_equals (mk_norm t) (cterm_of (ProofContext.theory_of ctxt') (Free(n,@{typ real}))))) lctab fxns
   361   val replace_conv = try_conv (rewrs_conv asl)
   362   val replace_rule = fconv_rule (funpow 2 arg_conv (replacenegnorms replace_conv))
   363   val ges' =
   364        fold_rev (fn th => fn ths => conjunct1(norm_mp th)::ths)
   365               asl (map replace_rule ges)
   366   val gts' = map replace_rule gts
   367   val nubs = map (conjunct2 o norm_mp) asl
   368   val th1 = real_vector_combo_prover ctxt' translator (nubs,ges',gts')
   369   val shs = filter (member (fn (t,th) => t aconvc cprop_of th) asl) (#hyps (crep_thm th1)) 
   370   val th11 = hd (Variable.export ctxt' ctxt [fold implies_intr shs th1])
   371   val cps = map (swap o dest_equals) (cprems_of th11)
   372   val th12 = instantiate ([], cps) th11
   373   val th13 = fold (C implies_elim) (map (reflexive o snd) cps) th12;
   374  in hd (Variable.export ctxt' ctxt [th13])
   375  end 
   376 in val real_vector_ineq_prover = real_vector_ineq_prover
   377 end;
   378 
   379 local
   380  val rawrule = fconv_rule (arg_conv (rewr_conv @{thm real_eq_0_iff_le_ge_0}))
   381  fun conj_pair th = (th RS @{thm conjunct1}, th RS @{thm conjunct2})
   382  fun simple_cterm_ord t u = TermOrd.term_ord (term_of t, term_of u) = LESS;
   383   (* FIXME: Lookup in the context every time!!! Fix this !!!*)
   384  fun splitequation ctxt th acc =
   385   let 
   386    val real_poly_neg_conv = #neg
   387        (Normalizer.semiring_normalizers_ord_wrapper ctxt
   388         (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) simple_cterm_ord)
   389    val (th1,th2) = conj_pair(rawrule th)
   390   in th1::fconv_rule (arg_conv (arg_conv real_poly_neg_conv)) th2::acc
   391   end
   392 in fun real_vector_prover ctxt translator (eqs,ges,gts) =
   393      real_vector_ineq_prover ctxt translator
   394          (fold_rev (splitequation ctxt) eqs ges,gts)
   395 end;
   396 
   397   fun init_conv ctxt = 
   398    Simplifier.rewrite (Simplifier.context ctxt 
   399      (HOL_basic_ss addsimps ([(*@{thm vec_0}, @{thm vec_1},*) @{thm vector_dist_norm}, @{thm diff_0_right}, @{thm right_minus}, @{thm diff_self}, @{thm norm_zero}] @ @{thms arithmetic_simps} @ @{thms norm_pths})))
   400    then_conv field_comp_conv 
   401    then_conv nnf_conv
   402 
   403  fun pure ctxt = RealArith.gen_prover_real_arith ctxt (real_vector_prover ctxt);
   404  fun norm_arith ctxt ct = 
   405   let 
   406    val ctxt' = Variable.declare_term (term_of ct) ctxt
   407    val th = init_conv ctxt' ct
   408   in equal_elim (Drule.arg_cong_rule @{cterm Trueprop} (symmetric th)) 
   409                 (pure ctxt' (rhs_of th))
   410  end
   411 
   412  fun norm_arith_tac ctxt = 
   413    clarify_tac HOL_cs THEN'
   414    ObjectLogic.full_atomize_tac THEN'
   415    CSUBGOAL ( fn (p,i) => rtac (norm_arith ctxt (Thm.dest_arg p )) i);
   416 
   417 end;