ErrorValidation.v 10.3 KB
Newer Older
Heiko Becker's avatar
Heiko Becker committed
1
Require Import Coq.QArith.QArith Coq.QArith.Qminmax Coq.QArith.Qabs Coq.QArith.Qreals.
2
(*Coq.QArith.Qcanon.*)
Heiko Becker's avatar
Heiko Becker committed
3 4 5
Require Import Coq.micromega.Psatz Coq.Reals.Reals .
Require Import Daisy.Infra.Abbrevs Daisy.Infra.RationalSimps Daisy.Infra.RealRationalProps.
Require Import Daisy.Expressions Daisy.Infra.RationalConstruction Daisy.Infra.RealSimps Daisy.IntervalArith Daisy.ErrorBounds.
6 7 8 9 10

Section ComputableErrors.

  Definition machineEpsilon := (1#(2^53)).

Heiko Becker's avatar
Heiko Becker committed
11 12 13 14 15 16 17 18
  Fixpoint toRExp (e:exp Q) :=
    match e with
    |Var v => Var R v
    |Param v => Param R v
    |Const n => Const (Q2R n)
    |Binop b e1 e2 => Binop b (toRExp e1) (toRExp e2)
    end.

19
  Fixpoint validErrorbound (e:exp Q) (env:analysisResult) :=
20 21 22 23
    let (intv, err) := (env e) in
      match e with
        |Var v => false
        |Param v => Qleb (maxAbs intv * machineEpsilon) err
Heiko Becker's avatar
Heiko Becker committed
24
        (* A constant will be a point intv. Take maxAbs never the less *)
25 26 27 28 29 30 31
        |Const n => Qleb (maxAbs intv * machineEpsilon) err
        |Binop b e1 e2 =>
         let (ive1, err1) := env e1 in
         let (ive2, err2) := env e2 in
         let rec := andb (validErrorbound e1 env) (validErrorbound e2 env) in
         let upperBoundE1 := maxAbs ive1 in
         let upperBoundE2 := maxAbs ive2 in
Heiko Becker's avatar
Heiko Becker committed
32 33 34 35 36
         (* WAS:
         let e1F := upperBoundE1 + (upperBoundE1 * err1) in
         let e2F := upperBoundE2 + (upperBoundE2 * err2) in *)
         let e1F := upperBoundE1 + err1 in
         let e2F := upperBoundE2 + err2 in
37 38
         let theVal :=
            match b with
Heiko Becker's avatar
Heiko Becker committed
39 40 41 42
            | Plus => Qleb (err1 + err2 + (Qabs e1F + Qabs e2F) * machineEpsilon) err
            (* TODO:Validity of next two computations *)
            | Sub => Qleb (err1 + err2 + (Qabs (e1F - e2F) * machineEpsilon)) err
            | Mult => Qleb ((upperBoundE1 * upperBoundE2) + (e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2) + (((e1F * e2F) + (e1F * e2F * err1) + (e1F * e2F * err2) + (e1F * e2F * err1 * err2))*machineEpsilon)) err
43 44 45 46
            | Div => false
            end
              in andb rec theVal
      end.
47

Heiko Becker's avatar
Heiko Becker committed
48 49
  Functional Scheme validErrorbound_ind := Induction for validErrorbound Sort Prop.

50 51 52 53 54 55 56 57 58 59
  Definition u:nat := 1.
  (** 1655/5 = 331; 0,4 = 2/5 **)
  Definition cst1:Q := 1657 # 5.

(** Define abbreviations **)
  Definition varU:exp Q := Param Q u.
  Definition valCst:exp Q := Const cst1.
  Definition valCstAddVarU:exp Q := Binop Plus valCst varU.

  (** Error values **)
60 61 62 63 64
  Definition errCst1 :Q := rationalFromNum (36792791036077685#1) 16 14.
  Definition errVaru := rationalFromNum (11102230246251565#1) 16 14.
  Definition lowerBoundAddUCst:Q := 1157 # 5.
  Definition upperBoundAddUCst:Q := 2157 # 5.
  Definition errAddUCst := rationalFromNum (9579004256465851#1) 15 14.
65 66 67 68 69 70

  (** The added assertion becomes the precondition for us **)
  Definition precondition := fun env:nat -> R => (- 100 <= env u)%R /\ (env u <= 100)%R.

  (** As stated, we need to define the absolute environment now as an inductive predicate
  Inductive absEnv : abs_env :=
Heiko Becker's avatar
Heiko Becker committed
71 72 73
    theCst: absEnv valCst (mkIntv cst1 cst1) errCst1
   |theVar: absEnv varU (mkIntv (- 100) (100)) errVaru
   |theAddition: absEnv valCstAddVarU (mkIntv lowerBoundAddUCst upperBoundAddUCst) errAddUCst. **)
74

75 76 77 78 79 80 81 82
  Definition validConstant := Eval compute in validErrorbound (valCst) (fun x => ((cst1,cst1),errCst1)).
  Definition validParam := Eval  compute in validErrorbound (varU) (fun x => (((-100) # 1,100 # 1),errVaru)).
  Definition validAddition := Eval compute in validErrorbound (valCstAddVarU)
                                                              (fun x => match x with
                                                                     |Param _ => ((cst1,cst1),errCst1)
                                                                     |Const _ => (((-100) # 1,100 # 1),errVaru)
                                                                     |_ => ((lowerBoundAddUCst,upperBoundAddUCst),errAddUCst) end).

Heiko Becker's avatar
Heiko Becker committed
83
  Definition envApproximatesPrecond (P:nat -> intv) (absenv:analysisResult) :=
84 85 86 87 88 89
    forall u:nat,
      let (ivlo,ivhi) := P u in
      let (absIv,err) := absenv (Param Q u) in
      let (abslo,abshi) := absIv in
      (abslo <= ivlo /\ ivhi <= abshi)%Q.

Heiko Becker's avatar
Heiko Becker committed
90
  Definition precondValidForExec (P:nat->intv) (cenv:nat->R) :=
91 92 93
    forall v:nat,
      let (ivlo,ivhi) := P v in
      (Q2R ivlo <= cenv v <= Q2R ivhi)%R.
94 95

  Lemma validErrorboundCorrectConstant:
Heiko Becker's avatar
Heiko Becker committed
96
    forall cenv absenv (n:Q) nR nF e nlo nhi,
97 98 99
      eval_exp 0%R cenv (Const (Q2R n)) nR ->
      eval_exp (Q2R machineEpsilon) cenv (Const (Q2R n)) nF ->
      validErrorbound (Const n) absenv = true ->
Heiko Becker's avatar
Heiko Becker committed
100
      absenv (Const n) = ((nlo,nhi),e) ->
101 102
      (Rabs (nR - nF) <= (Q2R e))%R.
  Proof.
Heiko Becker's avatar
Heiko Becker committed
103
    intros cenv absenv n nR nF e nlo nhi eval_real eval_float error_valid absenv_const.
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
    unfold validErrorbound in error_valid.
    rewrite absenv_const in error_valid.
    inversion eval_real; subst.
    rewrite perturb_0_val in eval_real; auto.
    rewrite perturb_0_val; auto.
    clear delta H0.
    inversion eval_float; subst.
    unfold perturb in *; simpl in *.
    rewrite Rabs_err_simpl.
    unfold Qleb in error_valid.
    apply Qle_bool_iff in error_valid.
    apply Qle_Rle in error_valid.
    eapply Rle_trans.
    rewrite Rabs_mult.
    eapply Rmult_le_compat_l; [ apply Rabs_pos | ].
    apply H0.
    rewrite Rabs_eq_Qabs.
    rewrite Q2R_mult in error_valid.
Heiko Becker's avatar
Heiko Becker committed
122 123 124
    eapply Rle_trans.
    eapply Rmult_le_compat_r.
    (*
125
    apply error_valid.
Heiko Becker's avatar
Heiko Becker committed
126 127
  Qed. *)
  Admitted.
128

129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
  Lemma validErrorboundCorrectParam:
    forall cenv absenv (v:nat) nR nF e P plo phi,
      envApproximatesPrecond P absenv ->
      precondValidForExec P cenv ->
      eval_exp 0%R cenv (Param R v) nR ->
      eval_exp (Q2R machineEpsilon) cenv (Param R v) nF ->
      validErrorbound (Param Q v) absenv = true ->
      absenv (Param Q v)  = ((plo,phi),e) ->
      (Rabs (nR - nF) <= (Q2R e))%R.
  Proof.
    intros cenv absenv v nR nF e P plo phi absenv_approx_p cenv_approx_p eval_real eval_float error_valid absenv_param.
    unfold validErrorbound in error_valid.
    rewrite absenv_param in error_valid.
    inversion eval_real; subst.
    rewrite perturb_0_val in eval_real; auto.
    rewrite perturb_0_val; auto.
    inversion eval_float; subst.
    unfold envApproximatesPrecond in absenv_approx_p.
    unfold precondValidForExec in cenv_approx_p.
    specialize (absenv_approx_p v).
    specialize (cenv_approx_p v).
    assert (exists ivlo ivhi,  (ivlo,ivhi) = P v) by (destruct (P v) as [ivlo ivhi]; repeat eexists).
    destruct H as [ivlo [ivhi P_eq]].
    rewrite <- P_eq in absenv_approx_p, cenv_approx_p.
    rewrite absenv_param in absenv_approx_p.
    destruct absenv_approx_p as [plo_le_ivlo ivhi_le_phi].
    destruct cenv_approx_p as [ivlo_le_cenv cenv_le_ivhi].
    apply Qle_Rle in plo_le_ivlo; apply Qle_Rle in ivhi_le_phi.
    pose proof (Rle_trans (Q2R plo) (Q2R ivlo) (cenv v) plo_le_ivlo ivlo_le_cenv).
    pose proof (Rle_trans (cenv v) (Q2R ivhi) (Q2R phi) cenv_le_ivhi ivhi_le_phi).
    pose proof (RmaxAbs (Q2R plo) (cenv v) (Q2R phi) H H2).
    unfold perturb.
    rewrite Rabs_err_simpl.
    rewrite Rabs_mult.
    eapply Rle_trans.
    eapply Rmult_le_compat; [ apply Rabs_pos | apply Rabs_pos |  | ].
    (* TODO: Does that work out? *)
    apply H3.
    apply H1.
    pose proof (maxAbs_impl_RmaxAbs (plo,phi)).
    unfold RmaxAbsFun in H4.
    simpl in H4.
    rewrite H4.
    apply Qle_bool_iff in error_valid.
    apply Qle_Rle in error_valid.
    rewrite Q2R_mult in error_valid.
    apply error_valid.
  Qed.

Heiko Becker's avatar
Heiko Becker committed
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
  Lemma validErrorboundCorrectAddition cenv absenv (e1:exp Q) (e2:exp Q) (nR nR1 nR2 nF nF1 nF2 :R) (e err1 err2 :error) (alo ahi e1lo e1hi e2lo e2hi:Q) :
      eval_exp 0%R cenv (toRExp e1) nR1 ->
      eval_exp 0%R cenv (toRExp e2) nR2 ->
      eval_exp 0%R cenv (toRExp (Binop Plus e1 e2)) nR ->
      eval_exp (Q2R machineEpsilon) cenv (toRExp e1) nF1 ->
      eval_exp (Q2R machineEpsilon) cenv (toRExp e2) nF2 ->
      eval_exp (Q2R machineEpsilon) (updEnv 1 nF1 (updEnv 2 nF2 cenv)) (toRExp (Binop Plus (Var Q 1) (Var Q 2))) nF ->
      validErrorbound (Binop Plus e1 e2) absenv = true ->
      absenv e1 = ((e1lo,e1hi),err1) ->
      absenv e2 = ((e2lo, e2hi),err2) ->
      absenv (Binop Plus e1 e2) = ((alo,ahi),e)->
      (Rabs (nR1 - nF1) <= (Q2R err1))%R ->
      (Rabs (nR2 - nF2) <= (Q2R err2))%R ->
      (Rabs (nR - nF) <= (Q2R e))%R.
  Proof.
    intros e1_real e2_real eval_real e1_float e2_float eval_float valid_error absenv_e1 absenv_e2 absenv_add err1_bounded err2_bounded.
    eapply Rle_trans.
    eapply add_abs_err_bounded.
    apply e1_real.
    admit.
    apply e2_real.
    admit.
    (* TODO: Fix machine epsilons *)
    apply eval_real.
    admit.
    apply err1_bounded.
    apply err2_bounded.
    unfold validErrorbound in valid_error at 1.
    rewrite absenv_add, absenv_e1, absenv_e2 in valid_error.
    apply Is_true_eq_left in valid_error.
    apply andb_prop_elim in valid_error.
    destruct valid_error as [ valid_rec valid_error].
    apply andb_prop_elim in valid_rec.
    destruct valid_rec as [valid_e1 valid_e2].
    apply Is_true_eq_true in valid_error.
    apply Qle_bool_iff in valid_error.
    apply Qle_Rle in valid_error.
    repeat rewrite Q2R_plus in valid_error.
    repeat rewrite Q2R_mult in valid_error.
    repeat rewrite Q2R_plus in valid_error.
    repeat rewrite Q2R_mult in valid_error.
    repeat rewrite <- Rabs_eq_Qabs in valid_error.
    repeat rewrite Q2R_plus in valid_error.
    repeat rewrite Q2R_mult in valid_error.
    repeat rewrite <- maxAbs_impl_RmaxAbs in valid_error.
    eapply Rle_trans.
    apply Rplus_le_compat_l.
    eapply Rmult_le_compat. (* FIXME: Make compat_r when machineEpsilons are fixed *)
    apply Rabs_pos.
    admit.
    Focus 3.
    apply valid_error.
    Focus 2.
    admit.
    eapply Rle_trans.
    eapply Rabs_triang.
    eapply Rplus_le_compat.
Heiko Becker's avatar
Heiko Becker committed
235
  Admitted.
Heiko Becker's avatar
Heiko Becker committed
236

Heiko Becker's avatar
Heiko Becker committed
237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252
  Lemma validErrorboundCorrect (e:exp Q):
    forall cenv absenv nR nF err P elo ehi,
      envApproximatesPrecond P absenv ->
      precondValidForExec P cenv ->
      eval_exp 0%R cenv (toRExp e) nR ->
      eval_exp (Q2R machineEpsilon) cenv (toRExp e) nF ->
      validErrorbound e absenv = true ->
      absenv e = ((elo,ehi),err) ->
      (Rabs (nR - nF) <= (Q2R err))%R.
  Proof.
    induction e.
    - intros; simpl in *.
      rewrite H4 in H3; inversion H3.
    - intros; eapply validErrorboundCorrectParam; eauto.
    - intros; eapply validErrorboundCorrectConstant; eauto.
  Admitted.
253 254

End ComputableErrors.