Commit fd8844a1 authored by Heiko Becker's avatar Heiko Becker
Browse files

Rework interval arithmetic to make it more understandable, still needs the multiplication result

parent 2366a153
......@@ -2,6 +2,7 @@
Infra/abbrevs.v
exps.v
daisy_lang.v
newIntervalArith.v
interval_arith.v
abs_err.v
simple_doppler.v
......@@ -3,41 +3,31 @@
used to verify analsysis result in the final theorem of a certificate.
**)
Require Import Coq.Reals.Reals.
Require Import Daisy.daisy_lang Daisy.exps Daisy.Infra.abbrevs Daisy.interval_arith.
Require Import Daisy.daisy_lang Daisy.exps Daisy.Infra.abbrevs Daisy.newIntervalArith.
Definition abs_env:Type := exp R -> intv -> err -> Prop.
(*
Definition ErrFromEnv (env:abs_env) (e:exp R) := (env e).
Definition IntvFromEnv (env:abs_env) (e:exp R) := fst (env e).
Arguments ErrFromEnv _ _ /.
Arguments IntvFromEnv _ _ /.
*)
(**
Define absolute error of an expression as an inductive relation
**)
First define when given an interval and an error, when another error is sound overapproximation of the absolute error
**)
Definition isSoundErr (error:err) (iv:intv) (propagatedErr:err) :=
(error >= m_eps * Rmax (Rabs (IVlo(iv) - propagatedErr)) (Rabs (IVhi(iv) + propagatedErr)))%R.
Definition isSupersetIntv (res:intv) (comp:intv) :=
(IVhi res >= IVhi comp)%R /\ (IVlo res <= IVlo comp)%R.
Definition maxAbs (IV:intv) := Rmax (Rabs (IVlo IV)) (Rabs (IVhi IV)).
Definition mult_err (e1:exp R) (val_e1:ann) (e2:exp R) (val_e2:ann) :=
(maxAbs
(intv_add
(intv_add
(intv_mult
(addInterval
(addInterval
(multInterval
(getIntv val_e1)
(getErr val_e2, getErr val_e2))
(intv_mult
(multInterval
(getIntv val_e2)
(getErr val_e1, getErr val_e1)))
(intv_mult
(intv_mult
(intv_mult
(multInterval
(multInterval
(multInterval
(getIntv val_e1)
(getIntv val_e2))
(getErr val_e1, getErr val_e1))
......@@ -46,16 +36,16 @@ Definition mult_err (e1:exp R) (val_e1:ann) (e2:exp R) (val_e2:ann) :=
Definition isSoundErrBin (op:binop) (e1:exp R) (val_e1:ann) (e2:exp R) (val_e2:ann) (val:ann):=
match op with
Plus
=> (isSupersetIntv (getIntv val) (intv_add_err (getIntv val_e1) (getErr val_e1) (getIntv val_e2) (getErr val_e2)) /\
=> (isSupersetInterval (getIntv val) (addInterval (widenInterval (getIntv val_e1) (getErr val_e1)) (widenInterval (getIntv val_e2) (getErr val_e2))) /\
isSoundErr (getErr val) (getIntv val) (getErr val_e1 + getErr val_e2)%R)
|Sub
=> (isSupersetIntv (getIntv val) (intv_sub_err (getIntv val_e1) (getErr val_e1) (getIntv val_e2) (getErr val_e2)) /\
=> (isSupersetInterval (getIntv val) (substractInterval (widenInterval (getIntv val_e1) (getErr val_e1)) (widenInterval (getIntv val_e2) (getErr val_e2))) /\
isSoundErr (getErr val) (getIntv val) (getErr val_e1 - getErr val_e2)%R)
|Mult
=> (isSupersetIntv (getIntv val) (intv_mult_err (getIntv val_e1) (getErr val_e1) (getIntv val_e2) (getErr val_e2)) /\
=> (isSupersetInterval (getIntv val) (multInterval (widenInterval (getIntv val_e1) (getErr val_e1)) (widenInterval (getIntv val_e2) (getErr val_e2))) /\
isSoundErr (getErr val) (getIntv val) (mult_err e1 val_e1 e2 val_e2))
|Div
=> (isSupersetIntv (getIntv val) (intv_div_err (getIntv val_e1) (getErr val_e1) (getIntv val_e2) (getErr val_e2)) /\
=> (isSupersetInterval (getIntv val) (intv_div_err (getIntv val_e1) (getErr val_e1) (getIntv val_e2) (getErr val_e2)) /\
isSoundErr (getErr val) (getIntv val) (mult_err e1 val_e1 e2 (((/ IVlo (getIntv val_e2))%R, (/ IVhi (getIntv val_e2))%R), getErr val_e2)))
end.
......
(**
Formalization of real valued interval arithmetic
TODO: Package this into a class or module that depends only on proving the existence of basic operators instead
**)
Require Import Coq.Reals.Reals Coq.micromega.Psatz.
Require Import Daisy.Infra.abbrevs.
(**
Intervals are a type, consisting of a pair of two real numbers
Additionally add some constructing and destructing definitions for encapsulation and define them to automatically unfold upon simplification.
**)
Definition interval:Type := R * R.
Definition mkInterval (ivlo:R) (ivhi:R) := (ivlo,ivhi).
Definition IVlo (intv:interval) := fst intv.
Definition IVhi (intv:interval) := snd intv.
Arguments mkInterval _ _/.
Arguments IVlo _ /.
Arguments IVhi _ /.
(**
Define validity of an interval, requiring that the lower bound is less than or equal to the upper bound.
Containement is defined such that if x is contained in the interval, it must lie between the lower and upper bound.
**)
Definition valid (intv:interval) :Prop :=
(IVlo intv <= IVhi intv)%R.
Definition contained (x:R) (intv:interval) :=
(IVlo intv <= x <= IVhi intv)%R.
Lemma containedImpliesValid (x:R) (intv:interval) :
contained x intv ->
valid intv.
Proof.
unfold contained, valid; intros inIntv; apply (Rle_trans _ x _); destruct inIntv; auto.
Qed.
(**
Subset definition; when an interval is a subinterval of another
**)
Definition isSupersetInterval (iv1:interval) (iv2:interval) :=
(IVlo iv2 <= IVlo iv1)%R /\ (IVhi iv1 <= IVhi iv2)%R.
Definition pointInterval (x:R) := mkInterval x x.
Lemma containedImpliesSubset (x:R) (intv:interval):
contained x intv ->
isSupersetInterval (pointInterval x) intv.
Proof.
intros containedIntv.
unfold contained, isSupersetInterval, pointInterval in *; destruct containedIntv; split; auto.
Qed.
(**
Definition of validity conditions for interval operations, Addition, Subtraction, Multiplication and Division
**)
Definition validIntervalAdd (iv1:interval) (iv2:interval) (iv3:interval) :=
forall a b, contained a iv1 -> contained b iv2 ->
contained (a + b) iv3.
Definition validIntervalSub (iv1:interval) (iv2:interval) (iv3:interval) :=
forall a b, contained a iv1 -> contained b iv2 ->
contained (a - b) iv3.
Definition validIntervalMult (iv1:interval) (iv2:interval) (iv3:interval) :=
forall a b, contained a iv1 -> contained b iv2 ->
contained (a * b) iv3.
Definition validIntervalDiv (iv1:interval) (iv2:interval) (iv3:interval) :=
forall a b, contained a iv1 -> contained b iv2 ->
contained (a / b) iv3.
(**
Now comes the old part with the computational definitions.
Where possible within time, they are shown sound with respect to the definitions from before, where not, we leave this as proof obligation for daisy.
**)
(**
TODO: Refactor into a list manipulating function instead
**)
Definition min4 (a:R) (b:R) (c:R) (d:R) := Rmin a (Rmin b (Rmin c d)).
Definition max4 (a:R) (b:R) (c:R) (d:R) := Rmax a (Rmax b (Rmax c d)).
(**
Asbtract interval update function, parametric by actual operator applied.
**)
Definition absIntvUpd (op:R->R->R) (I1:interval) (I2:interval) :=
(
min4 (op (IVlo I1) (IVlo I2))
(op (IVlo I1) (IVhi I2))
(op (IVhi I1) (IVlo I2))
(op (IVhi I1) (IVhi I2)),
max4 (op (IVlo I1) (IVlo I2))
(op (IVlo I1) (IVhi I2))
(op (IVhi I1) (IVlo I2))
(op (IVhi I1) (IVhi I2))
).
Definition widenInterval (Iv:interval) (v:R) := mkInterval (IVlo Iv - v) (IVhi Iv + v).
Definition negateInterval (intv:interval) := mkInterval (- IVhi intv) (- IVlo intv).
Lemma negationIsValid (intv:interval) (a:R) :
contained a intv-> contained (- a) (negateInterval intv).
Proof.
unfold contained; destruct intv as [ivlo ivhi]; simpl; intros [ivlo_le_a a_le_ivhi].
split; apply Ropp_ge_le_contravar; apply Rle_ge; auto.
Qed.
Definition addInterval (iv1:interval) (iv2:interval) :=
absIntvUpd Rplus iv1 iv2.
Lemma additionIsValid iv1 iv2:
validIntervalAdd iv1 iv2 (addInterval iv1 iv2).
Proof.
unfold validIntervalAdd.
intros a b.
unfold contained, addInterval, absIntvUpd, min4, max4.
intros [lo_leq_a a_leq_hi] [lo_leq_b b_leq_hi].
simpl; split.
(** Lower Bound **)
- assert ( fst iv1 + fst iv2 <= a + b)%R as lower_bound by (apply Rplus_le_compat; auto).
apply (Rle_trans _ (fst iv1 + fst iv2) _); [apply Rmin_l | auto].
(** Upper Bound **)
- assert (a + b <= snd iv1 + snd iv2)%R as upper_bound by (apply Rplus_le_compat; auto).
apply (Rle_trans _ (snd iv1 + snd iv2) _); [ auto | ].
apply (Rle_trans _ (Rmax (fst iv1 + snd iv2) (Rmax (snd iv1 + fst iv2) (snd iv1 + snd iv2))) _);
[ | apply Rmax_r].
apply (Rle_trans _ (Rmax (snd iv1 + fst iv2) (snd iv1 + snd iv2)) _ ); [ | apply Rmax_r].
apply (Rle_trans _ (snd iv1 + snd iv2) _); [ apply Req_le; auto | apply Rmax_r].
Qed.
Definition substractInterval (I1:interval) (I2:interval) :=
addInterval I1 (negateInterval I2).
Corollary subtractionIsValid iv1 iv2:
validIntervalSub iv1 iv2 (substractInterval iv1 iv2).
Proof.
unfold substractInterval.
intros a b.
intros contained_1 contained_I2.
assert (a - b = a + (- b))%R as simpl_eq by (field_simplify; reflexivity).
rewrite simpl_eq.
apply additionIsValid; auto.
apply negationIsValid; auto.
Qed.
(*
FIXME: This needs to be beautified
*)
Definition interval_mult_err (I1:interval) (e1:err) (I2:intv) (e2:err) :=
addInterval
(addInterval
(addInterval
(absIntvUpd Rmult I1 I2)
(absIntvUpd Rmult I1 (e1,e1)))
(absIntvUpd Rmult I2 (e2,e2)))
(Rmult e1 e2, Rmult e1 e2).
Definition multInterval (iv1:interval) (iv2:interval) :=
absIntvUpd Rmult iv1 iv2.
Definition absOpMin (op:R -> R -> R) (I1:interval) (I2:interval) :=
min4 (op (IVlo I1) (IVlo I2))
(op (IVlo I1) (IVhi I2))
(op (IVhi I1) (IVlo I2))
(op (IVhi I1) (IVhi I2)).
Definition absOpMax (op:R -> R -> R) (I1:interval) (I2:interval) :=
max4 (op (IVlo I1) (IVlo I2))
(op (IVlo I1) (IVhi I2))
(op (IVhi I1) (IVlo I2))
(op (IVhi I1) (IVhi I2)).
Lemma min_gt_0 (iv1:interval) (iv2:interval) :
valid iv1 ->
valid iv2 ->
(0 <= IVlo iv1)%R ->
(0 <= IVlo iv2)%R ->
absOpMin (Rmult) iv1 iv2 = (IVlo iv1 * IVlo iv2)%R.
Proof.
intros valid_iv1 valid_iv2 geq0_1 geq0_2.
unfold absOpMin, min4, Rmin; destruct Rle_dec; auto.
exfalso.
apply n.
destruct Rle_dec.
- apply Rmult_le_compat; auto.
apply Req_le; auto.
- destruct Rle_dec; apply Rmult_le_compat; auto.
apply Req_le; auto.
Qed.
Lemma max_gt_0 (iv1:interval) (iv2:interval) :
valid iv1 ->
valid iv2 ->
(0 <= IVlo iv1)%R ->
(0 <= IVlo iv2)%R ->
absOpMax (Rmult) iv1 iv2 = (IVhi iv1 * IVhi iv2)%R.
Proof.
intros valid_iv1 valid_iv2 geq0_1 geq0_2.
unfold absOpMax, max4, Rmax; destruct Rle_dec.
- destruct Rle_dec.
+ destruct Rle_dec; auto.
exfalso.
apply n.
apply Rmult_le_compat; auto.
* apply (Rle_trans _ (IVlo iv1) _); auto.
* apply Req_le; auto.
+ exfalso; apply n.
destruct Rle_dec.
* apply Rmult_le_compat; auto; [apply (Rle_trans _ (IVlo iv2) _); auto | apply Req_le; auto].
* exfalso; apply n0; apply Rmult_le_compat; auto; [ apply (Rle_trans _ (IVlo iv1) _); auto | apply Req_le; auto].
- exfalso.
apply n.
destruct Rle_dec.
+ destruct Rle_dec.
* apply Rmult_le_compat; auto.
* apply Rmult_le_compat; auto.
apply Req_le; auto.
+ apply Rmult_le_compat; auto.
apply Req_le; auto.
Qed.
Definition intv_div_err (I1:interval) (e1:err) (I2:interval) (e2:err) :=
interval_mult_err I1 e1 ( (/ IVlo (I2))%R, (/ IVhi (I2))%R) e2.
(**
Useful Lemmata:
Rinv_mult_simpl:
forall r1 r2 r3 : R,
r1 <> 0%R -> (r1 * / r2 * (r3 * / r1))%R = (r3 * / r2)%R
Rmult_lt_compat_r:
forall r r1 r2 : R, (0 < r)%R -> (r1 < r2)%R -> (r1 * r < r2 * r)%R
Rmult_gt_compat_r:
forall r r1 r2 : R, (r > 0)%R -> (r1 > r2)%R -> (r1 * r > r2 * r)%R
Rmult_gt_compat_l:
forall r r1 r2 : R, (r > 0)%R -> (r1 > r2)%R -> (r * r1 > r * r2)%R
Rmult_le_compat_l:
forall r r1 r2 : R, (0 <= r)%R -> (r1 <= r2)%R -> (r * r1 <= r * r2)%R
Rmult_le_compat_r:
forall r r1 r2 : R, (0 <= r)%R -> (r1 <= r2)%R -> (r1 * r <= r2 * r)%R
Rmult_ge_compat_l:
forall r r1 r2 : R, (r >= 0)%R -> (r1 >= r2)%R -> (r * r1 >= r * r2)%R
Rmult_ge_compat_r:
forall r r1 r2 : R, (r >= 0)%R -> (r1 >= r2)%R -> (r1 * r >= r2 * r)%R
Rmult_le_compat:
forall r1 r2 r3 r4 : R,
(0 <= r1)%R ->
(0 <= r3)%R -> (r1 <= r2)%R -> (r3 <= r4)%R -> (r1 * r3 <= r2 * r4)%R
Rmult_ge_compat:
forall r1 r2 r3 r4 : R,
(r2 >= 0)%R ->
(r4 >= 0)%R -> (r1 >= r2)%R -> (r3 >= r4)%R -> (r1 * r3 >= r2 * r4)%R
Rmult_gt_0_lt_compat:
forall r1 r2 r3 r4 : R,
(r3 > 0)%R ->
(r2 > 0)%R -> (r1 < r2)%R -> (r3 < r4)%R -> (r1 * r3 < r2 * r4)%R
Rmult_le_0_lt_compat:
forall r1 r2 r3 r4 : R,
(0 <= r1)%R ->
(0 <= r3)%R -> (r1 < r2)%R -> (r3 < r4)%R -> (r1 * r3 < r2 * r4)%R
Rmult_lt_0_compat:
forall r1 r2 : R, (0 < r1)%R -> (0 < r2)%R -> (0 < r1 * r2)%R
Rmult_gt_0_compat:
forall r1 r2 : R, (r1 > 0)%R -> (r2 > 0)%R -> (r1 * r2 > 0)%R
Rmult_le_compat_neg_l:
forall r r1 r2 : R, (r <= 0)%R -> (r1 <= r2)%R -> (r * r2 <= r * r1)%R
Rmult_le_ge_compat_neg_l:
forall r r1 r2 : R, (r <= 0)%R -> (r1 <= r2)%R -> (r * r1 >= r * r2)%R
Rmult_lt_gt_compat_neg_l:
forall r r1 r2 : R, (r < 0)%R -> (r1 < r2)%R -> (r * r1 > r * r2)%R
**)
(**
Lemma IntervalMultMono (I1:interval) (I2:interval) (a:R) (b:R) :
contained I1 a ->
contained I2 b ->
contained (interval_mult I1 I2) (a * b).
Proof.
unfold contained, interval_mult, absIntervalUpd, IVlo, IVhi;
destruct I1 as [ivlo1 ivhi1]; destruct I2 as [ivlo2 ivhi2]; simpl.
intros [ lo_leq_a a_leq_hi] [lo_leq_b b_leq_hi].
assert (ivlo1 <= ivhi1)%R by (apply (Rle_trans _ a _); auto).
assert (ivlo2 <= ivhi2)%R by (apply (Rle_trans _ b _); auto).
destruct (Rle_dec ivhi1 0); destruct (Rle_dec ivhi2 0);
[(* Both less do no further case distinction *)
assert (ivlo1 <= 0)%R by (apply (Rle_trans _ (ivhi1) _); auto);
assert (ivlo2 <= 0)%R by (apply (Rle_trans _ (ivhi2) _); auto)
| (* Fst one less, destruct second one *)
assert (ivlo1 <= 0)%R by (apply (Rle_trans _ (ivhi1) _); auto);
destruct (Rle_dec ivlo2 0)
| (* Opposite case *)
assert (ivlo2 <= 0)%R by (apply (Rle_trans _ (ivhi2) _); auto);
destruct (Rle_dec ivlo1 0)
| (* Both greater -> both destruct *)
destruct (Rle_dec ivlo2 0); destruct (Rle_dec ivlo1 0) ].
-
(*
First subcase, all intervals negative --> ivhi1 * ivhi2 is maximum
Rmult_le_compat_neg_l b a ivhi1
Rmult_le_compat_neg_l ivlo1 b ivhi2
*)
assert (b <= 0)%R as b_leq_0 by (apply (Rle_trans _ (ivhi2) _); auto).
split.
+ pose proof (Rmult_le_compat_neg_l b a ivhi1 b_leq_0 a_leq_hi).
pose proof (Rmult_le_compat_neg_l ivhi1 b ivhi2 r b_leq_hi).
rewrite Rmult_comm in H3.
pose proof (Rle_trans _ _ _ H4 H3).
rewrite (Rmult_comm a b).
unfold min4.
assert
(Rmin (ivlo1 * ivlo2)
(Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) <= ivhi1 * ivhi2)%R.
* assert (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) <= ivhi1 * ivhi2)%R.
{
assert (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2) <= ivhi1 * ivhi2)%R by apply Rmin_r.
apply (Rle_trans _ (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) _); auto.
apply Rmin_r.
}
{
apply (Rle_trans _ (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) _); auto.
apply Rmin_r.
}
* apply (Rle_trans _ (ivhi1 * ivhi2) _); auto.
+ unfold max4.
apply Rge_le.
pose proof (Rmult_le_ge_compat_neg_l ivlo1 ivlo2 b H1 lo_leq_b).
pose proof (Rmult_le_ge_compat_neg_l b ivlo1 a b_leq_0 lo_leq_a).
rewrite Rmult_comm in H4.
pose proof (Rge_trans _ _ _ H3 H4).
rewrite (Rmult_comm a b).
apply (Rge_trans _ (ivlo1 * ivlo2) _); auto.
apply Rle_ge, Rmax_l.
- assert (a <= 0)%R as a_leq_0 by (apply (Rle_trans _ (ivhi1) _); auto).
split.
+ apply Rnot_le_gt in n.
apply Rgt_ge in n.
pose proof (Rmult_le_compat_neg_l a b ivhi2).
pose proof (Rmult_le_compat_neg_l ivhi1).
rewrite Rmult_comm in H3.
pose proof (Rle_trans _ _ _ H4 H3).
rewrite (Rmult_comm a b).
unfold min4.
assert
(Rmin (ivlo1 * ivlo2)
(Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) <= ivhi1 * ivhi2)%R.
* assert (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) <= ivhi1 * ivhi2)%R.
{
assert (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2) <= ivhi1 * ivhi2)%R by apply Rmin_r.
apply (Rle_trans _ (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2)) _); auto.
apply Rmin_r.
}
{
apply (Rle_trans _ (Rmin (ivlo1 * ivhi2) (Rmin (ivhi1 * ivlo2) (ivhi1 * ivhi2))) _); auto.
apply Rmin_r.
}
* apply (Rle_trans _ (ivhi1 * ivhi2) _); auto.
+ unfold max4.
apply Rge_le.
pose proof (Rmult_le_ge_compat_neg_l ivlo1 ivlo2 b H1 lo_leq_b).
pose proof (Rmult_le_ge_compat_neg_l b ivlo1 a b_leq_0 lo_leq_a).
rewrite Rmult_comm in H4.
pose proof (Rge_trans _ _ _ H3 H4).
rewrite (Rmult_comm a b).
apply (Rge_trans _ (ivlo1 * ivlo2) _); auto.
apply Rle_ge, Rmax_l.
- admit.
- admit.
- admit.
- admit.
- admit.
- admit.
- admit.
**)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment