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