formal_lp/hypermap/ineqs/lp_approx_ineqs.hl

   1 needs "../formal_lp/ineqs/constants_approx.hl";;
   2
   3 module Lp_approx_ineqs = struct
   4
   5 open Constants_approx;;
   6
   7 (* Integral approximation of inequalities and equalities *)
   8 let mk_decimal =
   9   let decimal_const = `DECIMAL` and
  10       neg_op_real = `(--):real->real` in
  11     fun (n,m) ->
  12       let tm = mk_comb(mk_comb(decimal_const, mk_numeral (abs_num n)), mk_numeral m) in
  13         if (n </ Int 0) then mk_comb(neg_op_real, tm) else tm;;
  14
  15
  16 (*
  17 Given a rational number term `r` and a natural number k,
  18 returns a pair of theorems:
  19 |- low <= r, |- r <= high,
  20 with decimal low and high such that
  21 |low - r| <= 10^(-k) and |high - r| <= 10^(-k)
  22 *)
  23 let real_rat_approx =
  24   let le_op_real = `(<=):real->real->bool` in
  25     fun tm precision ->
  26       let n = rat_of_term tm in
  27       let m = Int 10 **/ (Int precision) in
  28       let n1 = n */ m in
  29       let low, high = floor_num n1, ceiling_num n1 in
  30       let l_tm, h_tm =
  31         if precision = 0 then
  32           term_of_rat low, term_of_rat high
  33         else
  34           mk_decimal (low, m), mk_decimal (high, m) in
  35       let l_th = EQT_ELIM (REAL_RAT_LE_CONV (mk_binop le_op_real l_tm tm)) in
  36       let h_th = EQT_ELIM (REAL_RAT_LE_CONV (mk_binop le_op_real tm h_tm)) in
  37         l_th, h_th;;
  38
  39
  40 (* Splits a sum into two parts: with and without free variables *)
  41 let split_sum_conv =
  42   let sum_th0 = REAL_ARITH `x = x + &0 + &0` and
  43       sum_th1 = REAL_ARITH `(x + y) + b + c = y + (x + b) + c:real` and
  44       sum_th2 = REAL_ARITH `(x + y) + b + c = y + b + (x + c):real` and
  45       sum_th1' = REAL_ARITH `x + b + c = (x + b) + c:real` and
  46       sum_th2' = REAL_ARITH `x + b + c = b + (x + c):real` in
  47   let x_var_real = `x:real` and
  48       y_var_real = `y:real` and
  49       b_var_real = `b:real` and
  50       c_var_real = `c:real` and
  51       add_op_real = `(+):real->real->real` in
  52   let rec split_sum_raw_conv tm =
  53     let xy_tm, bc_tm = dest_binop add_op_real tm in
  54     let b_tm, c_tm = dest_binop add_op_real bc_tm in
  55       if (is_binop add_op_real xy_tm) then
  56         let x_tm, y_tm = dest_binop add_op_real xy_tm in
  57         let inst_th = INST[x_tm, x_var_real; y_tm, y_var_real; b_tm, b_var_real; c_tm, c_var_real] in
  58         let th1 = if (frees x_tm <> []) then inst_th sum_th1 else inst_th sum_th2 in
  59         let th2 = split_sum_raw_conv (rand(concl th1)) in
  60           TRANS th1 th2
  61       else
  62         let inst_th = INST[xy_tm, x_var_real; b_tm, b_var_real; c_tm, c_var_real] in
  63           if (frees xy_tm <> []) then inst_th sum_th1' else inst_th sum_th2' in
  64
  65     fun tm ->
  66       let th0 = INST[tm, x_var_real] sum_th0 in
  67       let th1 = split_sum_raw_conv (rand(concl th0)) in
  68         TRANS th0 th1;;
  69
  70
  71
  72 let rearrange_mul_conv =
  73   let mul_op_real = `( * ):real->real->real` in
  74   let rec dest_mul tm =
  75     if (is_binop mul_op_real tm) then
  76       let lhs, rhs = dest_binop mul_op_real tm in
  77       let cs, vars = dest_mul rhs in
  78         if (frees lhs = []) then
  79           lhs::cs, vars
  80         else
  81           cs, lhs::vars
  82     else
  83       if (frees tm = []) then
  84         [tm], []
  85       else
  86         [], [tm] in
  87   let mk_mul list =
  88     if (list = []) then failwith "rearrange_mul: empty list"
  89     else itlist (fun l r -> mk_binop mul_op_real l r) (tl list) (hd list) in
  90
  91     fun tm ->
  92       let cs, vars = dest_mul tm in
  93       let cs_mul, vars_mul = mk_mul cs, mk_mul vars in
  94       let t = mk_eq(tm, mk_binop mul_op_real cs_mul vars_mul) in
  95         EQT_ELIM (REWRITE_CONV[REAL_MUL_AC] t);;
  96
  97
  98
  99 (* Moves everything with free variables on the left and performs basic reductions *)
 100 let ineq_rewrite_conv =
 101   let le_add_th = REAL_ARITH `x + y <= b + c <=> x - b <= c - y:real` in
 102     REWRITE_CONV[real_ge; real_div; DECIMAL] THENC
 103       LAND_CONV REAL_POLY_CONV THENC RAND_CONV REAL_POLY_CONV THENC
 104       LAND_CONV split_sum_conv THENC RAND_CONV split_sum_conv THENC
 105       ONCE_REWRITE_CONV[le_add_th] THENC
 106       LAND_CONV REAL_POLY_CONV THENC RAND_CONV REAL_POLY_CONV THENC
 107       REWRITE_CONV[GSYM real_div] THENC
 108       ONCE_DEPTH_CONV rearrange_mul_conv;;
 109
 110
 111 (* Approximation *)
 112
 113 let le_mul1_th = REAL_ARITH `!x. &1 * x <= x`;;
 114 let ge_mul1_th = REAL_ARITH `!x. x <= &1 * x`;;
 115 let INTERVAL_LO = prove(`interval_arith x (a,b) ==> a <= x`, SIMP_TAC[interval_arith]);;
 116 let INTERVAL_HI = prove(`interval_arith x (a,b) ==> x <= b`, SIMP_TAC[interval_arith]);;
 117
 118 let add_op_real = `(+):real->real->real`;;
 119
 120 let rec low_approx tm precision =
 121   let low_approx1 tm precision =
 122     if (is_binop mul_op_real tm && frees tm <> []) then
 123       let c, var = dest_binop mul_op_real tm in
 124       let interval_th = approx_interval (create_interval c) precision in
 125       let l_th = MATCH_MP INTERVAL_LO interval_th in
 126       let a, b = dest_binop le_op_real (concl l_th) in
 127         (PROVE_HYP l_th o UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPECL[a; b; var]) REAL_LE_RMUL
 128     else
 129       if (frees tm = []) then
 130         MATCH_MP INTERVAL_LO (approx_interval (create_interval tm) precision)
 131       else
 132         SPEC tm le_mul1_th in
 133
 134   if (is_binop add_op_real tm && frees tm <> []) then
 135     let lhs, rhs = dest_binop add_op_real tm in
 136     let l_th = low_approx1 lhs precision in
 137     let r_th = low_approx rhs precision in
 138       MATCH_MP REAL_LE_ADD2 (CONJ l_th r_th)
 139   else
 140     low_approx1 tm precision;;
 141
 142
 143 let rec hi_approx tm precision =
 144   let hi_approx1 tm precision =
 145     if (is_binop mul_op_real tm && frees tm <> []) then
 146       let c, var = dest_binop mul_op_real tm in
 147       let interval_th = approx_interval (create_interval c) precision in
 148       let h_th = MATCH_MP INTERVAL_HI interval_th in
 149       let a, b = dest_binop le_op_real (concl h_th) in
 150         (PROVE_HYP h_th o UNDISCH_ALL o REWRITE_RULE[GSYM IMP_IMP] o SPECL[a; b; var]) REAL_LE_RMUL
 151     else
 152       if (frees tm = []) then
 153         MATCH_MP INTERVAL_HI (approx_interval (create_interval tm) precision)
 154       else
 155         SPEC tm ge_mul1_th in
 156
 157   if (is_binop add_op_real tm && frees tm <> []) then
 158     let lhs, rhs = dest_binop add_op_real tm in
 159     let l_th = hi_approx1 lhs precision in
 160     let r_th = hi_approx rhs precision in
 161       MATCH_MP REAL_LE_ADD2 (CONJ l_th r_th)
 162   else
 163     hi_approx1 tm precision;;
 164
 165
 166
 167 let approx_le_ineq precision tm =
 168   let lhs, rhs = dest_binop le_op_real tm in
 169   let lhs_th = low_approx lhs precision in
 170   let rhs_th = hi_approx rhs precision in
 171   let ll, lr = dest_binop le_op_real (concl lhs_th) in
 172   let rl, rr = dest_binop le_op_real (concl rhs_th) in
 173
 174   let th0 = ASSUME tm in
 175   let th1 = SPECL[ll; lr; rhs] REAL_LE_TRANS in
 176   let th2 = SPECL[ll; rhs; rr] REAL_LE_TRANS in
 177   let s1 = MATCH_MP th1 (CONJ lhs_th th0) in
 178     MATCH_MP th2 (CONJ s1 rhs_th);;
 179
 180
 181 let integer_approx_le_ineq precision ineq =
 182   let lhs, rhs = dest_binop le_op_real ineq in
 183   let m = (Int 10 **/ Int precision) in
 184   let m_num, m_real = mk_numeral m, term_of_rat m in
 185   let m_pos = SPEC m_num REAL_POS in
 186   let mul_th = SPECL[m_real; lhs; rhs] REAL_LE_LMUL in
 187   let th0 = MATCH_MP mul_th (CONJ m_pos (ASSUME ineq)) in
 188   let th1 = (CONV_RULE ineq_rewrite_conv) th0 in
 189   let approx = approx_le_ineq 0 (concl th1) in
 190     PROVE_HYP th1 approx;;
 191
 192
 193 let create_approximations precision_list ineq =
 194   let ineq_th = ineq_rewrite_conv ineq in
 195   let rhs = rand(concl ineq_th) in
 196   let th0 = CONV_RULE (LAND_CONV (REWRITE_CONV[DECIMAL] THENC REAL_POLY_CONV))
 197     (approx_le_ineq 8 rhs) in
 198
 199   let int_approx p =
 200     let th1 = integer_approx_le_ineq p (concl th0) in
 201     let th2 = PROVE_HYP th0 th1 in
 202     let th3 = DISCH rhs th2 in
 203       REWRITE_RULE[GSYM ineq_th] th3 in
 204
 205     map int_approx precision_list;;
 206
 207 (**********************************)
 208
 209
 210 (*
 211 let tm = `(&34/ &13 + pi/sol0)* x + &2 + -- &14 / &3 * z <= pi + rho218 + z * sol0 - u / &1000`;;
 212 create_approximations [3;4;5] tm;;
 213 *)
 214
 215
 216 (*************************)
 217
 218 (* Additional step for generalizing hypotheses *)
 219
 220 let LIST_SUM_LMUL = prove(`!(f:A->real) c n. list_sum n (\x. c * f x) = c * list_sum n f`,
 221                           REWRITE_TAC[Seq2.list_sum_lmul]);;
 222 let le_hyp_gen = prove(`!f y. (!x. &0 <= f x) ==> &0 <= f y`, SIMP_TAC[]);;
 223 let le_list_sum_hyp_gen = prove(`!(f:A->real) n. (!x. &0 <= f x) ==> &0 <= list_sum n f`,
 224                                 REPEAT STRIP_TAC THEN MATCH_MP_TAC Seq2.list_sum_ge0 THEN ASM_REWRITE_TAC[]);;
 225
 226
 227 let generalize_hyp th =
 228   let gen_hyp tm =
 229     let fn, arg = dest_comb (rand tm) in
 230       if (is_comb fn && is_const (rator fn) && (fst o dest_const o rator) fn = "list_sum") then
 231         let f, n = arg, rand fn in
 232           UNDISCH_ALL (ISPECL[f; n] le_list_sum_hyp_gen)
 233       else
 234         UNDISCH_ALL (ISPECL[fn; arg] le_hyp_gen) in
 235
 236   let hyp_ths = map gen_hyp (hyp th) in
 237     itlist PROVE_HYP hyp_ths th;;
 238
 239
 240 (*
 241 let tm = `list_sum n (\x. s x * r) + &1 / &3 * (azim_dart (V,E) o g) x <= pi`;;
 242 let ths = create_approximations [3] tm;;
 243 map generalize_hyp ths;;
 244 *)
 245
 246
 247 (*******************************)
 248
 249
 250 let generate_ineqs =
 251   let imp_th = TAUT `(A ==> B) ==> ((P ==> A) ==> (P ==> B))` in
 252   let and_imp_th = TAUT `(A ==> B /\ C) <=> ((A ==> B) /\ (A ==> C))` in
 253   let p_bool = `P:bool` in
 254   let strip_imp = splitlist dest_imp in
 255   let list_mk_imp = itlist (curry mk_imp) in
 256
 257   let create_approxs pos_ths original_th precision_list =
 258     let var, ineq' = (dest_abs o rand o rator o concl) original_th in
 259     let cond_tms, ineq = strip_imp ineq' in
 260     let final_step = fun approx_th ->
 261       let p_tm', q_tm' = dest_imp (concl approx_th) in
 262       let p_tm, q_tm = mk_abs(var, list_mk_imp cond_tms p_tm'), mk_abs(var, list_mk_imp cond_tms q_tm') in
 263       let list_tm = (rand o concl) original_th in
 264       let mono_th = BETA_RULE (ISPECL[p_tm; q_tm; list_tm] (GEN_ALL MONO_ALL)) in
 265       let approx_th2 = itlist (fun p_tm th -> MATCH_MP (INST [p_tm, p_bool] imp_th) th) cond_tms approx_th in
 266       let s1 = MATCH_MP mono_th (GEN var approx_th2) in
 267         MATCH_MP s1 original_th in
 268
 269     let approx_ths0 = map generalize_hyp (create_approximations precision_list ineq) in
 270     let approx_ths1 = map (itlist PROVE_HYP pos_ths) approx_ths0 in
 271     let approx_ths = map (REWRITE_RULE[GSYM LIST_SUM_LMUL]) approx_ths1 in
 272       map final_step approx_ths in
 273
 274     fun pos_ths precision_list ineq_th ->
 275       let _, ineq = (strip_imp o snd o dest_abs o rand o rator o concl) ineq_th in
 276         if (is_eq ineq) then
 277           let eq_th = PURE_REWRITE_RULE[GSYM REAL_LE_ANTISYM; and_imp_th; GSYM AND_ALL] ineq_th in
 278           let ths1 = create_approxs pos_ths (CONJUNCT1 eq_th) precision_list in
 279           let ths2 = create_approxs pos_ths (CONJUNCT2 eq_th) precision_list in
 280             ths1, ths2
 281         else
 282           create_approxs pos_ths ineq_th precision_list, [];;
 283
 284
 285
 286 (*
 287 generate_ineqs [] [3;4] (ASSUME `ALL (\n. list_sum n (\x. f x) - &1 <= pi) s`);;
 288 *)
 289
 290 end;;
 291
 292