added quickcheck[approximation]
authorimmler
Wed Nov 12 17:37:43 2014 +0100 (2014-11-12)
changeset 589886ebf918128b9
parent 58987 119680ebf37c
child 58989 99831590def5
added quickcheck[approximation]
src/HOL/Decision_Procs/Approximation.thy
src/HOL/Decision_Procs/Decision_Procs.thy
src/HOL/Decision_Procs/approximation_generator.ML
src/HOL/Decision_Procs/ex/Approximation_Quickcheck_Ex.thy
     1.1 --- a/src/HOL/Decision_Procs/Approximation.thy	Wed Nov 12 17:37:43 2014 +0100
     1.2 +++ b/src/HOL/Decision_Procs/Approximation.thy	Wed Nov 12 17:37:43 2014 +0100
     1.3 @@ -3494,7 +3494,8 @@
     1.4      | term_of_bool false = @{term False};
     1.5  
     1.6    val mk_int = HOLogic.mk_number @{typ int} o @{code integer_of_int};
     1.7 -  val dest_int = @{code int_of_integer} o snd o HOLogic.dest_number;
     1.8 +  fun dest_int (@{term int_of_integer} $ j) = @{code int_of_integer} (snd (HOLogic.dest_number j))
     1.9 +    | dest_int i = @{code int_of_integer} (snd (HOLogic.dest_number i));
    1.10  
    1.11    fun term_of_float (@{code Float} (k, l)) =
    1.12      @{term Float} $ mk_int k $ mk_int l;
    1.13 @@ -3706,4 +3707,11 @@
    1.14  
    1.15  ML_file "approximation.ML"
    1.16  
    1.17 +
    1.18 +section "Quickcheck Generator"
    1.19 +
    1.20 +ML_file "approximation_generator.ML"
    1.21 +
    1.22 +setup "Approximation_Generator.setup"
    1.23 +
    1.24  end
     2.1 --- a/src/HOL/Decision_Procs/Decision_Procs.thy	Wed Nov 12 17:37:43 2014 +0100
     2.2 +++ b/src/HOL/Decision_Procs/Decision_Procs.thy	Wed Nov 12 17:37:43 2014 +0100
     2.3 @@ -10,6 +10,7 @@
     2.4    Commutative_Ring_Complete
     2.5    "ex/Commutative_Ring_Ex"
     2.6    "ex/Approximation_Ex"
     2.7 +  "ex/Approximation_Quickcheck_Ex"
     2.8    "ex/Dense_Linear_Order_Ex"
     2.9  begin
    2.10  
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/src/HOL/Decision_Procs/approximation_generator.ML	Wed Nov 12 17:37:43 2014 +0100
     3.3 @@ -0,0 +1,211 @@
     3.4 +(*  Title:      HOL/Decision_Procs/approximation_generator.ML
     3.5 +    Author:     Fabian Immler, TU Muenchen
     3.6 +*)
     3.7 +
     3.8 +signature APPROXIMATION_GENERATOR =
     3.9 +sig
    3.10 +  val custom_seed: int Config.T
    3.11 +  val precision: int Config.T
    3.12 +  val epsilon: real Config.T
    3.13 +  val approximation_generator:
    3.14 +    Proof.context ->
    3.15 +    (term * term list) list ->
    3.16 +    bool -> int list -> (bool * term list) option * Quickcheck.report option
    3.17 +  val setup: theory -> theory
    3.18 +end;
    3.19 +
    3.20 +structure Approximation_Generator : APPROXIMATION_GENERATOR =
    3.21 +struct
    3.22 +
    3.23 +val custom_seed = Attrib.setup_config_int @{binding quickcheck_approximation_custom_seed} (K ~1)
    3.24 +
    3.25 +val precision = Attrib.setup_config_int @{binding quickcheck_approximation_precision} (K 30)
    3.26 +
    3.27 +val epsilon = Attrib.setup_config_real @{binding quickcheck_approximation_epsilon} (K 0.0)
    3.28 +
    3.29 +val random_float = @{code "random_class.random::_ \<Rightarrow> _ \<Rightarrow> (float \<times> (unit \<Rightarrow> term)) \<times> _"}
    3.30 +
    3.31 +fun nat_of_term t =
    3.32 +  (HOLogic.dest_nat t handle TERM _ => snd (HOLogic.dest_number t)
    3.33 +    handle TERM _ => raise TERM ("nat_of_term", [t]));
    3.34 +
    3.35 +fun int_of_term t = snd (HOLogic.dest_number t) handle TERM _ => raise TERM ("int_of_term", [t]);
    3.36 +
    3.37 +fun real_of_man_exp m e = Real.fromManExp {man = Real.fromInt m, exp = e}
    3.38 +
    3.39 +fun mapprox_float (@{term Float} $ m $ e) = real_of_man_exp (int_of_term m) (int_of_term e)
    3.40 +  | mapprox_float t = Real.fromInt (snd (HOLogic.dest_number t))
    3.41 +      handle TERM _ => raise TERM ("mapprox_float", [t]);
    3.42 +
    3.43 +(* TODO: define using compiled terms? *)
    3.44 +fun mapprox_floatarith (@{term Add} $ a $ b) xs = mapprox_floatarith a xs + mapprox_floatarith b xs
    3.45 +  | mapprox_floatarith (@{term Minus} $ a) xs = ~ (mapprox_floatarith a xs)
    3.46 +  | mapprox_floatarith (@{term Mult} $ a $ b) xs = mapprox_floatarith a xs * mapprox_floatarith b xs
    3.47 +  | mapprox_floatarith (@{term Inverse} $ a) xs = 1.0 / mapprox_floatarith a xs
    3.48 +  | mapprox_floatarith (@{term Cos} $ a) xs = Math.cos (mapprox_floatarith a xs)
    3.49 +  | mapprox_floatarith (@{term Arctan} $ a) xs = Math.atan (mapprox_floatarith a xs)
    3.50 +  | mapprox_floatarith (@{term Abs} $ a) xs = abs (mapprox_floatarith a xs)
    3.51 +  | mapprox_floatarith (@{term Max} $ a $ b) xs =
    3.52 +      Real.max (mapprox_floatarith a xs, mapprox_floatarith b xs)
    3.53 +  | mapprox_floatarith (@{term Min} $ a $ b) xs =
    3.54 +      Real.min (mapprox_floatarith a xs, mapprox_floatarith b xs)
    3.55 +  | mapprox_floatarith @{term Pi} _ = Math.pi
    3.56 +  | mapprox_floatarith (@{term Sqrt} $ a) xs = Math.sqrt (mapprox_floatarith a xs)
    3.57 +  | mapprox_floatarith (@{term Exp} $ a) xs = Math.exp (mapprox_floatarith a xs)
    3.58 +  | mapprox_floatarith (@{term Ln} $ a) xs = Math.ln (mapprox_floatarith a xs)
    3.59 +  | mapprox_floatarith (@{term Power} $ a $ n) xs =
    3.60 +      Math.pow (mapprox_floatarith a xs, Real.fromInt (nat_of_term n))
    3.61 +  | mapprox_floatarith (@{term Var} $ n) xs = nth xs (nat_of_term n)
    3.62 +  | mapprox_floatarith (@{term Num} $ m) _ = mapprox_float m
    3.63 +  | mapprox_floatarith t _ = raise TERM ("mapprox_floatarith", [t])
    3.64 +
    3.65 +fun mapprox_atLeastAtMost eps x a b xs =
    3.66 +    let
    3.67 +      val x' = mapprox_floatarith x xs
    3.68 +    in
    3.69 +      mapprox_floatarith a xs + eps <= x' andalso x' + eps <= mapprox_floatarith b xs
    3.70 +    end
    3.71 +
    3.72 +fun mapprox_form eps (@{term Bound} $ x $ a $ b $ f) xs =
    3.73 +    (not (mapprox_atLeastAtMost eps x a b xs)) orelse mapprox_form eps f xs
    3.74 +| mapprox_form eps (@{term Assign} $ x $ a $ f) xs =
    3.75 +    (Real.!= (mapprox_floatarith x xs, mapprox_floatarith a xs)) orelse mapprox_form eps f xs
    3.76 +| mapprox_form eps (@{term Less} $ a $ b) xs = mapprox_floatarith a xs + eps < mapprox_floatarith b xs
    3.77 +| mapprox_form eps (@{term LessEqual} $ a $ b) xs = mapprox_floatarith a xs + eps <= mapprox_floatarith b xs
    3.78 +| mapprox_form eps (@{term AtLeastAtMost} $ x $ a $ b) xs = mapprox_atLeastAtMost eps x a b xs
    3.79 +| mapprox_form eps (@{term Conj} $ f $ g) xs = mapprox_form eps f xs andalso mapprox_form eps g xs
    3.80 +| mapprox_form eps (@{term Disj} $ f $ g) xs = mapprox_form eps f xs orelse mapprox_form eps g xs
    3.81 +| mapprox_form _ t _ = raise TERM ("mapprox_form", [t])
    3.82 +
    3.83 +fun dest_interpret_form (@{const "interpret_form"} $ b $ xs) = (b, xs)
    3.84 +  | dest_interpret_form t = raise TERM ("dest_interpret_form", [t])
    3.85 +
    3.86 +fun optionT t = Type (@{type_name "option"}, [t])
    3.87 +fun mk_Some t = Const (@{const_name "Some"}, t --> optionT t)
    3.88 +
    3.89 +fun random_float_list size xs seed =
    3.90 +  fold (K (apsnd (random_float size) #-> (fn c => apfst (fn b => b::c)))) xs ([],seed)
    3.91 +
    3.92 +fun real_of_Float (@{code Float} (m, e)) =
    3.93 +    real_of_man_exp (@{code integer_of_int} m) (@{code integer_of_int} e)
    3.94 +
    3.95 +fun is_True @{term True} = true
    3.96 +  | is_True _ = false
    3.97 +
    3.98 +val postproc_form_eqs =
    3.99 +  @{lemma
   3.100 +    "real (Float 0 a) = 0"
   3.101 +    "real (Float (numeral m) 0) = numeral m"
   3.102 +    "real (Float 1 0) = 1"
   3.103 +    "real (Float (- 1) 0) = - 1"
   3.104 +    "real (Float 1 (numeral e)) = 2 ^ numeral e"
   3.105 +    "real (Float 1 (- numeral e)) = 1 / 2 ^ numeral e"
   3.106 +    "real (Float a 1) = a * 2"
   3.107 +    "real (Float a (-1)) = a / 2"
   3.108 +    "real (Float (- a) b) = - real (Float a b)"
   3.109 +    "real (Float (numeral m) (numeral e)) = numeral m * 2 ^ (numeral e)"
   3.110 +    "real (Float (numeral m) (- numeral e)) = numeral m / 2 ^ (numeral e)"
   3.111 +    "- (c * d::real) = -c * d"
   3.112 +    "- (c / d::real) = -c / d"
   3.113 +    "- (0::real) = 0"
   3.114 +    "int_of_integer (numeral k) = numeral k"
   3.115 +    "int_of_integer (- numeral k) = - numeral k"
   3.116 +    "int_of_integer 0 = 0"
   3.117 +    "int_of_integer 1 = 1"
   3.118 +    "int_of_integer (- 1) = - 1"
   3.119 +    by auto
   3.120 +  }
   3.121 +
   3.122 +fun rewrite_with ctxt thms = Simplifier.rewrite (put_simpset HOL_basic_ss ctxt addsimps thms)
   3.123 +fun conv_term thy conv r = cterm_of thy r |> conv |> Thm.prop_of |> Logic.dest_equals |> snd
   3.124 +
   3.125 +fun approx_random ctxt prec eps frees e xs genuine_only size seed =
   3.126 +  let
   3.127 +    val (rs, seed') = random_float_list size xs seed
   3.128 +    fun mk_approx_form e ts =
   3.129 +      @{const "approx_form"} $
   3.130 +        HOLogic.mk_number @{typ nat} prec $
   3.131 +        e $
   3.132 +        (HOLogic.mk_list @{typ "(float * float) option"}
   3.133 +          (map (fn t => mk_Some @{typ "float * float"} $ HOLogic.mk_prod (t, t)) ts)) $
   3.134 +        @{term "[] :: nat list"}
   3.135 +  in
   3.136 +    (if mapprox_form eps e (map (real_of_Float o fst) rs)
   3.137 +    then
   3.138 +      let
   3.139 +        val ts = map (fn x => snd x ()) rs
   3.140 +        val ts' = map
   3.141 +          (AList.lookup op = (map dest_Free xs ~~ ts)
   3.142 +            #> the_default Term.dummy
   3.143 +            #> curry op $ @{term "real::float\<Rightarrow>_"}
   3.144 +            #> conv_term (Proof_Context.theory_of ctxt) (rewrite_with ctxt postproc_form_eqs))
   3.145 +          frees
   3.146 +      in
   3.147 +        if approximate ctxt (mk_approx_form e ts) |> is_True
   3.148 +        then SOME (true, ts')
   3.149 +        else (if genuine_only then NONE else SOME (false, ts'))
   3.150 +      end
   3.151 +    else NONE, seed')
   3.152 +  end
   3.153 +
   3.154 +val preproc_form_eqs =
   3.155 +  @{lemma
   3.156 +    "(a::real) \<in> {b .. c} \<longleftrightarrow> b \<le> a \<and> a \<le> c"
   3.157 +    "a = b \<longleftrightarrow> a \<le> b \<and> b \<le> a"
   3.158 +    "(p \<longrightarrow> q) \<longleftrightarrow> \<not>p \<or> q"
   3.159 +    "(p \<longleftrightarrow> q) \<longleftrightarrow> (p \<longrightarrow> q) \<and> (q \<longrightarrow> p)"
   3.160 +    "\<not> (a < b) \<longleftrightarrow> b \<le> a"
   3.161 +    "\<not> (a \<le> b) \<longleftrightarrow> b < a"
   3.162 +    "\<not> (p \<and> q) \<longleftrightarrow> \<not> p \<or> \<not> q"
   3.163 +    "\<not> (p \<or> q) \<longleftrightarrow> \<not> p \<and> \<not> q"
   3.164 +    "\<not> \<not> q \<longleftrightarrow> q"
   3.165 +    by auto
   3.166 +  }
   3.167 +
   3.168 +fun reify_goal ctxt t =
   3.169 +  HOLogic.mk_not t
   3.170 +    |> conv_term (Proof_Context.theory_of ctxt) (rewrite_with ctxt preproc_form_eqs)
   3.171 +    |> conv_term (Proof_Context.theory_of ctxt) (Reification.conv ctxt form_equations)
   3.172 +    |> dest_interpret_form
   3.173 +    ||> HOLogic.dest_list
   3.174 +
   3.175 +fun approximation_generator_raw ctxt t =
   3.176 +  let
   3.177 +    val iterations = Config.get ctxt Quickcheck.iterations
   3.178 +    val prec = Config.get ctxt precision
   3.179 +    val eps = Config.get ctxt epsilon
   3.180 +    val cs = Config.get ctxt custom_seed
   3.181 +    val seed = (Code_Numeral.natural_of_integer (cs + 1), Code_Numeral.natural_of_integer 1)
   3.182 +    val run = if cs < 0
   3.183 +      then (fn f => fn seed => (Random_Engine.run f, seed))
   3.184 +      else (fn f => fn seed => f seed)
   3.185 +    val frees = Term.add_frees t []
   3.186 +    val (e, xs) = reify_goal ctxt t
   3.187 +    fun single_tester b s =
   3.188 +      approx_random ctxt prec eps frees e xs b s |> run
   3.189 +    fun iterate _ _ 0 _ = NONE
   3.190 +      | iterate genuine_only size j seed =
   3.191 +        case single_tester genuine_only size seed of
   3.192 +          (NONE, seed') => iterate genuine_only size (j - 1) seed'
   3.193 +        | (SOME q, _) => SOME q
   3.194 +  in
   3.195 +    fn genuine_only => fn size => (iterate genuine_only size iterations seed, NONE)
   3.196 +  end
   3.197 +
   3.198 +fun approximation_generator ctxt [(t, _)] =
   3.199 +  (fn genuine_only =>
   3.200 +    fn [_, size] =>
   3.201 +      approximation_generator_raw ctxt t genuine_only
   3.202 +        (Code_Numeral.natural_of_integer size))
   3.203 +  | approximation_generator _ _ =
   3.204 +      error "Quickcheck-approximation does not support type variables (or finite instantiations)"
   3.205 +
   3.206 +val test_goals =
   3.207 +  Quickcheck_Common.generator_test_goal_terms
   3.208 +    ("approximation", (fn _ => fn _ => false, approximation_generator))
   3.209 +
   3.210 +val active = Attrib.setup_config_bool @{binding quickcheck_approximation_active} (K false)
   3.211 +
   3.212 +val setup = Context.theory_map (Quickcheck.add_tester ("approximation", (active, test_goals)))
   3.213 +
   3.214 +end
     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/src/HOL/Decision_Procs/ex/Approximation_Quickcheck_Ex.thy	Wed Nov 12 17:37:43 2014 +0100
     4.3 @@ -0,0 +1,39 @@
     4.4 +theory Approximation_Quickcheck_Ex
     4.5 +imports "../Approximation"
     4.6 +begin
     4.7 +
     4.8 +lemma
     4.9 +  fixes x::real and y::real
    4.10 +  shows "sin x \<le> tan x"
    4.11 +  using [[quickcheck_approximation_custom_seed = 1]]
    4.12 +  quickcheck[approximation, expect=counterexample]
    4.13 +  oops
    4.14 +
    4.15 +lemma "x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x"
    4.16 +  using [[quickcheck_approximation_custom_seed = 1]]
    4.17 +  quickcheck[approximation, expect=counterexample]
    4.18 +  oops
    4.19 +
    4.20 +lemma "0 < x \<Longrightarrow> x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x"
    4.21 +  using [[quickcheck_approximation_custom_seed = 1]]
    4.22 +  quickcheck[approximation, expect=no_counterexample]
    4.23 +  by (rule arctan_divide_mono)
    4.24 +
    4.25 +lemma
    4.26 +  fixes x::real
    4.27 +  shows "exp (exp x + exp y + sin x * sin y) - 0.4 > 0 \<or> 0.98 - sin x / (sin x * sin y + 2)^2 > 0"
    4.28 +  using [[quickcheck_approximation_custom_seed = 1]]
    4.29 +  quickcheck[approximation, expect=counterexample, size=10, iterations=1000, verbose]
    4.30 +  oops
    4.31 +
    4.32 +lemma
    4.33 +  fixes x::real
    4.34 +  shows "x > 1 \<Longrightarrow> x \<le> 2 powr 20 * log 2 x + 1 \<and> (sin x)\<^sup>2 + (cos x)\<^sup>2 = 1"
    4.35 +  using [[quickcheck_approximation_custom_seed = 1]]
    4.36 +  using [[quickcheck_approximation_epsilon = 0.00000001]]
    4.37 +    --\<open>avoids spurious counterexamples in approximate computation of @{term "(sin x)\<^sup>2 + (cos x)\<^sup>2"}
    4.38 +      and therefore avoids expensive failing attempts for certification\<close>
    4.39 +  quickcheck[approximation, expect=counterexample, size=20]
    4.40 +  oops
    4.41 +
    4.42 +end